[CFD-Lover]

Hi,

I have been reading the source code of Gmsh and got some questions that I wasn't able to understand the meaning of some of the parameters. Here is the source code;

typedef struct {
  int Num;
  // t is the local coordinate of the point
  // lc is x'(t)/h(x(t))
  // p is the value of the primitive
  // xp is the norm of the tangent vector
  double t, lc, p, xp;
} IntPoint;

static double smoothPrimitive(GEdge *ge, double alpha,
                              std::vector<IntPoint> &Points)
{
  int ITER = 0;
  while (1){
    int count1 = 0;
    int count2 = 0;

    // use a gauss-seidel iteration; iterate forward and then backward;
    // convergence is usually very fast
    for (unsigned int i = 1; i < Points.size(); i++){
      double dh = (Points[i].xp/Points[i].lc - Points[i-1].xp/Points[i-1].lc);
      double dt = Points[i].t - Points[i-1].t;
      double dhdt =  dh/dt;
      if (dhdt / Points[i].xp > (alpha - 1.)*1.01){
    double hnew = Points[i-1].xp / Points[i-1].lc + dt * (alpha-1) * Points[i].xp;
    Points[i].lc = Points[i].xp / hnew;
    count1++;
      }
    }

    for (int i = Points.size() - 1; i > 0; i--){
      double dh = (Points[i-1].xp/Points[i-1].lc - Points[i].xp/Points[i].lc);
      double dt = fabs(Points[i-1].t - Points[i].t);
      double dhdt =  dh/dt;
      if (dhdt / Points[i-1].xp > (alpha-1.)*1.01){
          double hnew = Points[i].xp / Points[i].lc + dt * (alpha-1) * Points[i].xp;
          Points[i-1].lc = Points[i-1].xp / hnew;
          count2 ++;
      }
    }

    ++ITER;
    if (ITER > 2000) break;
    if (!(count1 + count2)) break;
  }

  // recompute the primitive
  for (int i = 1; i < (int)Points.size(); i++){
    IntPoint &pt2 = Points[i];
    IntPoint &pt1 = Points[i-1];
    pt2.p = pt1.p + (pt2.t - pt1.t) * 0.5 * (pt2.lc + pt1.lc);
  }
  return Points[Points.size() - 1].p;
}

static double F_Lc(GEdge *ge, double t)
{
  GPoint p = ge->point(t);
  double lc_here;

  Range<double> bounds = ge->parBounds(0);
  double t_begin = bounds.low();
  double t_end = bounds.high();

  if(t == t_begin)
    lc_here = BGM_MeshSize(ge->getBeginVertex(), t, 0, p.x(), p.y(), p.z());
  else if(t == t_end)
    lc_here = BGM_MeshSize(ge->getEndVertex(), t, 0, p.x(), p.y(), p.z());
  else
    lc_here = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());

  SVector3 der = ge->firstDer(t);
  const double d = norm(der);

  return d / lc_here;
}

static double F_Lc_aniso(GEdge *ge, double t)
{
#if defined(HAVE_ANN)
  FieldManager *fields = ge->model()->getFields();
  BoundaryLayerField *blf = 0;
  Field *bl_field = fields->get(fields->getBoundaryLayerField());
  blf = dynamic_cast<BoundaryLayerField*> (bl_field);
#else
  bool blf = false;
#endif


  GPoint p = ge->point(t);
  SMetric3 lc_here;

  Range<double> bounds = ge->parBounds(0);
  double t_begin = bounds.low();
  double t_end = bounds.high();

  if(t == t_begin)
    lc_here = BGM_MeshMetric(ge->getBeginVertex(), t, 0, p.x(), p.y(), p.z());
  else if(t == t_end)
    lc_here = BGM_MeshMetric(ge->getEndVertex(), t, 0, p.x(), p.y(), p.z());
  else
    lc_here = BGM_MeshMetric(ge, t, 0, p.x(), p.y(), p.z());

#if defined(HAVE_ANN)
  if (blf && !blf->isEdgeBL(ge->tag())){
    SMetric3 lc_bgm;
    blf->computeFor1dMesh ( p.x(), p.y(), p.z() , lc_bgm );
    lc_here = intersection_conserveM1 (lc_here, lc_bgm );
  }
#endif

  SVector3 der = ge->firstDer(t);
  double lSquared = dot(der, lc_here, der);
  return sqrt(lSquared);
}

static double F_Transfinite(GEdge *ge, double t_)
{
  double length = ge->length();
  if(length == 0.0){
    Msg::Error("Zero-length curve %d in transfinite mesh", ge->tag());
    return 1.;
  }

  SVector3 der = ge->firstDer(t_) ;
  double d = norm(der);
  double coef = ge->meshAttributes.coeffTransfinite;
  int type = ge->meshAttributes.typeTransfinite;
  int nbpt = ge->meshAttributes.nbPointsTransfinite;

  if(CTX::instance()->mesh.flexibleTransfinite && CTX::instance()->mesh.lcFactor)
    nbpt /= CTX::instance()->mesh.lcFactor;

  Range<double> bounds = ge->parBounds(0);
  double t_begin = bounds.low();
  double t_end = bounds.high();
  double t = (t_ - t_begin)/(t_end-t_begin);

  double val;

  if(coef <= 0.0 || coef == 1.0) {
    // coef < 0 should never happen
    val = d * coef / ge->length();
  }
  else {
    switch (std::abs(type)) {

It says in the code that the parameter "t" is the local coordinate of the point. I don't what does it mean if I had a 2x2 domain with 4x4 nodes. Same with the parameter d, which is double d = norm(der);. Furthermroe, it says that coef = ge->meshAttributes.coeffTransfinite;, so what meshAttributes.coeffTransfinite means. Any idea what do these parameters mean?

Thanks for the help,