|
[Sponsors] |
July 28, 2017, 13:10 |
Gmsh source code-Reading issue!
|
#1 |
Member
OpenFoam
Join Date: Jun 2016
Posts: 82
Rep Power: 10 |
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)) { Thanks for the help, |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
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 |