From 6dc860f92fac902a5e0d7fad5ee2ce9c7b7cbf96 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 19 Nov 2008 21:22:45 +0000
Subject: [PATCH] fix high order extruded meshes

---
 Geo/MVertex.cpp    | 17 +++++++++--------
 Geo/MVertex.h      |  4 ++--
 Geo/gmshFace.cpp   |  2 +-
 Mesh/HighOrder.cpp |  4 ++--
 4 files changed, 14 insertions(+), 13 deletions(-)

diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index a16d3e84b4..18515fba75 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -241,8 +241,8 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params
   }
 }
 
-bool reparamMeshVerticesOnFace(MVertex *v1, MVertex *v2, GFace *gf, 
-                               SPoint2 &param1, SPoint2 &param2)
+bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf, 
+                           SPoint2 &param1, SPoint2 &param2)
 {
   std::vector<SPoint2> p1, p2;
   getAllParameters(v1, gf, p1);
@@ -250,7 +250,6 @@ bool reparamMeshVerticesOnFace(MVertex *v1, MVertex *v2, GFace *gf,
   if (p1.size() == 1 && p2.size() == 1){
     param1 = p1[0];
     param2 = p2[0];
-    return true;
   }
   else if (p1.size() == 1 && p2.size() == 2){
     double d1 = 
@@ -261,7 +260,6 @@ bool reparamMeshVerticesOnFace(MVertex *v1, MVertex *v2, GFace *gf,
       (p1[0].x() - p2[1].y()) * (p1[0].y() - p2[1].y());
     param1 = p1[0];
     param2 = d2 < d1 ? p2[1] : p2[0];
-    return true;
   }  
   else if (p2.size() == 1 && p1.size() == 2){
     double d1 = 
@@ -272,9 +270,12 @@ bool reparamMeshVerticesOnFace(MVertex *v1, MVertex *v2, GFace *gf,
       (p2[0].x() - p1[1].y()) * (p2[0].y() - p1[1].y());
     param1 = d2 < d1 ? p1[1] : p1[0];
     param2 = p2[0];
-    return true;
   }  
-  return false;
+  else{
+    param1 = gf->parFromPoint(SPoint3(v1->x(), v1->y(), v1->z()));
+    param2 = gf->parFromPoint(SPoint3(v2->x(), v2->y(), v2->z()));
+  }
+  return true;
 }
 
 bool reparamMeshVertexOnFace(MVertex *v, GFace *gf, SPoint2 &param)
@@ -295,7 +296,7 @@ bool reparamMeshVertexOnFace(MVertex *v, GFace *gf, SPoint2 &param)
     GVertex *gv = (GVertex*)v->onWhat();
     param = gv->reparamOnFace(gf, 1);
 
-    // abort if we could be on a seam
+    // shout if we could be on a seam
     std::list<GEdge*> ed = gv->edges();
     for(std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); it++)
       if((*it)->isSeam(gf)) return false;
@@ -306,7 +307,7 @@ bool reparamMeshVertexOnFace(MVertex *v, GFace *gf, SPoint2 &param)
     v->getParameter(0, t);
     param = ge->reparamOnFace(gf, t, 1);
 
-    // abort if we are on a seam (todo: try dir=-1 and compare)
+    // shout if we are on a seam
     if(ge->isSeam(gf))
       return false;
   }
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index bb90b48963..ab4cf5fa8e 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -165,8 +165,8 @@ class MFaceVertex : public MVertex{
   }
 };
 
-bool reparamMeshVerticesOnFace(MVertex *v1, MVertex *v2, GFace *gf, 
-                               SPoint2 &param1, SPoint2 &param2);
+bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf, 
+                           SPoint2 &param1, SPoint2 &param2);
 bool reparamMeshVertexOnFace(MVertex *v, GFace *gf, SPoint2 &param);
 bool reparamMeshVertexOnEdge(MVertex *v, GEdge *ge, double &param);
 
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 9094059d3c..f72229c98f 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -248,7 +248,7 @@ SPoint2 gmshFace::parFromPoint(const SPoint3 &qp) const
   if(s->Typ == MSH_SURF_PLAN){
     double x, y, z, VX[3], VY[3];
     getMeanPlaneData(VX, VY, x, y, z);
-    double u, v, vec[3] = {qp.x()-x, qp.y()-y, qp.z()-z};
+    double u, v, vec[3] = {qp.x() - x, qp.y() - y, qp.z() - z};
     prosca(vec, VX, &u);
     prosca(vec, VY, &v);
     return SPoint2(u, v); 
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index c3042e4b0c..6c8fbba459 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -317,7 +317,7 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
       if(!linear && 
          gf->geomType() != GEntity::DiscreteSurface &&
          gf->geomType() != GEntity::BoundaryLayerSurface){
-        reparamOK = reparamMeshVerticesOnFace(v0, v1, gf, p0, p1);
+        reparamOK = reparamMeshEdgeOnFace(v0, v1, gf, p0, p1);
       }
       double US[100], VS[100];
       if(reparamOK && !linear && gf->geomType() != GEntity::DiscreteSurface){
@@ -336,7 +336,7 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
 	  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()));
+	    displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
 	  }
         }
         temp.push_back(v);
-- 
GitLab