hack to take into account periodicity when optimizing high order meshes....

hack to take into account periodicity when optimizing high order meshes. Periodic node positions are corrected after each optimization pass. Optimization is executed until no more everything is ok. 

Mesh/HighOrder.cpp

@@ -21,6 +21,7 @@
 #include "Context.h"
 #include "fullMatrix.h"
 #include "BasisFactory.h"
+#include "MVertexRTree.h"
 
 
 // --------- Functions that help optimizing placement of points on geometry -----------
@@ -1307,6 +1308,73 @@ static void updateHighOrderVertices(GEntity *e,
   e->deleteVertexArrays();
 }
 
+static void updatePeriodicEdgesAndFaces(GModel *m)
+{
+  for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it) {
+    GEdge *slave = *it;
+    GEdge *master = dynamic_cast<GEdge*>(slave->meshMaster());
+    if (master == slave) continue;
+
+    std::map<MVertex*,MVertex*> &v2v = slave->correspondingVertices;
+    v2v.clear();
+
+    MVertexRTree rtree = MVertexRTree(1e-8 * CTX::instance()->lc);
+    for (unsigned int i = 0; i < slave->getNumMeshVertices(); ++i)
+      rtree.insert(slave->getMeshVertex(i));
+
+    std::vector<double> tfo = slave->affineTransform;
+
+    for (unsigned int i = 0; i < master->getNumMeshVertices(); ++i) {
+      MVertex *vs = master->getMeshVertex(i);
+      double ps[4] = {vs->x(), vs->y(), vs->z(), 1.};
+      double res[4] = {0., 0., 0., 0.};
+      int idx = 0;
+      for(int i = 0; i < 4; i++)
+        for(int j = 0; j < 4; j++)
+          res[i] +=  tfo[idx++] * ps[j];
+
+      SPoint3 p3 (res[0], res[1], res[2]);
+      double u = slave->parFromPoint(p3);
+      GPoint gp = slave->point(u);
+      MVertex *vt = rtree.find(gp.x(), gp.y(), gp.z());
+      if (!vt) Msg::Error("Couldn't find a vertex for updating periodicity");
+      else v2v[vt] = vs;
+    }
+  }
+
+  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
+    GFace *slave = *it;
+    GFace *master = dynamic_cast<GFace*>(slave->meshMaster());
+    if (master == slave) continue;
+
+    std::map<MVertex*,MVertex*> &v2v = slave->correspondingVertices;
+    v2v.clear();
+
+    MVertexRTree rtree = MVertexRTree(1e-8 * CTX::instance()->lc);
+    for (unsigned int i = 0; i < slave->getNumMeshVertices(); ++i)
+      rtree.insert(slave->getMeshVertex(i));
+
+    std::vector<double> tfo = slave->affineTransform;
+
+    for (unsigned int i = 0; i < master->getNumMeshVertices(); ++i) {
+      MVertex *vs = master->getMeshVertex(i);
+      double ps[4] = {vs->x(), vs->y(), vs->z(), 1.};
+      double res[4] = {0., 0., 0., 0.};
+      int idx = 0;
+      for(int i = 0; i < 4; i++)
+        for(int j = 0; j < 4; j++)
+          res[i] +=  tfo[idx++] * ps[j];
+
+      SPoint3 p3 (res[0], res[1], res[2]);
+      SPoint2 p2 = slave->parFromPoint(p3);
+      GPoint gp = slave->point(p2);
+      MVertex *vt = rtree.find(gp.x(), gp.y(), gp.z());
+      if (!vt) Msg::Error("Couldn't find a vertex for updating periodicity");
+      else v2v[vt] = vs;
+    }
+  }
+}
+
 void SetOrder1(GModel *m,  bool onlyVisible)
 {
   m->destroyMeshCaches();
@@ -1334,6 +1402,8 @@ void SetOrder1(GModel *m,  bool onlyVisible)
     updateHighOrderVertices(*it, newHOVert, onlyVisible);
   for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
     updateHighOrderVertices(*it, newHOVert, onlyVisible);
+
+  updatePeriodicEdgesAndFaces(m);
 }
 
 void checkHighOrderTriangles(const char* cc, GModel *m,
@@ -1500,6 +1570,8 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
   for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
     updateHighOrderVertices(*it, newHOVert[*it], onlyVisible);
 
+  updatePeriodicEdgesAndFaces(m);
+
   double t2 = Cpu();
 
   std::vector<MElement*> bad;
@@ -1559,6 +1631,8 @@ void SetHighOrderComplete(GModel *m, bool onlyVisible)
     (*it)->mesh_vertices.clear();
     (*it)->mesh_vertices.insert((*it)->mesh_vertices.begin(), newV.begin(), newV.end());
   }
+
+  updatePeriodicEdgesAndFaces(m);
 }
 
 void SetHighOrderInComplete (GModel *m, bool onlyVisible)
@@ -1609,4 +1683,6 @@ void SetHighOrderInComplete (GModel *m, bool onlyVisible)
     }
     (*it)->mesh_vertices = newV;
   }
+
+  updatePeriodicEdgesAndFaces(m);
 }

contrib/HighOrderMeshOptimizer/CMakeLists.txt

@@ -13,6 +13,7 @@ set(SRC
   SuperEl.cpp 
   OptHomElastic.cpp
   OptHomFastCurving.cpp
+  OptHomPeriodicity.cpp
   CADDistances.cpp
 )
 

contrib/HighOrderMeshOptimizer/OptHomPeriodicity.cpp

+// Copyright (C) 2013 ULg-UCL
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use, copy,
+// modify, merge, publish, distribute, and/or sell copies of the
+// Software, and to permit persons to whom the Software is furnished
+// to do so, provided that the above copyright notice(s) and this
+// permission notice appear in all copies of the Software and that
+// both the above copyright notice(s) and this permission notice
+// appear in supporting documentation.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE
+// COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR
+// ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY
+// DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
+// WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
+// ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
+// OF THIS SOFTWARE.
+//
+// Please report all bugs and problems to the public mailing list
+// <gmsh@onelab.info>.
+//
+// Contributors: Amaury Johnen
+
+#include "OptHomPeriodicity.h"
+#include "GEntity.h"
+#include "GFace.h"
+#include "GEdge.h"
+#include "MVertex.h"
+#include "fullMatrix.h"
+
+OptHomPeriodicity::OptHomPeriodicity(std::vector<GEntity*> &entities)
+{
+  for (unsigned int i = 0; i < entities.size(); ++i) {
+    // MVertex on GVertex cannot move
+    if (entities[i]->dim() == 0) continue;
+
+    GEntity *master = entities[i]->meshMaster();
+    if (master != entities[i]) {
+      _master2slave.insert(std::make_pair(master, entities[i]));
+    }
+  }
+}
+
+void OptHomPeriodicity::fixPeriodicity()
+{
+  _relocateMasterVertices();
+  _copyBackMasterVertices();
+}
+
+void OptHomPeriodicity::_relocateMasterVertices()
+{
+  std::multimap<GEntity*, GEntity*>::iterator it;
+  for (it = _master2slave.begin(); it != _master2slave.end(); ++it) {
+    if (it->first->dim() == 2) {
+      GFace *master = dynamic_cast<GFace*>(it->first);
+      GFace *slave = dynamic_cast<GFace*>(it->second);
+      std::vector<double> tfo = _inverse(slave->affineTransform);
+
+      std::map<MVertex*, MVertex*> &vertS2M = slave->correspondingVertices;
+      std::map<MVertex*, MVertex*>::iterator vit;
+      for (vit = vertS2M.begin(); vit != vertS2M.end(); ++vit) {
+        GPoint p = _transform(vit->first, master, tfo);
+        MFaceVertex *v = dynamic_cast<MFaceVertex*>(vit->second);
+        SPoint3 p3 ((v->x()+p.x())/2, (v->y()+p.y())/2, (v->z()+p.z())/2);
+        SPoint2 p2 = master->parFromPoint(p3);
+        GPoint gp = master->point(p2);
+        v->setXYZ(gp.x(), gp.y(), gp.z());
+        v->setParameter(0, gp.u());
+        v->setParameter(1, gp.v());
+      }
+    }
+    else {
+      GEdge *master = dynamic_cast<GEdge*>(it->first);
+      int numSlave = _master2slave.count(master);
+
+      for (int i = 0; i < numSlave; ++i) {
+        if (i > 0) ++it;
+        GEntity *slave = it->second;
+        std::vector<double> tfo = _inverse(slave->affineTransform);
+
+        std::map<MVertex*, MVertex*> &vertS2M = slave->correspondingVertices;
+        std::map<MVertex*, MVertex*>::iterator vit;
+        for (vit = vertS2M.begin(); vit != vertS2M.end(); ++vit) {
+          GPoint p = _transform(vit->first, master, tfo);
+          MVertex *v = vit->second;
+          v->setXYZ(v->x()+p.x(), v->y()+p.y(), v->z()+p.z());
+        }
+      }
+
+      double coeff = 1./(1+numSlave);
+      for (unsigned int k = 0; k < master->getNumMeshVertices(); ++k) {
+        MEdgeVertex *v = dynamic_cast<MEdgeVertex*>(master->getMeshVertex(k));
+        SPoint3 p3 (v->x()*coeff, v->y()*coeff, v->z()*coeff);
+        double u = master->parFromPoint(p3);
+        GPoint gp = master->point(u);
+        v->setXYZ(gp.x(), gp.y(), gp.z());
+        v->setParameter(0, gp.u());
+      }
+    }
+  }
+}
+
+void OptHomPeriodicity::_copyBackMasterVertices()
+{
+  std::multimap<GEntity*, GEntity*>::iterator it;
+  for (it = _master2slave.begin(); it != _master2slave.end(); ++it) {
+    if (it->first->dim() == 2) {
+      GFace *master = dynamic_cast<GFace*>(it->first);
+      GFace *slave = dynamic_cast<GFace*>(it->second);
+      std::vector<double> tfo = slave->affineTransform;
+
+      std::map<MVertex*, MVertex*> &vertS2M = slave->correspondingVertices;
+      std::map<MVertex*, MVertex*>::iterator vit;
+      for (vit = vertS2M.begin(); vit != vertS2M.end(); ++vit) {
+        GPoint p = _transform(vit->second, slave, tfo);
+        MFaceVertex *v = dynamic_cast<MFaceVertex*>(vit->first);
+        v->setXYZ(p.x(), p.y(), p.z());
+        v->setParameter(0, p.u());
+        v->setParameter(1, p.v());
+      }
+    }
+    else {
+      GEdge *master = dynamic_cast<GEdge*>(it->first);
+      GEdge *slave = dynamic_cast<GEdge*>(it->second);
+      std::vector<double> tfo = slave->affineTransform;
+
+      std::map<MVertex*, MVertex*> &vertS2M = slave->correspondingVertices;
+      std::map<MVertex*, MVertex*>::iterator vit;
+      for (vit = vertS2M.begin(); vit != vertS2M.end(); ++vit) {
+        GPoint p = _transform(vit->second, slave, tfo);
+        MEdgeVertex *v = dynamic_cast<MEdgeVertex*>(vit->first);
+        v->setXYZ(p.x(), p.y(), p.z());
+        v->setParameter(0, p.u());
+      }
+    }
+  }
+}
+
+GPoint OptHomPeriodicity::_transform(MVertex *vsource, GEntity *target, std::vector<double> &tfo)
+{
+  double ps[4] = {vsource->x(), vsource->y(), vsource->z(), 1.};
+  double res[4] = {0., 0., 0., 0.};
+  int idx = 0;
+  for(int i = 0; i < 4; i++)
+    for(int j = 0; j < 4; j++)
+      res[i] +=  tfo[idx++] * ps[j];
+
+  SPoint3 p3 (res[0], res[1], res[2]);
+  if (target->dim() == 2) {
+    SPoint2 p2 = dynamic_cast<GFace*>(target)->parFromPoint(p3);
+    return dynamic_cast<GFace*>(target)->point(p2);
+  }
+  else if (target->dim() == 1) {
+    double u = dynamic_cast<GEdge*>(target)->parFromPoint(p3);
+    return dynamic_cast<GEdge*>(target)->point(u);
+  }
+  else {
+    Msg::Error("Expected a face or an edge for computing "
+               "periodicity transformation.");
+    return GPoint();
+  }
+}
+
+std::vector<double> OptHomPeriodicity::_inverse(std::vector<double> &tfo)
+{
+  // Note that the last row of tfo must be (0 0 0 1)...
+  std::vector<double> result(16);
+  fullMatrix<double> mat(4, 4), inv;
+  int idx = 0;
+  for(int i = 0; i < 4; i++)
+    for(int j = 0; j < 4; j++)
+      mat(i, j) = tfo[idx++];
+  mat.invert(inv);
+
+  idx = 0;
+  for(int i = 0; i < 4; i++)
+    for(int j = 0; j < 4; j++)
+      result[idx++] = inv(i, j);
+  return result;
+}

contrib/HighOrderMeshOptimizer/OptHomPeriodicity.h

+// Copyright (C) 2013 ULg-UCL
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use, copy,
+// modify, merge, publish, distribute, and/or sell copies of the
+// Software, and to permit persons to whom the Software is furnished
+// to do so, provided that the above copyright notice(s) and this
+// permission notice appear in all copies of the Software and that
+// both the above copyright notice(s) and this permission notice
+// appear in supporting documentation.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE
+// COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR
+// ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY
+// DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
+// WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
+// ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
+// OF THIS SOFTWARE.
+//
+// Please report all bugs and problems to the public mailing list
+// <gmsh@onelab.info>.
+//
+// Contributors: Amaury Johnen
+
+#ifndef _OPTHOMPERIODICITY_H_
+#define _OPTHOMPERIODICITY_H_
+
+#include <map>
+#include <vector>
+
+class GEntity;
+class GFace;
+class MVertex;
+class GPoint;
+
+class OptHomPeriodicity {
+private:
+  std::multimap<GEntity*, GEntity*> _master2slave;
+  //std::map<GEntity*, std::map<MVertex*, MVertex*> > _ent2vv;
+
+public:
+  OptHomPeriodicity(std::vector<GEntity*>&);
+
+  void fixPeriodicity();
+
+private:
+  void _relocateMasterVertices();
+  void _copyBackMasterVertices();
+
+  static GPoint _transform(MVertex*, GEntity*, std::vector<double>&);
+  static std::vector<double> _inverse(std::vector<double>&);
+};
+
+#endif

contrib/HighOrderMeshOptimizer/OptHomRun.cpp

@@ -34,6 +34,7 @@
 #include "GmshConfig.h"
 #include "OptHOM.h"
 #include "OptHomRun.h"
+#include "OptHomPeriodicity.h"
 #include "GModel.h"
 #include "Gmsh.h"
 #include "MTriangle.h"
@@ -437,6 +438,7 @@ static void optimizeConnectedBlobs(const vertElVecMap &vertex2elements,
   std::vector<elSetVertSetPair> toOptimize =
                           getConnectedBlobs(vertex2elements, badasses, p.nbLayers,
                                     p.distanceFactor, weakMerge, p.optPrimSurfMesh);
+  p.numBlobs = static_cast<int>(toOptimize.size());
 
   //#pragma omp parallel for schedule(dynamic, 1)
   for (int i = 0; i < toOptimize.size(); ++i) {
@@ -507,6 +509,7 @@ static void optimizeOneByOne
 
   const int initNumBadElts = badElts.size();
   Msg::Info("%d badasses, starting to iterate...", initNumBadElts);
+  p.numBlobs = initNumBadElts;
 
   elElSetMap element2elements;                                          // Element to element connectivity, built progressively
 
@@ -630,7 +633,6 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
 
   std::map<MVertex*, std::vector<MElement *> > vertex2elements;
   std::map<MElement*,GEntity*> element2entity;
-  elSet badasses;
   double maxdist = 0.;                                                  // TODO: To be cleaned?
 
   std::map<MElement*,double> distances;
@@ -643,33 +645,37 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
               entity->tag());
     calcVertex2Elements(p.dim,entity,vertex2elements);
     if (p.optPrimSurfMesh) calcElement2Entity(entity,element2entity);
-    for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) {       // Detect bad elements
-      double jmin, jmax;
+  }
+
+  OptHomPeriodicity periodicity = OptHomPeriodicity(entities);
+  do {
+    // Detect bad elements
+    elSet badasses;
+    for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
+      GEntity* &entity = entities[iEnt];
+      if (entity->dim() != p.dim || (p.onlyVisible && !entity->getVisibility())) continue;
+      for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) {
         MElement *el = entity->getMeshElement(iEl);
-      if (el->getDim() == p.dim) {
-	// FIXME TEST
-	//        badasses.insert(el);
         if (p.optCAD) {
-	  //          const double DISTE =computeBndDist(el,2,fabs(p.discrTolerance));
           const double DISTE =distances[el];
-	  //	  if (DISTE > 0)printf("El %d dist %12.5E vs %12.5E\n",iEl,DISTE,p.optCADDistMax);
           maxdist = std::max(DISTE, maxdist);
           if (DISTE > p.optCADDistMax) badasses.insert(el);
         }
+        double jmin, jmax;
         el->scaledJacRange(jmin, jmax, p.optPrimSurfMesh ? entity : 0);
         if (p.BARRIER_MIN_METRIC > 0) jmax = jmin;
         if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX) badasses.insert(el);
       }
     }
-  }
     printf("maxdist = %g badasses size = %lu\n", maxdist, badasses.size());
-
     if (p.strategy == 0)
       optimizeConnectedBlobs(vertex2elements, element2entity, badasses, p, samples, false);
     else if (p.strategy == 2)
       optimizeConnectedBlobs(vertex2elements, element2entity, badasses, p, samples, true);
     else
       optimizeOneByOne(vertex2elements, element2entity, badasses, p, samples);
+    if (p.numBlobs) periodicity.fixPeriodicity();
+  } while (p.numBlobs);
 
   if (p.SUCCESS == 1)
     Msg::Info("Optimization succeeded");

contrib/HighOrderMeshOptimizer/OptHomRun.h

@@ -60,6 +60,7 @@ struct OptHomParameters {
 
   // OUTPUT ------>
   int SUCCESS; // 0 --> success , 1 --> Not converged
+  int numBlobs;
   double minJac, maxJac; // after optimization, range of jacobians
   double CPU; // Time for optimization