fix high order mesh for boundary layers

Mesh/HighOrder.cpp

@@ -226,6 +226,10 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
                             int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
                             gmshHighOrderSmoother *displ3D = 0)
 {
+  if(ge->geomType() == GEntity::DiscreteCurve ||
+     ge->geomType() == GEntity::BoundaryLayerCurve)
+    linear = true;
+
   for(int i = 0; i < ele->getNumEdges(); i++){
     MEdge edge = ele->getEdge(i);
     std::pair<MVertex*, MVertex*> p(edge.getMinVertex(), edge.getMaxVertex());
@@ -237,42 +241,37 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
     }
     else{
       MVertex *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);            
-      double u0 = 0., u1 = 0.;
+      double u0 = 0., u1 = 0., US[100];
       bool reparamOK = true;
-      if(!linear && ge->geomType() != GEntity::DiscreteCurve &&
-         ge->geomType() != GEntity::BoundaryLayerCurve){
+      if(!linear){
         reparamOK &= reparamMeshVertexOnEdge(v0, ge, u0);
-        if (ge->periodic(0) && v1 == ge->getEndVertex()->mesh_vertices[0]){
-          Range<double> par = ge->parBounds(0);
-          u1 = par.high();
-        }         
+        if(ge->periodic(0) && v1 == ge->getEndVertex()->mesh_vertices[0])
+          u1 = ge->parBounds(0).high();
         else
           reparamOK &= reparamMeshVertexOnEdge(v1, ge, u1);
+        if(reparamOK){
+          double relax = 1.;
+          while (1){
+            if(computeEquidistantParameters(ge, u0, u1, nPts + 2, US, relax)) 
+              break;
+            relax /= 2.0;
+            if(relax < 1.e-2) 
+              break;
+          } 
+          if(relax < 1.e-2)
+            Msg::Warning("Failed to compute equidistant parameters (relax = %g)",
+                         relax);
+        }
       }
-      double US[100];
-      if(reparamOK && !linear && ge->geomType() != GEntity::DiscreteCurve){
-        double relax = 1.;
-        while (1){
-          if(computeEquidistantParameters(ge, u0, u1, nPts + 2, US, relax)) 
-            break;
-          relax /= 2.0;
-          if (relax < 1.e-2) 
-            break;
-        } 
-        if (relax < 1.e-2)
-          Msg::Warning("Failed to compute equidistant parameters (relax = %g)",
-                       relax);
-      }
-      std::vector<MVertex*> temp;      
+      std::vector<MVertex*> temp;
       for(int j = 0; j < nPts; j++){
         const double t = (double)(j + 1)/(nPts + 1);
-        double uc = (1. - t) * u0 + t * u1;
+        double uc = (1. - t) * u0 + t * u1; // can be wrong, that's ok
         MVertex *v;
-        if(!reparamOK || linear || ge->geomType() == GEntity::DiscreteCurve || 
-           uc < u0 || uc > u1){ // need to treat periodic curves properly!
+        if(linear || !reparamOK || uc < u0 || uc > u1){ 
+          // we don't have a (valid) parameter on the curve
           SPoint3 pc = edge.interpolate(t);
           v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
-          v->setParameter(0, t);
         }
         else {
           GPoint pc = ge->point(US[j + 1]);
@@ -284,6 +283,7 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
 	  }
         }
         temp.push_back(v);
+        // this destroys the ordering of the mesh vertices on the edge
         ge->mesh_vertices.push_back(v);
         ve.push_back(v);
       }
@@ -300,6 +300,10 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
                             int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
                             gmshHighOrderSmoother *displ3D = 0)
 {
+  if(gf->geomType() == GEntity::DiscreteSurface ||
+     gf->geomType() == GEntity::BoundaryLayerSurface)
+    linear = true;
+
   for(int i = 0; i < ele->getNumEdges(); i++){
     MEdge edge = ele->getEdge(i);    
     std::pair<MVertex*, MVertex*> p(edge.getMinVertex(), edge.getMaxVertex());
@@ -312,27 +316,26 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
     else{
       MVertex *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);
       SPoint2 p0, p1;
+      double US[100], VS[100];
       bool reparamOK = true;
-      if(!linear && 
-         gf->geomType() != GEntity::DiscreteSurface &&
-         gf->geomType() != GEntity::BoundaryLayerSurface){
+      if(!linear){
         reparamOK = reparamMeshEdgeOnFace(v0, v1, gf, p0, p1);
-      }
-      double US[100], VS[100];
-      if(reparamOK && !linear && gf->geomType() != GEntity::DiscreteSurface){
-        computeEquidistantParameters(gf, p0[0], p1[0], p0[1], p1[1], nPts + 2, US, VS);
+        if(reparamOK)
+          computeEquidistantParameters(gf, p0[0], p1[0], p0[1], p1[1], nPts + 2,
+                                       US, VS);
       }
       std::vector<MVertex*> temp;
       for(int j = 0; j < nPts; j++){
         const double t = (double)(j + 1) / (nPts + 1);
         MVertex *v;
-        if(!reparamOK || linear || gf->geomType() == GEntity::DiscreteSurface){
+        if(linear || !reparamOK){
+          // we don't have (valid) parameters on the surface
           SPoint3 pc = edge.interpolate(t);
           v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
         }
         else{
-          GPoint pc = gf->point(US[j+1], VS[j+1]);
-	  v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, US[j+1], VS[j+1]);
+          GPoint pc = gf->point(US[j + 1], VS[j + 1]);
+	  v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, US[j + 1], VS[j + 1]);
 	  if (displ3D){
 	    SPoint3 pc2 = edge.interpolate(t);          
 	    displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
@@ -390,6 +393,10 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
                             bool linear, int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
                             gmshHighOrderSmoother *displ3D = 0)
 {
+  if(gf->geomType() == GEntity::DiscreteSurface ||
+     gf->geomType() == GEntity::BoundaryLayerSurface)
+    linear = true;
+
   Double_Matrix points;
   int start = 0;
 
@@ -421,19 +428,16 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
       std::vector<MVertex*> &vtcs = faceVertices[face];
       SPoint2 pts[20];
       bool reparamOK = true;
-      if(!linear && 
-         gf->geomType() != GEntity::DiscreteSurface &&
-         gf->geomType() != GEntity::BoundaryLayerSurface){
-	for (int k = 0; k < incomplete->getNumVertices(); k++){
+      if(!linear){
+	for(int k = 0; k < incomplete->getNumVertices(); k++)
 	  reparamOK &= reparamMeshVertexOnFace(incomplete->getVertex(k), gf, pts[k]);
-	}
       }
       if(face.getNumVertices() == 3 && nPts > 1){ // tri face
         for(int k = start; k < points.size1(); k++){
           MVertex *v;
           const double t1 = points(k, 0);
           const double t2 = points(k, 1);
-          if(linear || gf->geomType() == GEntity::DiscreteSurface){
+          if(linear){
             SPoint3 pc = face.interpolate(t1, t2);
             v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
           }
@@ -463,7 +467,7 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
 	    else{
 	      v = new MVertex(X, Y, Z, gf);
 	    }
-	    if (displ3D){
+	    if(displ3D){
 	      SPoint3 pc2 = face.interpolate(t1, t2);
 	      displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
 	    }	    
@@ -481,7 +485,7 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
             // parameters are between -1 and 1
             double t1 = 2. * (double)(j + 1) / (nPts + 1) - 1.;
             double t2 = 2. * (double)(k + 1) / (nPts + 1) - 1.;
-            if(!reparamOK || linear || gf->geomType() == GEntity::DiscreteSurface){
+            if(linear || !reparamOK){
               SPoint3 pc = face.interpolate(t1, t2);
               v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
             }
@@ -697,7 +701,7 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
   for(unsigned int i = 0; i < gf->triangles.size(); i++){
     MTriangle *t = gf->triangles[i];
     std::vector<MVertex*> ve, vf;
-    getEdgeVertices(gf, t, ve, edgeVertices, linear, nPts,displ2D,displ3D);
+    getEdgeVertices(gf, t, ve, edgeVertices, linear, nPts, displ2D, displ3D);
     if(nPts == 1){
       triangles2.push_back
         (new MTriangle6(t->getVertex(0), t->getVertex(1), t->getVertex(2),

Mesh/Refine.cpp

@@ -12,6 +12,17 @@
 #include "GmshMessage.h"
 #include "OS.h"
 
+class MVertexLessParam{
+ public:
+  bool operator()(const MVertex *v1, const MVertex *v2) const
+  {
+    double u1 = 0., u2 = 1.;
+    v1->getParameter(0, u1);
+    v2->getParameter(0, u2);
+    return u1 < u2;
+  }
+};
+
 static void Subdivide(GEdge *ge)
 {
   std::vector<MLine*> lines2;
@@ -25,6 +36,9 @@ static void Subdivide(GEdge *ge)
   }
   ge->lines = lines2;
 
+  // 2nd order meshing destroyed the ordering of the vertices on the edge
+  std::sort(ge->mesh_vertices.begin(), ge->mesh_vertices.end(), 
+            MVertexLessParam());
   for(unsigned int i = 0; i < ge->mesh_vertices.size(); i++)
     ge->mesh_vertices[i]->setPolynomialOrder(1);
   ge->deleteVertexArrays();

Mesh/meshGFaceTransfinite.cpp

@@ -104,10 +104,8 @@ int MeshTransfiniteSurface(GFace *gf)
 
   // get the indices of the interpolation corners as well as the u,v
   // coordinates of all the boundary vertices
-  int iCorner = 0;
-  int N[4] = {0, 0, 0, 0};
-  std::vector<double> U;
-  std::vector<double> V;
+  int iCorner = 0, N[4] = {0, 0, 0, 0};
+  std::vector<double> U, V;
   for(unsigned int i = 0; i < m_vertices.size(); i++){
     MVertex *v = m_vertices[i];
     if(v == corners[0] || v == corners[1] || v == corners[2] || 
@@ -119,38 +117,15 @@ int MeshTransfiniteSurface(GFace *gf)
       }
     }
     SPoint2 param;
-    if(v->onWhat()->dim() == 0){
-      GVertex *gv = (GVertex*)v->onWhat();
-      param = gv->reparamOnFace(gf, 1);
-    }
-    else if(v->onWhat()->dim() == 1){
-      GEdge *ge = (GEdge*)v->onWhat();
-      double UU;
-      v->getParameter(0, UU);
-      param = ge->reparamOnFace(gf, UU, 1);
-    }
-    else{
-      double UU, VV;
-      if(v->onWhat() == gf && v->getParameter(0, UU) && v->getParameter(1, VV))
-        param = SPoint2(UU, VV);
-      else
-        param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
-    }
-    U.push_back(param.x());
-    V.push_back(param.y());
+    reparamMeshVertexOnFace(v, gf, param);
+    U.push_back(param[0]);
+    V.push_back(param[1]);
   }
 
-  int N1 = N[0];
-  int N2 = N[1];
-  int N3 = N[2];
-  int N4 = N[3];
-
-  int L = N2 - N1;
-  int H = N3 - N2;
-
+  int N1 = N[0], N2 = N[1], N3 = N[2], N4 = N[3];
+  int L = N2 - N1, H = N3 - N2;
   if(corners.size () == 4){
-    int Lb = N4 - N3;
-    int Hb = m_vertices.size() - N4;
+    int Lb = N4 - N3, Hb = m_vertices.size() - N4;
     if(Lb != L || Hb != H){
       Msg::Error("Surface %d cannot be meshed using the transfinite algo", 
 		 gf->tag());
@@ -166,10 +141,8 @@ int MeshTransfiniteSurface(GFace *gf)
     }      
   }
   
-  std::vector<double> lengths_i;
-  std::vector<double> lengths_j;
-  double L_i = 0;
-  double L_j = 0;
+  std::vector<double> lengths_i, lengths_j;
+  double L_i = 0, L_j = 0;
   lengths_i.push_back(0.);
   lengths_j.push_back(0.);
   for(int i = 0; i < L; i++){
@@ -230,12 +203,8 @@ int MeshTransfiniteSurface(GFace *gf)
     }
   }
 
-  double UC1 = U[N1];
-  double UC2 = U[N2];
-  double UC3 = U[N3];
-  double VC1 = V[N1];
-  double VC2 = V[N2];
-  double VC3 = V[N3];
+  double UC1 = U[N1], UC2 = U[N2], UC3 = U[N3];
+  double VC1 = V[N1], VC2 = V[N2], VC3 = V[N3];
 
   //create points using transfinite interpolation
   if(corners.size() == 4){