CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Gmsh source code-Reading issue!

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   July 28, 2017, 13:10
Unhappy Gmsh source code-Reading issue!
  #1
Member
 
OpenFoam
Join Date: Jun 2016
Posts: 82
Rep Power: 9
CFD-Lover is on a distinguished road
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;

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,
CFD-Lover is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Custom Thermophysical Properties wsmith02 OpenFOAM 4 June 1, 2023 14:30
[Other] Adding solvers from DensityBasedTurbo to foam-extend 3.0 Seroga OpenFOAM Community Contributions 9 June 12, 2015 17:18
Problem compiling a custom Lagrangian library brbbhatti OpenFOAM Programming & Development 2 July 7, 2014 11:32
"parabolicVelocity" in OpenFoam 2.1.0 ? sawyer86 OpenFOAM Running, Solving & CFD 21 February 7, 2012 11:44
DecomposePar links against liblamso0 with OpenMPI jens_klostermann OpenFOAM Bugs 11 June 28, 2007 17:51


All times are GMT -4. The time now is 23:34.