Added fast curving functionality to high-order tools

CMakeLists.txt

@@ -133,9 +133,9 @@ set(GMSH_API
   contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h 
     contrib/kbipack/gmp_blas.h contrib/kbipack/mpz.h
   contrib/DiscreteIntegration/Integration3D.h
-  contrib/HighOrderMeshOptimizer/OptHOM.h
-  contrib/HighOrderMeshOptimizer/OptHomMesh.h
-  contrib/HighOrderMeshOptimizer/ParamCoord.h
+  contrib/HighOrderMeshOptimizer/OptHOM.h contrib/HighOrderMeshOptimizer/OptHomMesh.h
+  contrib/HighOrderMeshOptimizer/OptHOMRun.h contrib/HighOrderMeshOptimizer/ParamCoord.h
+  contrib/HighOrderMeshOptimizer/OptHOMFastCurving.h contrib/HighOrderMeshOptimizer/SuperEl.h
   contrib/MathEx/mathex.h)
 
 execute_process(COMMAND date "+%Y%m%d" OUTPUT_VARIABLE DATE 

Fltk/highOrderToolsWindow.cpp

@@ -28,6 +28,7 @@
 #if defined(HAVE_OPTHOM)
 #include "OptHomRun.h"
 #include "OptHomElastic.h"
+#include "OptHomFastCurving.h"
 #endif
 
 static void change_completeness_cb(Fl_Widget *w, void *data)
@@ -75,25 +76,36 @@ static void highordertools_runp_cb(Fl_Widget *w, void *data)
 static void chooseopti_cb(Fl_Widget *w, void *data)
 {
   highOrderToolsWindow *o = FlGui::instance()->highordertools;
-  int elastic = o->choice[2]->value();
+  const int algo = o->choice[2]->value();
 
-  if (elastic == 1){
+  switch(algo) {
+  case 0: {                                           // Optimization
+    o->push[0]->activate();
+    o->choice[0]->activate();
+    o->choice[3]->activate();
+    for (int i=2;i<=11;i++) o->value[i]->activate();
+    //    o->push[1]->activate();
+    break;
+  }
+  case 1: {                                           // Elastic analogy
     o->choice[0]->deactivate();
     o->choice[3]->deactivate();
-    for (int i=3;i<=6;i++)
+    for (int i=2;i<=11;i++) o->value[i]->deactivate();
+    //   o->push[1]->deactivate();
+    break;
+  }
+  case 2: {                                           // Fast curving
+    o->choice[0]->deactivate();
+    o->choice[3]->deactivate();
+    for (int i=2;i<=6;i++)
       o->value[i]->deactivate();
     for (int i=9;i<=11;i++) o->value[i]->deactivate();
+    for (int i=7;i<=8;i++) o->value[i]->activate();
     //   o->push[1]->deactivate();
+    break;
   }
-  else {
-    o->push[0]->activate();
-    o->choice[0]->activate();
-    o->choice[3]->activate();
-    for (int i=3;i<=6;i++)
-      o->value[i]->activate();
-    for (int i=9;i<=11;i++) o->value[i]->activate();
-    //    o->push[1]->activate();
   }
+
 }
 
 static void chooseopti_strategy(Fl_Widget *w, void *data)
@@ -110,17 +122,15 @@ static void highordertools_runopti_cb(Fl_Widget *w, void *data)
   if(o->butt[3]->value())
     FlGui::instance()->graph[0]->showMessages();
 
-  bool elastic = o->choice[2]->value() == 1;
+  const int algo = o->choice[2]->value();
   double threshold_min = o->value[1]->value();
   bool onlyVisible = (bool)o->butt[1]->value();
   int nbLayers = (int) o->value[2]->value();
   double threshold_max = o->value[8]->value();
 
 #if defined(HAVE_OPTHOM)
-  if(elastic){
-    ElasticAnalogy(GModel::current(), threshold_min, onlyVisible);
-  }
-  else{
+  switch(algo) {
+  case 0: {                                                               // Optimization
     OptHomParameters p;
     p.nbLayers = nbLayers;
     p.BARRIER_MIN = threshold_min;
@@ -138,6 +148,22 @@ static void highordertools_runopti_cb(Fl_Widget *w, void *data)
     p.adaptBlobLayerFact = (int) o->value[10]->value();
     p.adaptBlobDistFact = o->value[11]->value();
     HighOrderMeshOptimizer(GModel::current(), p);
+    break;
+  }
+  case 1: {                                                               // Elastic analogy
+    ElasticAnalogy(GModel::current(), threshold_min, onlyVisible);
+    break;
+  }
+  case 2: {                                                               // Fast curving
+    FastCurvingParameters p;
+    p.BARRIER_MIN = threshold_min;
+    p.BARRIER_MAX = threshold_max;
+    p.onlyVisible = onlyVisible;
+    p.dim = GModel::current()->getDim();
+    p.distanceFactor =  o->value[7]->value();
+    HighOrderMeshFastCurving(GModel::current(), p);
+    break;
+  }
   }
 #else
   Msg::Error("High-order mesh optimization requires the OPTHOM module");
@@ -276,7 +302,8 @@ highOrderToolsWindow::highOrderToolsWindow(int deltaFontSize)
 
   static Fl_Menu_Item menu_method[] = {
     {"Optimization", 0, 0, 0},
-    {"ElasticAnalogy", 0, 0, 0},
+    {"Elastic Analogy", 0, 0, 0},
+    {"Fast Curving", 0, 0, 0},
     {0}
   };
 

contrib/HighOrderMeshOptimizer/CMakeLists.txt

@@ -8,7 +8,9 @@ set(SRC
   OptHOM.cpp 
   OptHomRun.cpp 
   ParamCoord.cpp 
+  SuperEl.cpp 
   OptHomElastic.cpp
+  OptHomFastCurving.cpp
 )
 
 file(GLOB_RECURSE HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.hpp)

contrib/HighOrderMeshOptimizer/OptHomFastCurving.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@geuz.org>.
+//
+// Contributors: Thomas Toulorge, Jonathan Lambrechts
+
+#include <stdio.h>
+#include <sstream>
+#include <fstream>
+#include <iterator>
+#include <string.h>
+#include "GmshConfig.h"
+#include "OptHomFastCurving.h"
+#include "GModel.h"
+#include "Gmsh.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MTetrahedron.h"
+#include "MHexahedron.h"
+#include "MPrism.h"
+#include "MLine.h"
+#include "OS.h"
+#include <stack>
+#include "SuperEl.h"
+#include "SVector3.h"
+#include "BasisFactory.h"
+
+
+
+namespace {
+
+
+
+void exportMeshToDassault(GModel *gm, const std::string &fn, int dim)
+{
+  FILE *f = fopen(fn.c_str(),"w");
+
+  int numVertices = gm->indexMeshVertices(true);
+  std::vector<GEntity*> entities;
+  gm->getEntities(entities);
+  fprintf(f,"%d %d\n", numVertices, dim);
+  for(unsigned int i = 0; i < entities.size(); i++)
+    for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
+      MVertex *v = entities[i]->mesh_vertices[j];
+      if (dim == 2)
+        fprintf(f,"%d %22.15E %22.15E\n", v->getIndex(), v->x(), v->y());
+      else if (dim == 3)
+        fprintf(f,"%d %22.15E %22.15E %22.5E\n", v->getIndex(), v->x(),
+                v->y(), v->z());
+    }
+
+  if (dim == 2){
+    int nt = 0;
+    int order  = 0;
+    for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf){
+      std::vector<MTriangle*> &tris = (*itf)->triangles;
+      nt += tris.size();
+      if (tris.size())order = tris[0]->getPolynomialOrder();
+    }
+    fprintf(f,"%d %d\n", nt,(order+1)*(order+2)/2);
+    int count = 1;
+    for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf){
+      std::vector<MTriangle*> &tris = (*itf)->triangles;
+      for (size_t i=0;i<tris.size();i++){
+  MTriangle *t = tris[i];
+  fprintf(f,"%d ", count++);
+  for (int j=0;j<t->getNumVertices();j++){
+    fprintf(f,"%d ", t->getVertex(j)->getIndex());
+  }
+  fprintf(f,"\n");
+      }
+    }
+    int ne = 0;
+    for (GModel::eiter ite = gm->firstEdge(); ite != gm->lastEdge(); ++ite){
+      std::vector<MLine*> &l = (*ite)->lines;
+      ne += l.size();
+    }
+    fprintf(f,"%d %d\n", ne,(order+1));
+    count = 1;
+    for (GModel::eiter ite = gm->firstEdge(); ite != gm->lastEdge(); ++ite){
+      std::vector<MLine*> &l = (*ite)->lines;
+      for (size_t i=0;i<l.size();i++){
+  MLine *t = l[i];
+  fprintf(f,"%d ", count++);
+  for (int j=0;j<t->getNumVertices();j++){
+    fprintf(f,"%d ", t->getVertex(j)->getIndex());
+  }
+  fprintf(f,"%d \n",(*ite)->tag());
+      }
+    }
+  }
+  fclose(f);
+}
+
+
+
+//double distMaxStraight(MElement *el)
+//{
+//  const polynomialBasis *lagrange = (polynomialBasis*)el->getFunctionSpace();
+//  const polynomialBasis *lagrange1 = (polynomialBasis*)el->getFunctionSpace(1);
+//  int nV = lagrange->points.size1();
+//  int nV1 = lagrange1->points.size1();
+//  SPoint3 sxyz[256];
+//  for (int i = 0; i < nV1; ++i) {
+//    sxyz[i] = el->getVertex(i)->point();
+//  }
+//  for (int i = nV1; i < nV; ++i) {
+//    double f[256];
+//    lagrange1->f(lagrange->points(i, 0), lagrange->points(i, 1),
+//                 lagrange->points(i, 2), f);
+//    for (int j = 0; j < nV1; ++j)
+//      sxyz[i] += sxyz[j] * f[j];
+//  }
+//
+//  double maxdx = 0.0;
+//  for (int iV = nV1; iV < nV; iV++) {
+//    SVector3 d = el->getVertex(iV)->point()-sxyz[iV];
+//    double dx = d.norm();
+//    if (dx > maxdx) maxdx = dx;
+//  }
+//
+//  return maxdx;
+//}
+
+
+
+void calcVertex2Elements(int dim, GEntity *entity,
+                                std::map<MVertex*, std::vector<MElement *> > &vertex2elements)
+{
+  for (size_t i = 0; i < entity->getNumMeshElements(); ++i) {
+    MElement *element = entity->getMeshElement(i);
+    if (element->getDim() == dim)
+      for (int j = 0; j < element->getNumPrimaryVertices(); ++j)
+        vertex2elements[element->getVertex(j)].push_back(element);
+  }
+}
+
+
+
+//void getNormals2D(const std::vector<MVertex*> &edgeVert, SVector3 &n0, SVector3 &n1) {
+//
+////  if (v0->onWhat() != 0) {
+//  if (false) {
+//    GEdge *ent = edgeVert.front()->onWhat()->cast2Edge();
+//    std::map<MVertex*, SVector3, std::less<MVertex*> > &allNormals = ent->getNormals();
+//    n0 = allNormals[edgeVert[0]];
+//    n1 = allNormals[edgeVert[1]];
+//  }
+//  else {
+//    const int &Nsf = edgeVert.size();
+//    const nodalBasis* fs = BasisFactory::getNodalBasis(nodalBasis::getTag(TYPE_LIN,Nsf-1,false));
+//    double dsf0[Nsf][3], dsf1[Nsf][3];
+//    fs->df(fs->points(0,0),0.,0.,dsf0);
+//    fs->df(fs->points(1,0),0.,0.,dsf1);
+//    double t0x = 0., t0y = 0., t1x = 0., t1y = 0.;
+//    for (int i=0; i<Nsf; ++i) {
+//      t0x += dsf0[i][0]*edgeVert[i]->x();
+//      t0y += dsf0[i][0]*edgeVert[i]->y();
+//      t1x += dsf1[i][0]*edgeVert[i]->x();
+//      t1y += dsf1[i][0]*edgeVert[i]->y();
+//    }
+//    n0[0] = -t0y; n0[1] = t0x; n0[2] = 0.;
+//    n0.normalize();
+//    n1[0] = -t1y; n1[1] = t1x; n1[2] = 0.;
+//    n1.normalize();
+//  }
+//
+//}
+
+
+
+//void getBaseVertsAndNormals(std::set<MElement*> badElts, std::map<MElement*,std::vector<MVertex*> > &faces,
+//                                 std::map<MVertex*,SVector3> &normals) {
+//
+//  typedef std::map<MVertex*,std::vector<SVector3> > normalMap;
+//  normalMap allNormals;
+//
+//  for (std::set<MElement*>::iterator itBE = badElts.begin(); itBE != badElts.end(); itBE++) {
+//    std::vector<MVertex*> baseVert;
+//    for (int i=0; i<(*itBE)->getNumEdges(); i++) {
+//      (*itBE)->getEdgeVertices(i,baseVert);
+//      if (isEdgeCurved(baseVert)) {
+//        const SPoint3 p0 = baseVert[0]->point(), p1 = baseVert[1]->point(), pC = (*itBE)->barycenter(true);
+//        SVector3 n(p0.y()-p1.y(), p1.x()-p0.x(), 0.), p0C(p0,pC);
+//        if (dot(n,p0C) < 0.) n *= -1.;
+//        n.normalize();
+//        allNormals[baseVert[0]].push_back(n);
+//        allNormals[baseVert[1]].push_back(n);
+//        faces[*itBE] = baseVert;
+//        break;
+//      }
+//    }
+//  }
+//
+//  for(normalMap::iterator itAN = allNormals.begin(); itAN != allNormals.end(); itAN++) {
+//    std::vector<SVector3> &normSet = itAN->second;
+////    std::cout << "DBGTT: Considering node " << itAN->first->getNum() << "\n";
+//    SVector3 norm(0.);
+//    for (std::vector<SVector3>::iterator itN = normSet.begin(); itN != normSet.end(); itN++) {
+////      std::cout << "DBGTT:  -> accumulating normal (" << itN->x() << ", "<< itN->y() << ", "<< itN->z() << ")\n";
+//      norm += *itN;
+//    }
+//    norm *= 1./double(normSet.size());
+//    normals[itAN->first] = norm;
+//  }
+//
+//}
+
+
+
+//bool isEdgeCurved(const std::vector<MVertex*> &edgeVert) {
+//
+//  static const double eps = 1.e-8;
+//
+//  const double n = edgeVert.size();
+//  const MVertex *vS = edgeVert[0];
+//  const double vSEx = edgeVert[1]->x()-vS->x(), vSEy = edgeVert[1]->y()-vS->y();
+//  const double normSE = vSEx*vSEx+vSEy*vSEy;
+//
+//  double distSq = 0.;
+//  for (int i=2; i<n; ++i) {
+//    const double vSix = edgeVert[i]->x()-vS->x(), vSiy = edgeVert[i]->y()-vS->y();
+//    const double fact = (vSix*vSEx+vSiy*vSEy)/normSE;
+//    const double dvx = vSix-fact*vSEx, dvy = vSiy-fact*vSEy;
+//    distSq += dvx*dvx+dvy*dvy;
+//  }
+//
+//  return (distSq > eps*normSE);
+//
+//}
+
+
+
+// Among edges connected to a given vertex, return the direction of the one that is closest to the given normal
+// Return the given normal if no edge is sufficiently close
+SVector3 getNormalEdge(MVertex *vert, const SVector3 &n,
+                       const std::map<MVertex*, std::vector<MElement*> > &vertex2elements) {
+
+//  static const double spLimit = 0.70711;                          // Limit in dot product below which we return the normal
+  static const double spLimit = 0.5;                          // Limit in dot product below which we return the normal
+
+  const std::vector<MElement*> &elts = (*vertex2elements.find(vert)).second;     // All elements connected to vertex
+
+  double spMax = 0.;
+  SVector3 normalEdge;
+//  std::cout << "DBGTT: Looking for normal edge at vertex " << vert->getNum()
+//              << " with n = (" << n.x() << ", " << n.y() << ", " << n.z() << ")\n";
+
+  for (std::vector<MElement*>::const_iterator itEl = elts.begin(); itEl != elts.end(); ++itEl)
+    for (int i=0; i<(*itEl)->getNumEdges(); i++) {
+      std::vector<MVertex*> edgeVert;
+      (*itEl)->getEdgeVertices(i,edgeVert);
+      SVector3 edge;
+      if (edgeVert[0] == vert) edge = SVector3(vert->point(),edgeVert[1]->point());
+      else if (edgeVert[1] == vert) edge = SVector3(vert->point(),edgeVert[0]->point());
+      else continue;
+      edge.normalize();
+      double sp = dot(edge,n);
+      if (sp > spMax) {                                           // Retain the edge giving max. dot product with normal
+        spMax = sp;
+        normalEdge = edge;
+      }
+//      std::cout << "DBGTT:   -> checking edge " << edgeVert[0]->getNum() << " - " << edgeVert[1]->getNum()
+//                  << ", sp = " << sp << "\n";
+    }
+
+  if (spMax < spLimit) { std::cout << "DBGTT: no normal edge\n"; normalEdge = n; }                            // If max. dot product is below limit, just take normal
+
+  return normalEdge;
+
+}
+
+
+
+//// Detect whether edge/face is curved, and give (non-oriented) normal
+//bool isCurvedAndNormal(int type, int order, const std::vector<MVertex*> &faceVert, SVector3 &normal) {
+//
+//  static const double eps = 1.e-8;
+//
+//  const nodalBasis *lagrange = BasisFactory::getNodalBasis(ElementType::getTag(type,order,false));
+//  const nodalBasis *lagrange1 = BasisFactory::getNodalBasis(ElementType::getTag(type,1,false));
+//  int nV = lagrange->points.size1();
+//  int nV1 = lagrange1->points.size1();
+//  SPoint3 sxyz[256];
+//  for (int i = 0; i < nV1; ++i) sxyz[i] = faceVert[i]->point();
+//  for (int i = nV1; i < nV; ++i) {
+//    double f[256];
+//    lagrange1->f(lagrange->points(i, 0), lagrange->points(i, 1), lagrange->points(i, 2), f);
+//    for (int j = 0; j < nV1; ++j)
+//      sxyz[i] += sxyz[j] * f[j];
+//  }
+//
+//  double sumDistSq = 0.;
+//  for (int iV = nV1; iV < nV; iV++) sumDistSq += SVector3(sxyz[iV],faceVert[iV]->point()).normSq();
+//
+//  double scale;
+//  if (type == TYPE_LIN) {
+//    const SPoint3 &p0 = sxyz[0], &p1 = sxyz[1];
+//    normal = SVector3(p0.y()-p1.y(),p1.x()-p0.x(),0.);
+//    scale = normal.normSq();
+//  }
+//  else {
+//    const SPoint3 &p0 = sxyz[0], &p1 = sxyz[1], &p2 = sxyz[2];
+//    SVector3 p01(p0,p1), p02(p0,p2);
+//    normal = crossprod(p01,p02);
+//    scale = normal.norm();
+//  }
+//  normal.normalize();
+//
+////  std::cout << "DBGTT: v0 is " << faceVert[0]->getNum() << ", v1 is " << faceVert[1]->getNum() << ", sumDistSq = " << sumDistSq << ", scale = " << scale << ", test = " << (sumDistSq > eps*double(nV)*scale) << "\n";
+//  return (sumDistSq > eps*double(nV)*scale);
+//
+//}
+//
+//
+//
+//void getBaseVertsAndNormals(std::set<MElement*> badElts,
+//                              const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
+//                              std::map<MElement*,int> &baseType,
+//                              std::map<MElement*,std::vector<MVertex*> > &baseVert,
+//                              std::map<MElement*,std::vector<SVector3> > &normals) {
+//
+//  for (std::set<MElement*>::iterator itBE = badElts.begin(); itBE != badElts.end(); itBE++) {
+//    const bool is2D = ((*itBE)->getDim() == 2);
+//    const int order = (*itBE)->getPolynomialOrder();
+//    const int numFaces = is2D ? (*itBE)->getNumEdges() : (*itBE)->getNumFaces();
+//    for (int iFace=0; iFace<numFaces; iFace++) {            // Loop on edges/faces
+//      std::vector<MVertex*> bv;                             // Vertices of edge/face
+//      int type, numPrimVert;
+//      if (is2D) {
+//        (*itBE)->getEdgeVertices(iFace,bv);
+//        type = TYPE_LIN;
+//        numPrimVert = 2;
+//      }
+//      else {
+//        (*itBE)->getFaceVertices(iFace,bv);
+//        MFace face = (*itBE)->getFace(iFace);
+//        numPrimVert = face.getNumVertices();
+//        type = (numPrimVert == 3) ? TYPE_TRI : TYPE_QUA;
+//      }
+//      SVector3 n;                                           // Normal to straight edge/face
+//      if (isCurvedAndNormal(type,order,bv,n)) {             // Check whether edge/face is curved and compute straight normal
+//        std::cout << "DBGTT: Curved edge/face in el. " << (*itBE)->getNum() << "\n";
+//        if (baseVert.find(*itBE) != baseVert.end()) {       // If more than 1 curved edge/face in bad el., drop it
+//          baseVert.erase(*itBE);
+//          normals.erase(*itBE);
+//          std::cout << "DBGTT: More than 1 curved edge/face detected in el. " << (*itBE)->getNum() << "\n";
+//          break;
+//        }
+//        baseType[*itBE] = type;
+//        baseVert[*itBE] = bv;
+//        const SPoint3 p0 = bv[0]->point(), pC = (*itBE)->barycenter(true);
+//        SVector3 p0C(p0,pC);
+//        if (dot(n,p0C) < 0.) n *= -1.;                                      // Make straight normal in-going
+//        for (int iV=0; iV<numPrimVert; iV++)                                // Compute normals to prim. vert. of edge/face
+//          normals[*itBE].push_back(getNormalEdge(bv[iV],n,vertex2elements));
+//      }
+//    }
+//  }
+//
+//}
+
+
+
+// Detect whether edge/face is curved, and give normal
+bool isCurvedAndNormal(int type, int order, const std::vector<MVertex*> &faceVert,
+                       const SPoint3 pBar, SVector3 &normal, double &maxDist) {
+
+//  static const double eps = 1.e-10;
+  static const double eps = 1.e-6;
+
+  // Compute HO points in straight edge/face
+  const nodalBasis *lagrange = BasisFactory::getNodalBasis(ElementType::getTag(type,order,false));
+  const nodalBasis *lagrange1 = BasisFactory::getNodalBasis(ElementType::getTag(type,1,false));
+  int nV = lagrange->points.size1();
+  int nV1 = lagrange1->points.size1();
+  SPoint3 sxyz[256];
+  for (int i = 0; i < nV1; ++i) sxyz[i] = faceVert[i]->point();
+  for (int i = nV1; i < nV; ++i) {
+    double f[256];
+    lagrange1->f(lagrange->points(i, 0), lagrange->points(i, 1), lagrange->points(i, 2), f);
+    for (int j = 0; j < nV1; ++j)
+      sxyz[i] += sxyz[j] * f[j];
+  }
+
+  // Compute (non-oriented) unit normal to straight edge/face and its scale [length]
+  double scale;
+  const SPoint3 &p0 = sxyz[0], &p1 = sxyz[1];
+  if (type == TYPE_LIN) {
+    normal = SVector3(p0.y()-p1.y(),p1.x()-p0.x(),0.);
+    scale = normal.normalize();
+  }
+  else {
+    const SPoint3 &p2 = sxyz[2];
+    SVector3 p01(p0,p1), p02(p0,p2);
+    normal = crossprod(p01,p02);
+    scale = sqrt(normal.normalize());
+  }
+
+  // Orient normal with help of pBar
+  SVector3 p0C(p0,pBar);
+  if (dot(normal,p0C) < 0.) normal *= -1.;                                      // Make straight normal in-going
+
+  // Calc max. normal dist. from straight to HO points
+  maxDist = 0.;
+  for (int iV = nV1; iV < nV; iV++) {
+    const double normalDisp = dot(SVector3(sxyz[iV],faceVert[iV]->point()),normal);
+    maxDist = std::max(maxDist,fabs(normalDisp));
+  }
+
+  std::cout << "DBGTT: v0 is " << faceVert[0]->getNum() << ", v1 is " << faceVert[1]->getNum()
+            << ", v2 is " << faceVert[2]->getNum() << ", maxDist = " << maxDist
+            << ", scale = " << scale << ", test = " << (maxDist > eps*scale) << "\n";
+  return (maxDist > eps*scale);
+
+}
+
+
+
+bool getCurvedFace(MElement* badEl,
+                   const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
+                   std::vector<MVertex*> &baseVert, std::vector<SVector3> &normals,
+                   int &baseType, double &baseMaxDist) {
+
+  bool found = false;
+  const int dim = badEl->getDim();
+  const int order = badEl->getPolynomialOrder();
+  const int numFaces = (dim == 2) ? badEl->getNumEdges() : badEl->getNumFaces();
+
+  for (int iFace=0; iFace<numFaces; iFace++) {            // Loop on edges/faces
+
+    // Get face info (type, vertices)
+    std::vector<MVertex*> bv;                             // Vertices of edge/face
+    int numPrimVert, type;
+    if (dim == 2) {
+      badEl->getEdgeVertices(iFace,bv);
+      type = TYPE_LIN;
+      numPrimVert = 2;
+    }
+    else {
+      badEl->getFaceVertices(iFace,bv);
+      MFace face = badEl->getFace(iFace);
+      numPrimVert = face.getNumVertices();
+      type = (numPrimVert == 3) ? TYPE_TRI : TYPE_QUA;
+    }
+
+    // Skip face if at least 1 vertex not on boundary
+    bool inDom = false;
+    for (int iV=0; iV<bv.size(); ++iV)
+      if (bv[iV]->onWhat()->dim() == dim) inDom = true;
+    if (inDom) continue;
+
+    //    std::cout << "DBGTT: Checking edge/face " << iFace << " in el. " << badEl->getNum() << ": "
+//              << numPrimVert << " vert., type " << type << "\n";
+    // If face curved, mark it and compute normals to its primary vertices
+    SVector3 n;                                                           // Normal to straight edge/face
+    double maxDist;                                                       // TOFIX: Max. normal. dist. to straight in curved face
+    const SPoint3 pBar = badEl->barycenter(true);                         // Bary. of el. to help orienting normal
+    if (isCurvedAndNormal(type,order,bv,pBar,n,maxDist)) {                // Check whether edge/face is curved and compute straight normal
+      std::cout << "DBGTT: Curved edge/face in el. " << badEl->getNum() << "\n";
+      if (found) {                                                        // If more than 1 curved edge/face in bad el., drop it
+        std::cout << "DBGTT: More than 1 curved edge/face detected in el. " << badEl->getNum() << "\n";
+        return false;
+      }
+      for (int iV=0; iV<numPrimVert; iV++)                                // Compute normals to prim. vert. of edge/face
+        normals.push_back(getNormalEdge(bv[iV],n,vertex2elements));
+      baseVert = bv;
+      baseType = type;
+      baseMaxDist = maxDist;
+      found = true;
+    }
+
+  }
+
+  return found;
+
+}
+
+
+
+void makeStraight(MElement *el, const std::set<MVertex*> &movedVert) {
+
+  const nodalBasis *nb = el->getFunctionSpace();
+  const fullMatrix<double> &pts = nb->points;
+
+  SPoint3 p;
+
+  for(int iPt = el->getNumPrimaryVertices(); iPt < el->getNumVertices(); ++iPt) {
+    MVertex *vert = el->getVertex(iPt);
+    if (movedVert.find(vert) == movedVert.end()) {
+      el->primaryPnt(pts(iPt,0),pts(iPt,1),pts(iPt,2),p);
+      vert->setXYZ(p.x(),p.y(),p.z());
+    }
+  }
+
+}
+
+
+
+std::set<MElement*> getSuperElBlob(MElement *el, const std::map<MVertex*,
+                                   std::vector<MElement*> > &vertex2elements,
+                                   const SuperEl *sEl)
+{
+
+  static const int depth = 100;
+
+  std::set<MElement*> blob;
+  std::list<MElement*> currentLayer, lastLayer;
+
+  blob.insert(el);
+  lastLayer.push_back(el);
+  for (int d = 0; d < depth; ++d) {
+    currentLayer.clear();
+    for (std::list<MElement*>::iterator it = lastLayer.begin();
+         it != lastLayer.end(); ++it) {
+      for (int i = 0; i < (*it)->getNumPrimaryVertices(); ++i) {
+        const std::vector<MElement*> &neighbours = vertex2elements.find
+          ((*it)->getVertex(i))->second;
+        for (std::vector<MElement*>::const_iterator itN = neighbours.begin();
+             itN != neighbours.end(); ++itN){
+          if (sEl->isPointIn((*itN)->barycenter(true))) {
+            // Assume that if an el is too far, its neighbours are too far as well
+            if (blob.insert(*itN).second) currentLayer.push_back(*itN);
+          }
+        }
+      }
+    }
+    lastLayer = currentLayer;
+    if (currentLayer.empty()) break;
+  }
+
+  return blob;
+}
+
+
+
+void curveMesh(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
+               std::set<MElement*> &badElements, FastCurvingParameters &p, int samples)
+{
+
+  const int nbBadElts = badElements.size();
+  std::vector<MElement*> badElts;
+//  std::vector<int> types;
+//  std::vector<std::vector<MVertex*> > faces;
+//  std::vector<std::vector<SVector3> > &normals;
+//  std::vector<double> &maxDist;
+  std::vector<SuperEl*> superElts;
+  badElts.reserve(nbBadElts);
+//  types.reserve(nbBadElts);
+//  faces.reserve(nbBadElts);
+//  normals.reserve(nbBadElts);
+//  maxDist.reserve(nbBadElts);
+  superElts.reserve(nbBadElts);
+
+  std::ofstream of("dum.pos");
+  of << "View \" \"{\n";
+
+  for (std::set<MElement*>::const_iterator itBE = badElements.begin(); itBE != badElements.end(); ++itBE) {
+    std::vector<MVertex*> faceBE;
+    std::vector<SVector3> normalsBE;
+    int typeBE;
+    double maxDistBE;
+    if (getCurvedFace(*itBE,vertex2elements,faceBE,normalsBE,typeBE,maxDistBE)) {
+      badElts.push_back(*itBE);
+//      types.push_back(typeBE);
+//      faces.push_back(faceBE);
+//      normals.push_back(normalsBE);
+//      maxDist.push_back(maxDistBE);
+      superElts.push_back(new SuperEl(*itBE,maxDistBE*p.distanceFactor,typeBE,faceBE,normalsBE));
+//      superElts.push_back(new SuperEl(*itBE,distMaxStraight(*itBE)*p.distanceFactor,typeBE,faceBE,normalsBE));
+//      movedVert.insert(faceBE.begin(),faceBE.end());                                                // Do not touch vert. of curved face
+      of << superElts.back()->printPOS();
+    }
+  }
+
+  of << "};\n";
+  of.close();
+
+  std::set<MVertex*> movedVert;
+  for (int iBE=0; iBE<badElts.size(); ++iBE) {
+//    std::set<MElement*> blob = getSurroundingBlob(badElts[iBE], p.nbLayers, vertex2elements, p.distanceFactor);
+    std::set<MElement*> blob = getSuperElBlob(badElts[iBE], vertex2elements, superElts[iBE]);
+    std::cout << "DBGTT: Blob of bad el. " << badElts[iBE]->getNum() << " contains elts.";
+    for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE) std::cout << " " << (*itE)->getNum();
+    std::cout << "\n";
+    makeStraight(badElts[iBE],movedVert);                                             // Make bad. el. straight
+    for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE)
+      for (int i = 0; i < (*itE)->getNumVertices(); ++i) {                            // For each vert. of each el. in blob
+        MVertex* vert = (*itE)->getVertex(i);
+        if (movedVert.find(vert) == movedVert.end()) {                                // If vert. not already moved
+          double xyzS[3] = {vert->x(), vert->y(), vert->z()}, xyzC[3];
+          if (superElts[iBE]->straightToCurved(xyzS,xyzC)) {
+            std::cout << "DBGTT: moving vertex " << vert->getNum() << " from (" << xyzS[0] << "," << xyzS[1] << "," << xyzS[2] << ") to (" << xyzC[0] << "," << xyzC[1] << "," << xyzC[2] << ")\n";
+            vert->setXYZ(xyzC[0],xyzC[1],xyzC[2]);
+            movedVert.insert(vert);
+          }
+          else std::cout << "DBGTT: Failed to move vertex " << vert->getNum() << " with bad. el " << badElts[iBE]->getNum() << "\n";
+        }
+//        else std::cout << "DBGTT: Already moved vertex " << vert->getNum() << " with bad. el " << badElts[iBE]->getNum() << "\n";
+      }
+  }
+
+}
+
+
+
+}
+
+
+
+//static void optimizeBL(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
+//                       std::set<MElement*> &badElements, OptHomParameters &p, int samples)
+//{
+//
+//  const int nbBadElts = badElements.size();
+//  std::vector<MElement*> &badElts();
+//  std::map<MElement*,int> types;                                 // Type (LIN, TRI, QUA) of curved face in bad. el.
+//  std::map<MElement*,std::vector<MVertex*> > faces;              // Vertices of curved face in bad. el.
+//  std::map<MElement*,std::vector<SVector3> > normals;            // Normal to each prim. vertex of curved face in bad. el.
+//  getBaseVertsAndNormals(badElements,vertex2elements,types,faces,normals);
+//
+//  std::ofstream of("dum.pos");
+//  of << "View \" \"{\n";
+//
+//  for (std::set<MElement*>::const_iterator itBE = badElements.begin(); itBE != badElements.end(); ++itBE) {
+//
+//    if (faces.find(*itBE) == faces.end()) continue;           // No base edge/face -> no super el. for this bad. el.
+//
+//    SuperEl se(*itBE,distMaxStraight(*itBE)*p.distanceFactor,types[*itBE],faces[*itBE],normals[*itBE]);
+//
+//    std::set<MElement*> blob = getSurroundingBlob(*itBE, p.nbLayers, vertex2elements, p.distanceFactor);
+//    std::set<MVertex*> toMove;
+//    for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE)
+//      for (int i = 0; i < (*itE)->getNumVertices(); ++i) toMove.insert((*itE)->getVertex(i));
+//
+////    std::cout << "DBGTT: Before makeStraight:\n";
+////    se.printCoord();
+//    makeStraight(*itBE);                                      // Make bad. el. straight
+////    for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE) makeStraight(*itE);
+////    std::cout << "DBGTT: After makeStraight:\n";
+////    se.printCoord();
+//
+//    // Apply transformation straight -> curved to all vertices in super el.
+//    for (std::set<MVertex*>::iterator itV = toMove.begin(); itV != toMove.end(); ++itV) {
+//      double xyzS[3] = {(*itV)->x(), (*itV)->y(), (*itV)->z()}, xyzC[3];
+//      if (se.straightToCurved(xyzS,xyzC)) {
+////        std::cout << "DBGTT: moving vertex " << (*itV)->getNum() << " from (" << xyzS[0] << "," << xyzS[1] << "," << xyzS[2] << ") to (" << xyzC[0] << "," << xyzC[1] << "," << xyzC[2] << ")\n";
+//        (*itV)->setXYZ(xyzC[0],xyzC[1],xyzC[2]);
+//      }
+//      else std::cout << "DBGTT: Failed to move vertex " << (*itV)->getNum() << " with bad. el " << (*itBE)->getNum() << "\n";
+//    }
+//
+////    std::cout << se.printPOS();
+//    of << se.printPOS();
+//
+//  }
+//
+//  of << "};\n";
+//  of.close();
+//
+//}
+
+
+
+void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p)
+{
+
+  double t1 = Cpu();
+
+  int samples = 30;
+
+  double tf1 = Cpu();
+
+  Msg::StatusBar(true, "Optimizing high order mesh...");
+  std::vector<GEntity*> entities;
+  gm->getEntities(entities);
+
+  std::map<MVertex*, std::vector<MElement *> > vertex2elements;
+  std::set<MElement*> badasses;
+  for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
+    GEntity* &entity = entities[iEnt];
+    if (entity->dim() != p.dim || (p.onlyVisible && !entity->getVisibility())) continue;
+    Msg::Info("Computing connectivity and bad elements for entity %d...",
+              entity->tag());
+    calcVertex2Elements(p.dim,entity,vertex2elements); // Compute vert. -> elt. connectivity
+//    std::cout << "DBGTT: num. el. = " << entity->getNumMeshElements() << "\n";
+    for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
+      double jmin, jmax;
+      MElement *el = entity->getMeshElement(iEl);
+      if (el->getDim() == p.dim) {
+        el->scaledJacRange(jmin, jmax);
+//        std::cout << "El. " << iEl << ", jmin = " << jmin << ", jmax = " << jmax << "\n";
+//        Msg::Info("El %i, jmin = %g, jmax = %g\n",iEl,jmin,jmax);
+        if (jmin < p.BARRIER_MIN || jmax > p.BARRIER_MAX) badasses.insert(el);
+      }
+    }
+  }
+
+  curveMesh(vertex2elements, badasses, p, samples);
+
+  double t2 = Cpu();
+
+  Msg::StatusBar(true, "Done curving high order mesh (%g s)", t2-t1);
+
+}

contrib/HighOrderMeshOptimizer/OptHomFastCurving.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@geuz.org>.
+//
+// Contributors: Thomas Toulorge, Jonathan Lambrechts
+
+#ifndef _OPTHOMFASTCURVING_H_
+#define _OPTHOMFASTCURVING_H_
+
+class GModel;
+
+struct FastCurvingParameters {
+  // INPUT ------>
+  double BARRIER_MIN ; // minimum scaled jcaobian
+  double BARRIER_MAX ; // maximum scaled jcaobian
+  int dim ; // which dimension to optimize
+  bool onlyVisible ; // apply optimization to visible entities ONLY
+  double distanceFactor; // filter elements such that no elements further away
+                         // than DistanceFactor times the max distance to
+                         // straight sided version of an element are optimized
+
+  FastCurvingParameters ()
+    : BARRIER_MIN(0.1), BARRIER_MAX(2.0), dim(3) , onlyVisible(true), distanceFactor(12)
+  {
+  }
+};
+
+void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p);
+
+#endif

contrib/HighOrderMeshOptimizer/SuperEl.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@geuz.org>.
+//
+// Contributor: Thomas Toulorge
+
+#include <iostream>
+#include <sstream>
+#include "GEdge.h"
+#include "MLine.h"
+#include "MQuadrangle.h"
+#include "MPrism.h"
+#include "MHexahedron.h"
+#include "SuperEl.h"
+
+
+
+SuperEl::SuperEl(MElement *badEl, double dist, int type, const std::vector<MVertex*> &baseVert,
+                 const std::vector<SVector3> &normals) {
+
+//  std::cout << "DBGTT: badEl = " << badEl->getNum() << ", _superEl = " << _superEl << std::endl;
+
+  switch (type) {
+    case TYPE_LIN:
+      createSuperElQuad(badEl, dist, baseVert, normals[0], normals[1]);
+      break;
+    case TYPE_TRI:
+      createSuperElPrism(badEl, dist, baseVert, normals[0], normals[1], normals[2]);
+      break;
+    case TYPE_QUA:
+      createSuperElHex(badEl, dist, baseVert, normals[0], normals[1], normals[2], normals[3]);
+      break;
+    default:
+      std::cout << "ERROR: SuperEl not implemented for element of type " << type << std::endl;
+      _superEl = 0;
+      break;
+  }
+
+  if (!_superEl) std::cout << "ERROR: SuperEl not created" << std::endl;
+
+}
+
+
+
+void SuperEl::createSuperElQuad(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+                                  const SVector3 &n0, const SVector3 &n1) {
+
+
+  if (baseVert.empty()) {
+    std::cout << "ERROR: Base edge not given for SuperEl linked to bad. el. " << badEl->getNum() << "\n";
+    _superEl = 0;
+    return;
+  }
+
+  // First-order vertices for super-element
+  MVertex *v0 = new MVertex(*baseVert[0]);
+  MVertex *v1 = new MVertex(*baseVert[1]);
+//  double dist = distFactor*distMaxStraight(badEl);
+  double v2x = v1->x()+n1.x()*dist, v2y = v1->y()+n1.y()*dist;
+  MVertex *v2 = new MVertex(v2x,v2y,0.);
+  double v3x = v0->x()+n0.x()*dist, v3y = v0->y()+n0.y()*dist;
+  MVertex *v3 = new MVertex(v3x,v3y,0.);
+//  std::cout << "DBGTT: createSuperElQuad: v0 = (" << v0->x() << "," << v0->y()
+//              << "), v1 = (" << v1->x() << "," << v1->y()
+//              << "), v2 = (" << v2->x() << "," << v2->y()
+//              << "), v3 = (" << v3->x() << "," << v3->y() << ")" << std::endl;
+
+  // Get basis functions
+  MQuadrangle quad1(v0,v1,v2,v3);
+  const int p = badEl->getPolynomialOrder();
+  const nodalBasis *quadBasis = quad1.getFunctionSpace(p);
+  const int Nv = quadBasis->getNumShapeFunctions();
+
+  // Store/create all vertices for super-element
+  _superVert.resize(Nv);
+  _superVert[0] = v0;                                                     // First-order vertices
+  _superVert[1] = v1;
+  _superVert[2] = v2;
+  _superVert[3] = v3;
+  for (int i=2; i<p+1; ++i) _superVert[2+i] = new MVertex(*baseVert[i]);  // HO vertices on base edge
+  for (int i=3+p; i<Nv; ++i) {                                            // Additional HO vertices
+    SPoint3 p;
+    quad1.pnt(quadBasis->points(i,0),quadBasis->points(i,1),0.,p);
+    _superVert[i] = new MVertex(p.x(),p.y(),0.);
+//    std::cout << "DBGTT: createSuperElQuad: add vertex (" << _superVert[i]->x()
+//                << "," << _superVert[i]->y()<< "," << _superVert[i]->z() << ")" << std::endl;
+  }
+
+  _superEl = new MQuadrangleN(_superVert,p);
+  _superEl0 = new MQuadrangle(_superVert[0],_superVert[1],_superVert[2],_superVert[3]);
+//  std::cout << "Created SuperEl at address " << _superEl << std::endl;
+
+}
+
+
+
+void SuperEl::createSuperElPrism(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+                                  const SVector3 &n0, const SVector3 &n1, const SVector3 &n2) {
+
+
+  if (baseVert.empty()) {
+    std::cout << "ERROR: Base edge not given for SuperEl linked to bad. el. " << badEl->getNum() << "\n";
+    _superEl = 0;
+    _superEl0 = 0;
+    return;
+  }
+
+  // First-order vertices for super-element
+  MVertex *v0 = new MVertex(*baseVert[0]);
+  MVertex *v1 = new MVertex(*baseVert[1]);
+  MVertex *v2 = new MVertex(*baseVert[2]);
+  double v3x = v0->x()+n0.x()*dist, v3y = v0->y()+n0.y()*dist,  v3z = v0->z()+n0.z()*dist;
+  MVertex *v3 = new MVertex(v3x,v3y,v3z);
+  double v4x = v1->x()+n1.x()*dist, v4y = v1->y()+n1.y()*dist,  v4z = v1->z()+n1.z()*dist;
+  MVertex *v4 = new MVertex(v4x,v4y,v4z);
+  double v5x = v2->x()+n2.x()*dist, v5y = v2->y()+n2.y()*dist,  v5z = v2->z()+n2.z()*dist;
+  MVertex *v5 = new MVertex(v5x,v5y,v5z);
+
+  // Get basis functions
+  MPrism prism1(v0,v1,v2,v3,v4,v5);
+  const int p = badEl->getPolynomialOrder();
+  const nodalBasis *prismBasis = prism1.getFunctionSpace(p);
+  const int Nv = prismBasis->getNumShapeFunctions();                        // Number of vertices in HO prism
+
+  // Store/create all vertices for super-element
+  if (badEl->getTypeForMSH() == MSH_PRI_18) {
+    _superVert.resize(18);
+    _superVert[0] = v0;                                                     // First-order vertices
+    _superVert[1] = v1;
+    _superVert[2] = v2;
+    _superVert[3] = v3;
+    _superVert[4] = v4;
+    _superVert[5] = v5;
+    _superVert[6] = new MVertex(*baseVert[3]);                              // HO vertices on base face
+    _superVert[9] = new MVertex(*baseVert[4]);
+    _superVert[7] = new MVertex(*baseVert[5]);
+    const int indAddVerts[9] = {8,10,11,12,13,14,15,16,17};                 // Additional HO vertices
+    SPoint3 p;
+    for (int i=0; i<9; ++i) {
+      const int &ind = indAddVerts[i];
+      prism1.pnt(prismBasis->points(ind,0),prismBasis->points(ind,1),prismBasis->points(ind,2),p);
+      _superVert[ind] = new MVertex(p.x(),p.y(),p.z());
+    }
+    _superEl = new MPrism18(_superVert);
+  }
+  else {
+    std::cout << "ERROR: Type of prism not supported for bad. el. " << badEl->getNum() << "\n";
+    _superEl = 0;
+    _superEl0 = 0;
+    return;
+  }
+
+  _superEl0 = new MPrism(_superVert[0],_superVert[1],_superVert[2],
+                         _superVert[3],_superVert[4],_superVert[5]);
+
+}
+
+
+
+void SuperEl::createSuperElHex(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+                               const SVector3 &n0, const SVector3 &n1, const SVector3 &n2, const SVector3 &n3) {
+
+
+  if (baseVert.empty()) {
+    std::cout << "ERROR: Base edge not given for SuperEl linked to bad. el. " << badEl->getNum() << "\n";
+    _superEl = 0;
+    _superEl0 = 0;
+    return;
+  }
+
+  // First-order vertices for super-element
+  MVertex *v0 = new MVertex(*baseVert[0]);
+  MVertex *v1 = new MVertex(*baseVert[1]);
+  MVertex *v2 = new MVertex(*baseVert[2]);
+  MVertex *v3 = new MVertex(*baseVert[3]);
+  double v4x = v0->x()+n0.x()*dist, v4y = v0->y()+n0.y()*dist,  v4z = v0->z()+n0.z()*dist;
+  MVertex *v4 = new MVertex(v4x,v4y,v4z);
+  double v5x = v1->x()+n1.x()*dist, v5y = v1->y()+n1.y()*dist,  v5z = v1->z()+n1.z()*dist;
+  MVertex *v5 = new MVertex(v5x,v5y,v5z);
+  double v6x = v2->x()+n2.x()*dist, v6y = v2->y()+n2.y()*dist,  v6z = v2->z()+n2.z()*dist;
+  MVertex *v6 = new MVertex(v6x,v6y,v6z);
+  double v7x = v3->x()+n3.x()*dist, v7y = v3->y()+n3.y()*dist,  v7z = v3->z()+n3.z()*dist;
+  MVertex *v7 = new MVertex(v7x,v7y,v7z);
+
+  // Get basis functions
+  MHexahedron *hex1 = new MHexahedron(v0,v1,v2,v3,v4,v5,v6,v7);
+  const int p = badEl->getPolynomialOrder();
+  const nodalBasis *prismBasis = hex1->getFunctionSpace(p);
+  const int Nv = prismBasis->getNumShapeFunctions();                         // Number of vertices in HO hex
+
+  // Store/create all vertices for super-element
+  if (badEl->getTypeForMSH() == MSH_HEX_27) {
+    _superVert.resize(27);
+    _superVert[0] = v0;                                                      // First-order vertices
+    _superVert[1] = v1;
+    _superVert[2] = v2;
+    _superVert[3] = v3;
+    _superVert[4] = v4;
+    _superVert[5] = v5;
+    _superVert[6] = v6;
+    _superVert[7] = v7;
+    _superVert[8] = new MVertex(*baseVert[4]);                                // HO vertices on base face
+    _superVert[11] = new MVertex(*baseVert[5]);                               // HO vertices on base face
+    _superVert[13] = new MVertex(*baseVert[6]);                               // HO vertices on base face
+    _superVert[9] = new MVertex(*baseVert[7]);
+    _superVert[20] = new MVertex(*baseVert[8]);
+    const int indAddVerts[14] = {10,12,14,15,16,17,18,19,21,22,23,24,25,26};  // Additional HO vertices
+    SPoint3 p;
+    for (int i=0; i<14; ++i) {
+      const int &ind = indAddVerts[i];
+      hex1->pnt(prismBasis->points(ind,0),prismBasis->points(ind,1),prismBasis->points(ind,2),p);
+      _superVert[ind] = new MVertex(p.x(),p.y(),p.z());
+    }
+    _superEl = new MHexahedron27(_superVert);
+  }
+  else {
+    std::cout << "ERROR: Type of prism not supported for bad. el. " << badEl->getNum() << "\n";
+    _superEl = 0;
+    _superEl0 = 0;
+    return;
+  }
+
+  _superEl0 = hex1;
+
+}
+
+
+
+bool SuperEl::isPointIn(const SPoint3 p) const {
+
+  double xyz[3] = {p.x(), p.y(), p.z()}, uvw[3];
+  _superEl0->xyz2uvw(xyz,uvw);
+  return _superEl0->isInside(uvw[0],uvw[1],uvw[2]);
+
+}
+
+
+
+bool SuperEl::straightToCurved(double *xyzS, double *xyzC) const {
+
+  double uvw[3];
+  _superEl0->xyz2uvw(xyzS,uvw);
+  if (!_superEl0->isInside(uvw[0],uvw[1],uvw[2])) return false;
+
+  SPoint3 pC;
+  _superEl->pnt(uvw[0],uvw[1],uvw[2],pC);
+  xyzC[0] = pC[0];
+  xyzC[1] = pC[1];
+  xyzC[2] = pC[2];
+
+  return true;
+
+}
+
+
+
+//std::string SuperEl::printPOS() {
+//
+//  std::vector<MVertex*> verts;
+//  _superEl0->getVertices(verts);
+//  std::string posStr(_superEl0->getStringForPOS());
+//  int n = _superEl0->getNumVertices();
+//
+//  std::ostringstream oss;
+//
+//  oss << posStr << "(";
+//  for(int i = 0; i < n; i++){
+//    if(i) oss << ",";
+//    oss << _superVert[i]->x() << "," <<  _superVert[i]->y() << "," <<  _superVert[i]->z();
+//  }
+//  oss << "){0";
+//  for(int i = 0; i < n; i++) oss << ",0.";
+//  oss << "};\n";
+//
+//  return oss.str();
+//
+//}
+
+
+
+std::string SuperEl::printPOS() {
+
+  std::vector<MVertex*> verts;
+  _superEl0->getVertices(verts);
+  std::string posStr(_superEl0->getStringForPOS());
+  int n = _superEl0->getNumVertices();
+
+  std::ostringstream oss;
+
+  oss << posStr << "(";
+  for(int i = 0; i < n; i++){
+    if(i) oss << ",";
+    oss << _superVert[i]->x() << "," <<  _superVert[i]->y() << "," <<  _superVert[i]->z();
+  }
+  oss << "){0";
+  for(int i = 0; i < n; i++) oss << ",0.";
+  oss << "};\n";
+
+  return oss.str();
+
+}

contrib/HighOrderMeshOptimizer/SuperEl.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@geuz.org>.
+//
+// Contributor: Thomas Toulorge
+
+#ifndef _SUPEREL_H_
+#define _SUPEREL_H_
+
+#include <string>
+#include "MElement.h"
+
+
+
+class SuperEl
+{
+public:
+
+  SuperEl(MElement *badEl, double dist, int type, const std::vector<MVertex*> &baseVert,
+          const std::vector<SVector3> &normals);
+  ~SuperEl() { _superVert.clear(); delete _superEl, _superEl0; }
+
+  bool isOK() const { return _superEl; }
+  bool isPointIn(const SPoint3 p) const;
+  bool straightToCurved(double *xyzS, double *xyzC) const;
+
+  std::string printPOS();
+
+  void printCoord() {
+    std::cout << "DBGTT: superEl -> ";
+    for(int i = 0; i < _superVert.size(); i++){
+      std::cout << "v" << i << " = (" << _superVert[i]->x() << "," <<  _superVert[i]->y() << "," <<  _superVert[i]->z() << ")";
+      if (i == _superVert.size()-1) std::cout << "\n"; else std::cout << ", ";
+    }
+
+  }
+
+private:
+
+  std::vector<MVertex*> _superVert;
+  MElement *_superEl, *_superEl0;
+
+  void createSuperElQuad(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+                         const SVector3 &n0, const SVector3 &n1);
+  void createSuperElPrism(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+                          const SVector3 &n0, const SVector3 &n1, const SVector3 &n2);
+  void createSuperElHex(MElement *badEl, double dist, const std::vector<MVertex*> &baseVert,
+                        const SVector3 &n0, const SVector3 &n1, const SVector3 &n2, const SVector3 &n3);
+
+
+};
+
+
+#endif