Quadratic simx meshes can be imported to Bezier

The inversion preserves the validity of the curved elements

Author: mortezah

crv/crvBezier.h

@@ -131,7 +131,7 @@ void getGregoryBlendedTransformationCoefficients(int blend, int type,
     apf::NewArray<double>& c);
 
 /** \brief a per entity version of above */
-void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e);
+void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e, bool isNew = false);
 
 /** \brief compute the matrix to transform between Bezier and Lagrange Points
     \details this is a support function, not actual ever needed.

crv/crvCurveMesh.cc

@@ -44,7 +44,7 @@ void convertInterpolationPoints(apf::Mesh2* m, apf::MeshEntity* e,
   apf::destroyElement(elem);
 }
 
-void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e)
+void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e, bool isNew)
 {
   PCU_ALWAYS_ASSERT(m->canSnap());
   int type = m->getType(e);
@@ -56,9 +56,18 @@ void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e)
     m->setPoint(e,0,pt);
     return;
   }
+  // e is an edge or a face
+  // either way, get a length-scale by computing
+  // the distance b/w first two downward verts
+  apf::MeshEntity* down[12];
+  m->getDownward(e, 0, down);
+  apf::Vector3 p0, p1;
+  m->getPoint(down[0], 0, p0);
+  m->getPoint(down[1], 0, p1);
+  double lengthScale = (p1 - p0).getLength();
   apf::FieldShape * fs = m->getShape();
   int non = fs->countNodesOn(type);
-  apf::Vector3 p, xi, pt(0,0,0);
+  apf::Vector3 p, xi, pt0, pt(0,0,0);
   for(int i = 0; i < non; ++i){
     apf::ModelEntity* g = m->toModel(e);
     fs->getNodeXi(type,i,xi);
@@ -67,7 +76,13 @@ void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e)
     else
       transferParametricOnTriSplit(m,e,xi,p);
     m->snapToModel(g,p,pt);
-    m->setPoint(e,i,pt);
+    if (isNew || !m->canGetClosestPoint()) {
+      m->setPoint(e,i,pt);
+      continue;
+    }
+    m->getPoint(e,i,pt0);
+    if (!m->isOnModel(g, pt0, lengthScale))
+      m->setPoint(e,i,pt);
   }
 }
 

crv/crvShapeHandler.cc

@@ -149,7 +149,7 @@ class BezierTransfer : public ma::SolutionTransfer
       for (int i = 0; i < ne; ++i){
         if(shouldSnap && midEdgeVerts[i] &&
             isBoundaryEntity(mesh,midEdgeVerts[i])){
-          snapToInterpolate(mesh,midEdgeVerts[i]);
+          snapToInterpolate(mesh,midEdgeVerts[i],true);
         }
       }
       int np = mesh->getShape()->getEntityShape(parentType)->countNodes();
@@ -173,7 +173,7 @@ class BezierTransfer : public ma::SolutionTransfer
 
           // do snapping here, inside refinement
           if (isBdryEnt && shouldSnap){
-            snapToInterpolate(mesh,newEntities[i]);
+            snapToInterpolate(mesh,newEntities[i],true);
           } else if (useLinear && !isBdryEnt) {
             // boundary entities that don't get snapped should interpolate
             if(childType == apf::Mesh::EDGE){
@@ -538,7 +538,7 @@ class BezierHandler : public ma::ShapeHandler
         if (newType != apf::Mesh::EDGE)
         {
           if (snap && P > 2){
-            snapToInterpolate(mesh,newEntities[i]);
+            snapToInterpolate(mesh,newEntities[i],true);
             convertInterpolationPoints(mesh,newEntities[i],n,ni,bt->coeffs[2]);
           } else {
             for (int j = 0; j < ni; ++j){
@@ -548,7 +548,7 @@ class BezierHandler : public ma::ShapeHandler
           }
         } else {
           if (snap){
-            snapToInterpolate(mesh,newEntities[i]);
+            snapToInterpolate(mesh,newEntities[i],true);
             convertInterpolationPoints(mesh,newEntities[i],n,ni,bt->coeffs[1]);
           } else {
             setLinearEdgePoints(mesh,newEntities[i]);