Skip to content
Snippets Groups Projects
Select Git revision
  • gmsh_4_10_4
  • master default protected
  • patch_releases_4_14
  • overlaps_tags_and_distributed_export
  • overlaps_tags_and_distributed_export_rebased
  • relaying
  • alphashapes
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
40 results

t19.geo

Blame
  • meshGEdge.cpp 21.93 KiB
    // Gmsh - Copyright (C) 1997-2016 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to the public mailing list <gmsh@onelab.info>.
    
    #include "GmshConfig.h"
    #include "GModel.h"
    #include "discreteEdge.h"
    #include "meshGEdge.h"
    #include "GEdge.h"
    #include "MLine.h"
    #include "BackgroundMeshTools.h"
    #include "boundaryLayersData.h"
    #include "Numeric.h"
    #include "GmshMessage.h"
    #include "Context.h"
    #include "STensor3.h"
    #include "Field.h"
    #include "OS.h"
    
    #define SQU(a)      ((a)*(a))
    
    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_LcB(GEdge *ge, double t)
    {
      /*  BoundaryLayerField *blf = 0;
    #if defined(HAVE_ANN)
      FieldManager *fields = ge->model()->getFields();
      Field *bl_field = fields->get(fields->getBoundaryLayerField());
      blf = dynamic_cast<BoundaryLayerField*> (bl_field);
    #endif
      */
      GPoint p = ge->point(t);
      double lc = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
    
      /*  if (blf){
        double lc2 = (*blf)( p.x(), p.y(), p.z() , ge);
        //    printf("p %g %g lc %g\n",p.x(),p.y(),lc2);
        lc = std::min(lc, lc2);
      }
      */
    
      return lc;
    }
    
    static double F_Lc(GEdge *ge, double t)
    {
      /*
      BoundaryLayerField *blf = 0;
    #if defined(HAVE_ANN)
      FieldManager *fields = ge->model()->getFields();
      Field *bl_field = fields->get(fields->getBoundaryLayerField());
      blf = dynamic_cast<BoundaryLayerField*> (bl_field);
    #endif
      */
      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());
      /*
      if (blf){
        double lc2 = (*blf)( p.x(), p.y(), p.z() , ge);
        //    printf("p %g %g lc %g\n",p.x(),p.y(),lc2);
        lc_here = std::min(lc_here, lc2);
      }
      */
      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
    
      //printf("coucou\n");
    
    
      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)) {
    
        case 1: // Geometric progression ar^i; Sum of n terms = length = a (r^n-1)/(r-1)
          {
            double r = (gmsh_sign(type) >= 0) ? coef : 1. / coef;
            double a = length * (r - 1.) / (pow(r, nbpt - 1.) - 1.);
            int i = (int)(log(t * length / a * (r - 1.) + 1.) / log(r));
            val = d / (a * pow(r, (double)i));
          }
          break;
    
        case 2: // Bump
          {
            double a;
            if(coef > 1.0) {
              a = -4. * sqrt(coef - 1.) * atan2(1., sqrt(coef - 1.)) /
                ((double)nbpt *  length);
            }
            else {
              a = 2. * sqrt(1. - coef) * log(fabs((1. + 1. / sqrt(1. - coef)) /
                                                  (1. - 1. / sqrt(1. - coef))))
                / ((double)nbpt * length);
            }
            double b = -a * length * length / (4. * (coef - 1.));
            val = d / (-a * SQU(t * length - (length) * 0.5) + b);
          }
          break;
    
        default:
          Msg::Warning("Unknown case in Transfinite Line mesh");
          val = 1.;
          break;
        }
      }
      return val;
    }
    
    static double F_One(GEdge *ge, double t)
    {
      SVector3 der = ge->firstDer(t) ;
      return norm(der);
    }
    
    static double trapezoidal(IntPoint * P1, IntPoint * P2)
    {
      return (0.5 * (P1->lc + P2->lc) * (P2->t - P1->t));
    }
    
    static void RecursiveIntegration(GEdge *ge, IntPoint *from, IntPoint *to,
                                     double (*f) (GEdge *e, double X),
                                     std::vector<IntPoint> &Points,
                                     double Prec, int *depth)
    {
      IntPoint P, p1;
    
      (*depth)++;
    
      P.t = 0.5 * (from->t + to->t);
      P.lc = f(ge, P.t);
    
      double val1 = trapezoidal(from, to);
      double val2 = trapezoidal(from, &P);
      double val3 = trapezoidal(&P, to);
      double err = fabs(val1 - val2 - val3);
    
      if(((err < Prec) && (*depth > 6)) || (*depth > 25)) {
        p1 = Points.back();
        P.p = p1.p + val2;
        Points.push_back(P);
    
        p1 = Points.back();
        to->p = p1.p + val3;
        Points.push_back(*to);
      }
      else {
        RecursiveIntegration(ge, from, &P, f, Points, Prec, depth);
        RecursiveIntegration(ge, &P, to, f, Points, Prec, depth);
      }
    
      (*depth)--;
    }
    
    static double Integration(GEdge *ge, double t1, double t2,
                              double (*f) (GEdge *e, double X),
                              std::vector<IntPoint> &Points, double Prec)
    {
      IntPoint from, to;
    
      int depth = 0;
    
      from.t = t1;
      from.lc = f(ge, from.t);
      from.p = 0.0;
      Points.push_back(from);
    
      to.t = t2;
      to.lc = f(ge, to.t);
    
      RecursiveIntegration(ge, &from, &to, f, Points, Prec, &depth);
    
      return Points.back().p;
    }
    
    void copyMesh(GEdge *from, GEdge *to, int direction)
    {
      Range<double> u_bounds = from->parBounds(0);
      double u_min = u_bounds.low();
      double u_max = u_bounds.high();
    
      Range<double> to_u_bounds = to->parBounds(0);
      double to_u_min = to_u_bounds.low();
    
      // include begin and end point to avoid conflicts when realigning
    
      MVertex* vt0 = to->getBeginVertex()->mesh_vertices[0];
      MVertex* vt1 = to->getEndVertex()->mesh_vertices[0];
    
      MVertex* vs0 = from->getBeginVertex()->mesh_vertices[0];
      MVertex* vs1 = from->getEndVertex()->mesh_vertices[0];
    
      to->correspondingVertices[vt0] = direction > 0 ? vs0 : vs1;
      to->correspondingVertices[vt1] = direction > 0 ? vs1 : vs0;
    
      for(unsigned int i = 0; i < from->mesh_vertices.size(); i++){
        int index = (direction < 0) ? (from->mesh_vertices.size() - 1 - i) : i;
        MVertex *v = from->mesh_vertices[index];
        double u; v->getParameter(0, u);
        double newu = (direction > 0) ? (u-u_min+to_u_min) : (u_max-u+to_u_min);
        GPoint gp = to->point(newu);
        MEdgeVertex *vv = new MEdgeVertex(gp.x(), gp.y(), gp.z(), to, newu);
        to->mesh_vertices.push_back(vv);
        to->correspondingVertices[vv] = v;
      }
      for(unsigned int i = 0; i < to->mesh_vertices.size() + 1; i++){
        MVertex *v0 = (i == 0) ?
          to->getBeginVertex()->mesh_vertices[0] : to->mesh_vertices[i - 1];
        MVertex *v1 = (i == to->mesh_vertices.size()) ?
          to->getEndVertex()->mesh_vertices[0] : to->mesh_vertices[i];
        to->lines.push_back(new MLine(v0, v1));
      }
    }
    
    void deMeshGEdge::operator() (GEdge *ge)
    {
      if(ge->geomType() == GEntity::DiscreteCurve && !CTX::instance()->meshDiscrete)
        return;
      ge->deleteMesh();
      ge->meshStatistics.status = GEdge::PENDING;
      ge->correspondingVertices.clear();
    }
    
    /*
    static void printFandPrimitive(int tag, std::vector<IntPoint> &Points)
    {
      char name[256];
      sprintf(name, "line%d.dat", tag);
      FILE *f = Fopen(name, "w");
      if(!f) return;
      double l = 0;
      for (unsigned int i = 0; i < Points.size(); i++){
        const IntPoint &P = Points[i];
        if (i) l +=(P.t - Points[i-1].t)*P.xp;
        fprintf(f, "%g %g %g %g %g\n", P.t, P.xp/P.lc, P.p,P.lc, l);
      }
      fclose(f);
    }
    */
    
    // new algo for recombining + splitting
    static int increaseN (int N)
    {
      if (((N+1)/2 - 1) % 2 != 0) return N+2;
      return N;
    }
    
    // ensure not to have points that are too close to each other.
    // can be caused by a coarse 1D mesh or by a noisy curve
    static void filterPoints(GEdge*ge, int nMinimumPoints)
    {
      if(ge->mesh_vertices.empty()) return;
      if(ge->meshAttributes.method == MESH_TRANSFINITE) return;
      //if (ge->mesh_vertices.size() <=3)return;
      bool forceOdd = false;
      if((ge->meshAttributes.method != MESH_TRANSFINITE ||
          CTX::instance()->mesh.flexibleTransfinite) &&
         CTX::instance()->mesh.algoRecombine != 0){
        if(CTX::instance()->mesh.recombineAll){
          forceOdd = true;
        }
      }
    
      MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
      std::vector<std::pair<double, MVertex*> > lengths;
      for (unsigned int i=0;i<ge->mesh_vertices.size();i++){
        MEdgeVertex *v = dynamic_cast<MEdgeVertex*> (ge->mesh_vertices[i]);
        if (!v){
          Msg::Error("in 1D mesh filterPoints");
          return;
        }
        double d = distance (v,v0);
        double t;
        v->getParameter(0,t);
        if (i != 0){
          double t0;
          v0->getParameter(0,t0);
          t=0.5*(t+t0);
        }
        double lc = F_LcB(ge, t);
        // double lc = v->getLc();
        if (d < lc * .3) {
          lengths.push_back(std::make_pair(lc/d,v));
        }
        else
          v0=v;
      }
      std::sort(lengths.begin(),lengths.end());
      int last = lengths.size();
      if (forceOdd) {
        while (last %2 != 0)last--;
      }
      /*
        if (CTX::instance()->mesh.algoRecombine == 2){
        if (last < 4)last = 0;
          while (last %4 != 0)last--;
          }
        else {
          while (last %2 != 0)last--;
          }
        }
      */
    
      bool filteringObservesMinimumN = (((int)ge->mesh_vertices.size() - last) >= nMinimumPoints);
      if (filteringObservesMinimumN){
        for (int i = 0; i < last; i++){
          std::vector<MVertex*>::iterator it = std::find(ge->mesh_vertices.begin(),
                                                         ge->mesh_vertices.end(),
                                                         lengths[i].second);
          ge->mesh_vertices.erase(it);
          delete lengths[i].second;
        }
      }
    }
    
    static void createPoints ( GVertex *gv, GEdge *ge, BoundaryLayerField *blf ,
    			   std::vector<MVertex*>& v, const SVector3 &dir){
      double L = blf->hwall_n;
    
      double LEdge = distance (ge->getBeginVertex()->mesh_vertices[0],
    			   ge->getEndVertex()->mesh_vertices[0]);
    
      while (1){
        if (L > blf->thickness || L > LEdge * .4) break;
        SPoint3 p (gv->x() + dir.x() * L, gv->y() + dir.y() * L, 0.0);
        v.push_back(new MEdgeVertex (p.x(), p.y(), p.z(), ge,  ge->parFromPoint(p), blf->hfar));
        int ith = v.size() ;
        L+= blf->hwall_n * pow (blf->ratio, ith);
        //    printf("parameter %g length %g\n",ge->parFromPoint(p),L);
      }
    }
    
    void addBoundaryLayerPoints (GEdge *ge,
    			     double &t_begin, // may change the left  parameter of the interval
    			     double &t_end,   // may change the right parameter of the interval
    			     std::vector<MVertex*> &_addBegin, // additional points @ left
    			     std::vector<MVertex*> &_addEnd)   // additional points @ right
    {
      BoundaryLayerField *blf = 0;
    #if defined(HAVE_ANN)
      FieldManager *fields = ge->model()->getFields();
      Field *bl_field = fields->get(fields->getBoundaryLayerField());
      blf = dynamic_cast<BoundaryLayerField*> (bl_field);
      if (blf) blf->setupFor1d(ge->tag());
    #endif
      if (!blf) return;
      if (blf->isEdgeBL(ge->tag()))return;
      SVector3 dir ( ge->getEndVertex()->x() - ge->getBeginVertex()->x(),
    		 ge->getEndVertex()->y() - ge->getBeginVertex()->y(),
    		 ge->getEndVertex()->z() - ge->getBeginVertex()->z());
      dir.normalize();
      GVertex *gvb = ge->getBeginVertex();
      GVertex *gve = ge->getEndVertex();
      if (blf->isEndNode(gvb->tag())){
        if (ge->geomType() != GEntity::Line){
          Msg::Error ("Boundary layer end point %d should lie on a straight line", gvb->tag());
          return;
        }
        createPoints (gvb, ge, blf, _addBegin, dir);
        if (!_addBegin.empty())_addBegin[_addBegin.size()-1]->getParameter(0,t_begin);
      }
      if (blf->isEndNode(gve->tag())){
        if (ge->geomType() != GEntity::Line){
          Msg::Error ("Boundary layer end point %d should lie on a straight line", gve->tag());
          return;
        }
        createPoints (gve, ge, blf, _addEnd, dir * -1.0);
        if (!_addEnd.empty())_addEnd[_addEnd.size()-1]->getParameter(0,t_end);
      }
    
      //  printf("Edge %d § %d %d points added (%g %g)-\n", ge->tag(), _addBegin.size(),_addEnd.size(),t_begin,t_end);
    }
    
    void meshGEdge::operator() (GEdge *ge)
    {
      /*
      BoundaryLayerField *blf = 0;
    #if defined(HAVE_ANN)
      FieldManager *fields = ge->model()->getFields();
      Field *bl_field = fields->get(fields->getBoundaryLayerField());
      blf = dynamic_cast<BoundaryLayerField*> (bl_field);
      if (blf) blf->setupFor1d(ge->tag());
    #endif
      */
    
      ge->model()->setCurrentMeshEntity(ge);
    
      //  if(ge->geomType() == GEntity::DiscreteCurve) return;
      if(ge->geomType() == GEntity::BoundaryLayerCurve) return;
      if(ge->meshAttributes.method == MESH_NONE) return;
      if(CTX::instance()->mesh.meshOnlyVisible && !ge->getVisibility()) return;
    
      // look if we are doing the STL triangulation
      std::vector<MVertex*> &mesh_vertices = ge->mesh_vertices ;
      std::vector<MLine*> &lines = ge->lines ;
    
      deMeshGEdge dem;
      dem(ge);
    
      if(MeshExtrudedCurve(ge)) return;
    
      if (ge->meshMaster() != ge){
        GEdge *gef = dynamic_cast<GEdge*> (ge->meshMaster());
        if (gef->meshStatistics.status == GEdge::PENDING) return;
        Msg::Info("Meshing curve %d (%s) as a copy of %d", ge->tag(),
                  ge->getTypeString().c_str(), ge->meshMaster()->tag());
        copyMesh(gef, ge, ge->masterOrientation);
        ge->meshStatistics.status = GEdge::DONE;
        return;
      }
    
      Msg::Info("Meshing curve %d (%s)", ge->tag(), ge->getTypeString().c_str());
    
      // compute bounds
      Range<double> bounds = ge->parBounds(0);
      double t_begin = bounds.low();
      double t_end = bounds.high();
    
      // if a BL is ending at one of the ends, then create specific points
      std::vector<MVertex*> _addBegin, _addEnd;
      addBoundaryLayerPoints (ge, t_begin, t_end, _addBegin, _addEnd);
    
      // first compute the length of the curve by integrating one
      double length;
      std::vector<IntPoint> Points;
      if(ge->geomType() == GEntity::Line &&
          ge->getBeginVertex() == ge->getEndVertex() &&
          //do not consider closed lines as degenerated
          (ge->position(0.5) - ge->getBeginVertex()->xyz()).norm() < CTX::instance()->geom.tolerance)
        length = 0.; // special case t avoid infinite loop in integration
      else
        length = Integration(ge, t_begin, t_end, F_One, Points, 1.e-8 * CTX::instance()->lc);
      ge->setLength(length);
      Points.clear();
    
      if(length < CTX::instance()->mesh.toleranceEdgeLength){
        ge->setTooSmall(true);
      }
    
      // Integrate detJ/lc du
      double a;
      int N;
      int filterMinimumN = 1;
      if(length == 0. && CTX::instance()->mesh.toleranceEdgeLength == 0.){
        Msg::Warning("Curve %d has a zero length", ge->tag());
        a = 0.;
        N = 1;
      }
      else if(ge->degenerate(0)){
        a = 0.;
        N = 1;
      }
      else if(ge->meshAttributes.method == MESH_TRANSFINITE){
        a = Integration(ge, t_begin, t_end, F_Transfinite, Points,
                        CTX::instance()->mesh.lcIntegrationPrecision);
        N = ge->meshAttributes.nbPointsTransfinite;
        if(CTX::instance()->mesh.flexibleTransfinite && CTX::instance()->mesh.lcFactor)
          N /= CTX::instance()->mesh.lcFactor;
      }
      else{
        if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG/* || blf*/){
          a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
                          CTX::instance()->mesh.lcIntegrationPrecision);
        }
        else{
           a = Integration(ge, t_begin, t_end, F_Lc, Points,
                          CTX::instance()->mesh.lcIntegrationPrecision);
        }
    
        // we should maybe provide an option to disable the smoothing
        for (unsigned int i = 0; i < Points.size(); i++){
          IntPoint &pt = Points[i];
          SVector3 der = ge->firstDer(pt.t);
          pt.xp = der.norm();
        }
        a = smoothPrimitive(ge, sqrt(CTX::instance()->mesh.smoothRatio), Points);
        filterMinimumN = ge->minimumMeshSegments() + 1;
        N = std::max(filterMinimumN, (int)(a + 1.99));
      }
    
      // force odd number of points if blossom is used for recombination
      if((ge->meshAttributes.method != MESH_TRANSFINITE ||
          CTX::instance()->mesh.flexibleTransfinite) &&
         CTX::instance()->mesh.algoRecombine != 0){
        if(CTX::instance()->mesh.recombineAll){
          if (N % 2 == 0) N++;
          if (CTX::instance()->mesh.algoRecombine == 2)
    	N = increaseN(N);
        }
        else{
          std::list<GFace*> faces = ge->faces();
          for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
            if((*it)->meshAttributes.recombine){
    	  if (N % 2 == 0) N ++;
    	  if (CTX::instance()->mesh.algoRecombine == 2)
    	    N = increaseN(N);
              break;
            }
          }
        }
      }
    
    
      //printFandPrimitive(ge->tag(),Points);
    
      // if the curve is periodic and if the begin vertex is identical to
      // the end vertex and if this vertex has only one model curve
      // adjacent to it, then the vertex is not connecting any other
      // curve. So, the mesh vertex and its associated geom vertex are not
      // necessary at the same location
      GPoint beg_p, end_p;
      if(ge->getBeginVertex() == ge->getEndVertex() &&
         ge->getBeginVertex()->edges().size() == 1){
        end_p = beg_p = ge->point(t_begin);
        Msg::Debug("Meshing periodic closed curve");
      }
      else{
        MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
        MVertex *v1 = ge->getEndVertex()->mesh_vertices[0];
        beg_p = GPoint(v0->x(), v0->y(), v0->z());
        end_p = GPoint(v1->x(), v1->y(), v1->z());
      }
    
      // do not consider the first and the last vertex (those are not
      // classified on this mesh edge)
      if(N > 1){
        const double b = a / (double)(N - 1);
        int count = 1, NUMP = 1;
        IntPoint P1, P2;
        mesh_vertices.resize(N - 2);
        while(NUMP < N - 1) {
          P1 = Points[count - 1];
          P2 = Points[count];
          const double d = (double)NUMP * b;
          if((fabs(P2.p) >= fabs(d)) && (fabs(P1.p) < fabs(d))) {
            double dt = P2.t - P1.t;
            double dlc = P2.lc - P1.lc;
            double dp = P2.p - P1.p;
            double t   = P1.t + dt / dp * (d - P1.p);
            SVector3 der = ge->firstDer(t);
            const double d = norm(der);
            double lc  = d/(P1.lc + dlc / dp * (d - P1.p));
            GPoint V = ge->point(t);
    	// printf("%d %g\n",NUMP-1,t);
            mesh_vertices[NUMP - 1] = new MEdgeVertex(V.x(), V.y(), V.z(), ge, t, lc);
            NUMP++;
          }
          else {
            count++;
          }
        }
        mesh_vertices.resize(NUMP - 1);
      }
    
      // Boundary Layer points are added
      {
        std::vector<MVertex*> vv;
        vv.insert(vv.end(), _addBegin.begin(), _addBegin.end());
        vv.insert(vv.end(), mesh_vertices.begin(), mesh_vertices.end());
        for (unsigned int i=0;i<_addEnd.size();i++)
          vv.push_back(_addEnd[_addEnd.size()-1-i]);
        //    vv.insert(vv.end(), _addEnd.rend(), _addEnd.rbegin());
        mesh_vertices = vv;
      }
    
      //  printf("%ld ----> ", ge->mesh_vertices.size());
      if (_addBegin.empty() && _addEnd.empty())
        filterPoints (ge, filterMinimumN - 2);
      //  printf("%ld \n", ge->mesh_vertices.size());
    
      for(unsigned int i = 0; i < mesh_vertices.size() + 1; i++){
        MVertex *v0 = (i == 0) ?
          ge->getBeginVertex()->mesh_vertices[0] : mesh_vertices[i - 1];
        MVertex *v1 = (i == mesh_vertices.size()) ?
          ge->getEndVertex()->mesh_vertices[0] : mesh_vertices[i];
        lines.push_back(new MLine(v0, v1));
      }
    
      if(ge->getBeginVertex() == ge->getEndVertex() &&
         ge->getBeginVertex()->edges().size() == 1){
        MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
        v0->x() = beg_p.x();
        v0->y() = beg_p.y();
        v0->z() = beg_p.z();
      }
    
      ge->meshStatistics.status = GEdge::DONE;
    }
    
    void orientMeshGEdge::operator()(GEdge *ge)
    {
      // apply user-specified mesh orientation constraints
      if(ge->meshAttributes.reverseMesh)
        for(unsigned int k = 0; k < ge->getNumMeshElements(); k++)
          ge->getMeshElement(k)->reverse();
    }