mesh drawing debugged

INRIA input corrected and enhanced

Geo/GFace.cpp

-// $Id: GFace.cpp,v 1.67 2008-07-04 12:03:50 geuzaine Exp $
+// $Id: GFace.cpp,v 1.68 2008-07-08 12:44:33 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -730,7 +730,21 @@ bool GFace::buildSTLTriangulation()
 // by default we assume that straight lines are geodesics
 SPoint2 GFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
 {
-  return pt1 + (pt2 - pt1) * t;
+  if(CTX.mesh.second_order_experimental){
+    // FIXME: this is buggy -- remove the CTX option once we do it in
+    // a robust manner
+    GPoint gp1 = point(pt1.x(), pt1.y());
+    GPoint gp2 = point(pt2.x(), pt2.y());
+    SPoint2 guess = pt1 + (pt2 - pt1) * t;
+    GPoint gp = closestPoint(SPoint3(gp1.x() + t * (gp2.x() - gp1.x()),
+				     gp1.y() + t * (gp2.y() - gp1.y()),
+				     gp1.z() + t * (gp2.z() - gp1.z())),
+			     (double*)guess);
+    return SPoint2(gp.u(), gp.v());
+  }
+  else{
+    return pt1 + (pt2 - pt1) * t;
+  }
 }
 
 // length of a curve drawn on a surface

Geo/GModelIO_Mesh.cpp

-// $Id: GModelIO_Mesh.cpp,v 1.57 2008-07-04 16:58:48 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.58 2008-07-08 12:44:33 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -1396,18 +1396,19 @@ int GModel::readMESH(const std::string &name)
   char str[256];
   int format;
   sscanf(buffer, "%s %d", str, &format);
-
   if(format != 1){
     Msg::Error("Medit mesh import only available for ASCII files");
     return 0;
   }
 
   std::vector<MVertex*> vertexVector;
-  std::map<int, std::vector<MElement*> > elements[3];
+  std::map<int, std::vector<MElement*> > elements[4];
+  std::vector<MVertex*> corners,ridges;
 
   while(!feof(fp)) {
-    if(!fgets(buffer, sizeof(buffer), fp)) break;
-    if(buffer[0] != '#'){ // skip comments
+    if(!fgets(buffer, 256, fp)) break;
+    if(buffer[0] != '#'){ // skip comments and empty lines
+      str[0]='\0';
       sscanf(buffer, "%s", str);
       if(!strcmp(str, "Dimension")){
         if(!fgets(buffer, sizeof(buffer), fp)) break;
@@ -1426,6 +1427,21 @@ int GModel::readMESH(const std::string &name)
           vertexVector[i] = new MVertex(x, y, z);
         }
       }
+      else if(!strcmp(str, "Edges")){
+        if(!fgets(buffer, sizeof(buffer), fp)) break;
+        int nbe;
+        sscanf(buffer, "%d", &nbe);
+        Msg::Info("%d edges", nbe);
+        for(int i = 0; i < nbe; i++) {
+          if(!fgets(buffer, sizeof(buffer), fp)) break;
+          int n[2], cl;
+          sscanf(buffer, "%d %d %d", &n[0], &n[1], &cl);
+          for(int j = 0; j < 2; j++) n[j]--;
+          std::vector<MVertex*> vertices;
+          if(!getVertices(2, n, vertexVector, vertices)) return 0;
+          elements[3][cl].push_back(new MLine(vertices));
+        }
+      }
       else if(!strcmp(str, "Triangles")){
         if(!fgets(buffer, sizeof(buffer), fp)) break;
         int nbe;
@@ -1441,6 +1457,36 @@ int GModel::readMESH(const std::string &name)
           elements[0][cl].push_back(new MTriangle(vertices));
         }
       }
+      else if(!strcmp(str, "Corners")){
+        if(!fgets(buffer, sizeof(buffer), fp)) break;
+        int nbe;
+        sscanf(buffer, "%d", &nbe);
+        Msg::Info("%d corners", nbe);
+        for(int i = 0; i < nbe; i++) {
+          if(!fgets(buffer, sizeof(buffer), fp)) break;
+          int  n[1];
+          sscanf(buffer, "%d",&n[0]);
+          for(int j = 0; j < 1; j++) n[j]--;
+	  //          std::vector<MVertex*> vertices;
+	  //          if(!getVertices(1, n, vertexVector, vertices)) return 0;
+	  //          corners.push_back(vertices[0]);
+        }
+      }
+      else if(!strcmp(str, "Ridges")){
+        if(!fgets(buffer, sizeof(buffer), fp)) break;
+        int nbe;
+        sscanf(buffer, "%d", &nbe);
+        Msg::Info("%d ridges", nbe);
+        for(int i = 0; i < nbe; i++) {
+          if(!fgets(buffer, sizeof(buffer), fp)) break;
+          int  n[1];
+          sscanf(buffer, "%d",&n[0]);
+          for(int j = 0; j < 1; j++) n[j]--;
+	  //          std::vector<MVertex*> vertices;
+	  //          if(!getVertices(1, n, vertexVector, vertices)) return 0;
+	  //          ridges.push_back(vertices[0]);
+        }
+      }
       else if(!strcmp(str, "Quadrilaterals")) {
         if(!fgets(buffer, sizeof(buffer), fp)) break;
         int nbe;

Geo/MElement.cpp

-// $Id: MElement.cpp,v 1.77 2008-07-03 17:06:01 geuzaine Exp $
+// $Id: MElement.cpp,v 1.78 2008-07-08 12:44:33 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -690,27 +690,29 @@ void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, double ww, SPo
   double sf[256];
   int nf = getNumFaceVertices();
 
+  //  printf("%d %d\n",nf,ord);
   if (!nf){
     switch(ord){
-    case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, ww,sf); break;
-    case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, ww,sf); break;
-    case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv, ww,sf); break;
-    case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv, ww,sf); break;
-    case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv, ww,sf); break;
+    case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv,sf); break;
+    case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv,sf); break;
+    case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv,sf); break;
+    case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv,sf); break;
+    case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv,sf); break;
     default: Msg::Error("Order %d triangle pnt not implemented", ord); break;
     }
   }
   else{
     switch(ord){
-    case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, ww,sf); break;
-    case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, ww,sf); break;
-    case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv, ww,sf); break;
-    case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv, ww,sf); break;
-    case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv, ww,sf); break;
+    case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv,sf); break;
+    case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv,sf); break;
+    case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv,sf); break;
+    case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv,sf); break;
+    case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv,sf); break;
     default: Msg::Error("Order %d triangle pnt not implemented", ord); break;
     }
   }
-  
+  //  printf("coucou\n");
+
   double x = 0 ; for(int i = 0; i < 3; i++) x += sf[i] * _v[i]->x();
   double y = 0 ; for(int i = 0; i < 3; i++) y += sf[i] * _v[i]->y();
   double z = 0 ; for(int i = 0; i < 3; i++) z += sf[i] * _v[i]->z();

Geo/OCCFace.cpp

-// $Id: OCCFace.cpp,v 1.43 2008-07-03 17:06:02 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.44 2008-07-08 12:44:33 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -157,7 +157,7 @@ GPoint OCCFace::closestPoint(const SPoint3 & qp, const double initialGuess[2]) c
   double pp[2] = {initialGuess[0],initialGuess[1]};
   proj.LowerDistanceParameters(pp[0], pp[1]);
 
-  Msg::Info("projection lower distance parameters %g %g",pp[0],pp[1]);
+  //  Msg::Info("projection lower distance parameters %g %g",pp[0],pp[1]);
 
   if((pp[0] < umin || umax < pp[0]) || (pp[1]<vmin || vmax<pp[1])){
     Msg::Error("Point projection is out of face bounds");

Geo/gmshFace.cpp

-// $Id: gmshFace.cpp,v 1.65 2008-07-03 17:06:02 geuzaine Exp $
+// $Id: gmshFace.cpp,v 1.66 2008-07-08 12:44:33 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -248,25 +248,6 @@ GEntity::GeomType gmshFace::geomType() const
   }
 }
 
-// by default we assume that straight lines are geodesics
-SPoint2 gmshFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
-{
-  if(isSphere && CTX.mesh.second_order_experimental){
-    // FIXME: this is buggy -- remove the CTX option once we do it in
-    // a robust manner
-    GPoint gp1 = point(pt1.x(), pt1.y());
-    GPoint gp2 = point(pt2.x(), pt2.y());
-    SPoint2 guess = pt1 + (pt2 - pt1) * t;
-    GPoint gp = closestPoint(SPoint3(gp1.x() + t * (gp2.x() - gp1.x()),
-				     gp1.y() + t * (gp2.y() - gp1.y()),
-				     gp1.z() + t * (gp2.z() - gp1.z())),
-			     (double*)guess);
-    return SPoint2(gp.u(), gp.v());
-  }
-  else{
-    return pt1 + (pt2 - pt1) * t;
-  }
-}
 
 bool gmshFace::containsPoint(const SPoint3 &pt) const
 { 

Geo/gmshFace.h

@@ -34,7 +34,6 @@ class gmshFace : public GFace {
   virtual ~gmshFace(){}
   Range<double> parBounds(int i) const; 
   void setModelEdges(std::list<GEdge*>&);
-  virtual SPoint2 geodesic(const SPoint2 &pt1, const SPoint2 &pt2, double t);
   virtual GPoint point(double par1, double par2) const; 
   virtual GPoint closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const; 
   virtual bool containsPoint(const SPoint3 &pt) const;  

Mesh/HighOrder.cpp

-// $Id: HighOrder.cpp,v 1.32 2008-07-03 17:06:03 geuzaine Exp $
+// $Id: HighOrder.cpp,v 1.33 2008-07-08 12:44:34 remacle Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -1358,7 +1358,7 @@ void printJacobians(GModel *m, const char *nm)
 {
   return;
 
-  const int n = 5;
+  const int n = 15;
   double D[n][n];
   double X[n][n];
   double Y[n][n];

Numeric/FunctionSpace.h

@@ -36,7 +36,6 @@ struct gmshFunctionSpace
       p[j][1] = pow(vv,monomials(j, 1));
     }
   }
-
   inline void f(double u, double v, double w, double *sf) const
   {
     for(int i = 0; i < coefficients.size1(); i++){