setting of the optimizer for P2 parameterization

Geo/discreteDiskFace.cpp

@@ -26,6 +26,8 @@
 #include "convexCombinationTerm.h"  // #FIXME
 
 #include "qualityMeasuresJacobian.h" // #temporary?
+#include "OptHomRun.h"
+#include "MeshQualityOptimizer.h"
 
 static inline void functionShapes(int p, double Xi[2], double* phi)
 {
@@ -286,13 +288,17 @@ discreteDiskFace::discreteDiskFace(int id, GFace *gf, triangulation* diskTriangu
   if (!checkOrientationUV()){
     Msg::Info("discreteDiskFace:: parametrization is not one-to-one; fixing "
               "the discrete system.");
-    parametrize(true);
+
+    //parametrize(true);
+    //buildOct(CAD); 
+    optimize();
     buildOct(CAD);    
-  }
-  
+    
+    }
   putOnView(true,false);
   printParamMesh();
   if(!checkOrientationUV()) Msg::Fatal("discreteDiskFace:: failing to fix the discrete system");
+  else Msg::Info("Parameterization done :-)");
 }
 /*end BUILDER*/
 
@@ -497,6 +503,87 @@ bool discreteDiskFace::checkOrientationUV()
   }
 }
 
+
+void discreteDiskFace::optimize(){
+
+  // parameters for mesh optimization
+  // -- high order
+  OptHomParameters optParams;
+  optParams.dim = 2;
+  optParams.optPassMax = 10;
+  optParams.TMAX = 300;
+  optParams.fixBndNodes = true;
+  optParams.strategy = 1;
+  // -- linear
+  MeshQualOptParameters opt;
+  opt.dim = 2;
+  opt.fixBndNodes = true;
+
+  //creating the "geometry" of the parametrization, and its corresponding mesh
+  // -- generation of parametric nodes
+  std::map<SPoint3,MVertex*> sp2mv;
+  std::vector<MElement*> paramTriangles;
+  for(std::map<MVertex*,SPoint3>::iterator it=coordinates.begin(); it!= coordinates.end(); ++it)
+    sp2mv[it->second] = new MVertex(it->second.x(),it->second.y(),0.);
+  // -- generation of parametric triangles    
+  paramTriangles.resize(discrete_triangles.size() - geoTriangulation->fillingHoles.size());
+  for(unsigned int i=0; i<discrete_triangles.size() -geoTriangulation->fillingHoles.size(); i++){
+    discreteDiskFaceTriangle* ct = &_ddft[i];
+    std::vector<MVertex*> mv;
+    mv.resize(ct->tri->getNumVertices());
+    for (int j=0; j<ct->tri->getNumVertices(); j++)
+      mv[j] = sp2mv[ct->p[j]];    
+    if(_order==1)
+      paramTriangles[i] = new MTriangle(mv);
+    else
+      paramTriangles[i] = new MTriangle6(mv);
+  }
+  // -- generation of the parametric topology #mark what about the GFace for the GModel ?
+  std::map<int,std::vector<MElement*> > e2e;
+  e2e[0] = paramTriangles;
+  std::vector<int> v;
+  v.push_back(0);
+  std::map<int,std::vector<int> > e2p;
+  e2p[0] = v;
+  GModel* paramDisk = GModel::createGModel(e2e,e2p);
+  discreteVertex dv(paramDisk,NEWPOINT());
+  dv.mesh_vertices.push_back(sp2mv[coordinates[_U0[0]]]);
+  sp2mv[coordinates[_U0[0]]]->setEntity(&dv);
+  paramDisk->add(&dv);
+  discreteEdge de(paramDisk,NEWLINE(),&dv,&dv);
+  paramDisk->add(&de);
+  for(unsigned int i=1; i<_U0.size(); i++){
+    sp2mv[coordinates[_U0[i]]]->setEntity(&de);
+    de.mesh_vertices.push_back(sp2mv[coordinates[_U0[i]]]);
+    de.lines.push_back(new MLine(sp2mv[coordinates[_U0[i-1]]],sp2mv[coordinates[_U0[i]]]));
+  }
+  de.createGeometry();
+
+  
+  // optimization
+  if(_order >1)
+    HighOrderMeshOptimizer(paramDisk, optParams);
+  else
+    MeshQualityOptimizer(paramDisk,opt);
+
+  // update the parametrization
+  paramTriangles = e2e[0];
+  for(unsigned int i=0; i< paramTriangles.size(); i++){
+    discreteDiskFaceTriangle* ct = &_ddft[i];
+    MElement* tri = paramTriangles[i];
+    for(int j=0; j<ct->tri->getNumVertices(); j++){
+      MVertex* mv = tri->getVertex(j);
+      SPoint3 p(mv->x(), mv->y(), 0);
+      coordinates[ct->tri->getVertex(j)] = p;
+      ct->p[j] = p;
+    }
+  }
+
+  // cleaning
+  delete paramDisk;
+  
+}
+
 // (u;v) |-> < (x;y;z); GFace; (u;v) >
 GPoint discreteDiskFace::point(const double par1, const double par2) const
 {
@@ -1010,7 +1097,7 @@ double triangulation::geodesicDistance () {
     addTo (v2t, tri[i]->getVertex(2),tri[i]);
   }
 
-  //  printf("computing geodesic distance with %d triangles and %d vertcie\n",
+  //  printf("computing geodesic distance with %d triangles and %d vertices\n",
   //	 tri.size(),v2t.size());
   
   std::map<MVertex*,double> Fixed;

Geo/discreteDiskFace.h

@@ -22,20 +22,7 @@
 #include "meshGFaceOptimize.h"
 #include "PView.h"
 #include "robustPredicates.h"
-
-inline double tri3Darea(MVertex* mv0, MVertex* mv1, MVertex* mv2)
-{
-
-  SVector3 v1(mv1->x()-mv0->x(),mv1->y()-mv0->y(),mv1->z()-mv0->z());
-  SVector3 v2(mv2->x()-mv0->x(),mv2->y()-mv0->y(),mv2->z()-mv0->z());
-  
-  SVector3 n(v1.y()*v2.z()-v2.y()*v1.z(),v2.x()*v1.z()-v1.x()*v2.z(),v1.x()*v2.y()-v2.x()*v1.y());
-
-  
-  return .5*n.norm();
-
-}
-
+#include "MLine.h"
 
 inline int nodeLocalNum(MElement* e, MVertex* v)
 {
@@ -231,6 +218,7 @@ class discreteDiskFace : public GFace {
   bool parametrize(bool one2one = false) const;
   void putOnView(bool,bool);
   bool checkOrientationUV();
+  void optimize();
   void printParamMesh();
 
  public:
@@ -274,6 +262,7 @@ class discreteDiskFace : public GFace {
   mutable discreteDiskFaceTriangle *_ddft;
   mutable Octree *oct;
   mutable std::vector<double> _coords;
+  
 };
 
 #endif

Geo/discreteFace.cpp

@@ -111,7 +111,7 @@ void discreteFace::secondDer(const SPoint2 &param,
 void discreteFace::createGeometry()
 {
   checkAndFixOrientation();
-  int order = 1;
+  int order = 2;
   int nPart = 2;
   double eta = 2*3.14/7;
 
@@ -154,12 +154,12 @@ void discreteFace::createGeometry()
 
   for(unsigned int i=0; i<toParam.size(); i++){
     //printf("MAP(%d) : aspect ratio = %12.5E\n",i,toParam[i]->aspectRatio());
-    char name[256];
+    //char name[256];
     //sprintf(name,"map%d.pos",i);
-    toParam[i]->print(name,i);
+    //toParam[i]->print(name,i);
     fillHoles(toParam[i]);
     //sprintf(name,"mapFilled%d.pos",i);
-    toParam[i]->print(name,i);
+    //toParam[i]->print(name,i);
   }
   for(unsigned int i=0; i<toParam.size(); i++){
     std::vector<MElement*> mytri = toParam[i]->tri;
@@ -167,6 +167,7 @@ void discreteFace::createGeometry()
     df->replaceEdges(toParam[i]->my_GEdges);
     _atlas.push_back(df);
   }
+
 #endif
 }
 
@@ -215,7 +216,6 @@ void discreteFace::gatherMeshes()
 
 void discreteFace::mesh(bool verbose)
 {
-  
 #if defined(HAVE_ANN) && defined(HAVE_SOLVER) && defined(HAVE_MESH)
   if (!CTX::instance()->meshDiscrete) return;
   for (unsigned int i=0;i<_atlas.size();i++){