Merge branch 'develop' of https://github.com/SCOREC/core into develop

Author: yangf4

apf/apf.h

@@ -466,6 +466,9 @@ double computeCosAngleInTri(Mesh* m, MeshEntity* tri,
 double computeCosAngleInTet(Mesh* m, MeshEntity* tet,
     MeshEntity* e1, MeshEntity* e2, const Matrix3x3& Q);
 
+Vector3 computeFaceNormalAtEdgeInTet(Mesh* m, MeshEntity* tet,
+    MeshEntity* face, MeshEntity* edge, Matrix3x3 Q);
+
 Vector3 computeFaceNormalAtVertex(Mesh* m, MeshEntity* face,
     MeshEntity* vert, const Matrix3x3& Q);
 

apf/apfSimplexAngleCalcs.cc

@@ -116,22 +116,23 @@ static double computeEdgeEdgeCosAngleInTri(Mesh* m, MeshEntity* tri,
   return t1*t2;
 }
 
-static Vector3 computeFaceNormalAtEdge(Mesh* m, /*MeshEntity* tet,*/
+Vector3 computeFaceNormalAtEdgeInTet(Mesh* m, MeshEntity* tet,
     MeshEntity* face, MeshEntity* edge, Matrix3x3 Q)
 {
+  PCU_ALWAYS_ASSERT(m->getType(tet) == Mesh::TET);
   PCU_ALWAYS_ASSERT(m->getType(face) == Mesh::TRIANGLE);
   PCU_ALWAYS_ASSERT(m->getType(edge) == Mesh::EDGE);
 
-  MeshEntity* fe[3];
-  m->getDownward(face, 1, fe);
+  MeshEntity* te[6];
+  m->getDownward(tet, 1, te);
   int index = -1;
-  for (int i = 0; i < 3; i++)
-    if (fe[i] == edge) {
+  for (int i = 0; i < 6; i++)
+    if (te[i] == edge) {
       index = i;
       break;
     }
-  PCU_ALWAYS_ASSERT(index > -1 && index < 3);
-  Vector3 xi; // corresponding parametric coord in tri of the mid-edge point
+  PCU_ALWAYS_ASSERT(index > -1 && index < 6);
+  Vector3 xi; // corresponding parametric coord in tet for the mid-edge point
   switch (index) {
     case 0:
       xi = Vector3(0.5, 0.0, 0.0);
@@ -142,17 +143,44 @@ static Vector3 computeFaceNormalAtEdge(Mesh* m, /*MeshEntity* tet,*/
     case 2:
       xi = Vector3(0.0, 0.5, 0.0);
       break;
+    case 3:
+      xi = Vector3(0.0, 0.0, 0.5);
+      break;
+    case 4:
+      xi = Vector3(0.5, 0.0, 0.5);
+      break;
+    case 5:
+      xi = Vector3(0.0, 0.5, 0.5);
+      break;
     default:
       break;
   }
-  MeshElement* faceElem = createMeshElement(m, face);
+  MeshElement* tetElem = createMeshElement(m, tet);
   Matrix3x3 J;
-  getJacobian(faceElem, xi, J);
-  destroyMeshElement(faceElem);
+  getJacobian(tetElem, xi, J);
+  destroyMeshElement(tetElem);
+
+  // get the two columns of Jacobian corresponding to face tangents
+  MeshEntity* tf[4];
+  m->getDownward(tet, 2, tf);
+  index = -1;
+  for (int i = 0; i < 4; i++)
+    if (tf[i] == face) {
+      index = i;
+      break;
+    }
+  PCU_ALWAYS_ASSERT(index > -1 && index < 4);
 
   // transform Jacobian into Metric space
   J = getJacobianInMetric(J, Q);
-  return cross(J[0], J[1]).normalize();
+  if (index == 0)
+    return cross(J[0], J[1]).normalize();
+  else if (index == 1)
+    return cross(J[2], J[0]).normalize();
+  else if (index == 2)
+    return cross(J[2]-J[0], J[1]-J[0]).normalize();
+  else
+    return cross(J[1], J[2]).normalize();
 }
 
 
@@ -180,8 +208,8 @@ static double computeFaceFaceCosAngleInTet(Mesh* m, MeshEntity* tet,
       }
 
   PCU_ALWAYS_ASSERT(sharedEdge);
-  Vector3 f1n = computeFaceNormalAtEdge(m, /*tet,*/ face1, sharedEdge, Q);
-  Vector3 f2n = computeFaceNormalAtEdge(m, /*tet,*/ face2, sharedEdge, Q);
+  Vector3 f1n = computeFaceNormalAtEdgeInTet(m, tet, face1, sharedEdge, Q);
+  Vector3 f2n = computeFaceNormalAtEdgeInTet(m, tet, face2, sharedEdge, Q);
   // turn f2n so you get the inside angle
   f2n = Vector3(0.0, 0.0, 0.0) - f2n;
   return f1n * f2n;

ma/maFaceSplit.cc

@@ -172,6 +172,8 @@ Entity* makeSplitVertOnFace(Adapt* a, Entity* face)
   Vector param(0,0,0);
   if (a->input->shouldTransferParametric)
     transferParametricOnTriSplit(m, face, xi, param);
+  if (a->input->shouldTransferToClosestPoint)
+    transferToClosestPointOnTriSplit(m, face, xi, param);
   Entity* vert = buildVertex(a, c, point, param);
   st->onVertex(me, xi, vert);
   sf->interpolate(me, xi, vert);

ma/maInput.cc

@@ -30,6 +30,7 @@ void setDefaultValues(Input* in)
   in->shouldCoarsen = true;
   in->shouldSnap = in->mesh->canSnap();
   in->shouldTransferParametric = in->mesh->canSnap();
+  in->shouldTransferToClosestPoint = false;
   in->shouldHandleMatching = in->mesh->hasMatching();
   in->shouldFixShape = true;
   in->shouldForceAdaptation = false;
@@ -46,6 +47,8 @@ void setDefaultValues(Input* in)
     in->goodQuality = 0.2;
     //this basically turns off short edge removal...
     in->maximumEdgeRatio = 100.0;
+    //2D mesh adapt performs better if forceAdapt is on
+    in->shouldForceAdaptation = true;
   }
   in->shouldCheckQualityForDoubleSplits = false;
   in->validQuality = 1e-10;
@@ -91,8 +94,13 @@ void validateInput(Input* in)
     &&( ! in->mesh->canSnap()))
     rejectInput("user requested parametric coordinate transfer "
                 "but the geometric model does not support it");
-  if (in->shouldSnap && ( ! in->shouldTransferParametric))
-    rejectInput("snapping requires parametric coordinate transfer");
+  if (in->shouldTransferToClosestPoint
+    &&( ! in->mesh->canSnap()))
+    rejectInput("user requested transfer to closest point on model"
+                "but the geometric model does not support it");
+  if (in->shouldSnap && ( ! (in->shouldTransferParametric ||
+			     in->shouldTransferToClosestPoint)))
+    rejectInput("snapping requires parametric coordinate transfer or transfer to closest point");
   if ((in->mesh->hasMatching())
     &&( ! in->shouldHandleMatching))
     rejectInput("the mesh has matching entities but matched support is off");

ma/maInput.h

@@ -51,6 +51,9 @@ class Input
 /** \brief whether to transfer parametric coordinates
   \details requires modeler support, see gmi_reparam */
     bool shouldTransferParametric;
+/** \brief whether to transfer to the parametric coords of the closest point
+  \details requires modeler support, see gmi_closest_point */
+    bool shouldTransferToClosestPoint;
 /** \brief whether to update matched entity info (limited support) */
     bool shouldHandleMatching;
 /** \brief whether to run shape correction (default true) */

ma/maLayerTemplates.cc

@@ -51,6 +51,8 @@ void splitQuad_4(Refine* r, Entity* q, Entity** v)
   Vector param(0,0,0); //prevents uninitialized values
   if (a->input->shouldTransferParametric)
     transferParametricOnQuadSplit(m, q, sv[0] ,sv[2], 0.5, param);
+  if (a->input->shouldTransferToClosestPoint)
+    transferToClosestPointOnQuadSplit(m, q, sv[0] ,sv[2], 0.5, param);
   Entity* cv = buildVertex(a, m->toModel(q), point, param);
   a->solutionTransfer->onVertex(me,xi,cv);
   a->sizeField->interpolate(me,xi,cv);

ma/maRefine.cc

@@ -142,6 +142,8 @@ Entity* makeSplitVert(Refine* r, Entity* edge)
   Vector param(0,0,0); //prevents uninitialized values
   if (a->input->shouldTransferParametric)
     transferParametricOnEdgeSplit(m,edge,0.5,param);
+  if (a->input->shouldTransferToClosestPoint)
+    transferToClosestPointOnEdgeSplit(m,edge,0.5,param);
   Entity* vert = buildVertex(a,c,point,param);
   st->onVertex(me,xi,vert);
   sf->interpolate(me,xi,vert);

ma/maSnap.cc

@@ -597,6 +597,79 @@ void transferParametricOnQuadSplit(
   ma::transferParametricBetween(m, g, v, y, p);
 }
 
+void getClosestPointParametricCoordinates(
+    apf::Mesh* m,
+    Model* g,
+    double t,
+    Vector const& a,
+    Vector const& b,
+    Vector& p)
+{
+  Vector testPt = a * (1 - t) + b * t;
+  Vector targetPt;
+  m->getClosestPoint(g, testPt, targetPt, p);
+  (void) targetPt;
+}
+
+void transferToClosestPointOnEdgeSplit(
+    Mesh* m,
+    Entity* e,
+    double t,
+    Vector& p)
+{
+  Model* g = m->toModel(e);
+  int modelDimension = m->getModelType(g);
+  if (m->getDimension()==3 && modelDimension==3) return;
+  Entity* ev[2];
+  m->getDownward(e,0,ev);
+  Vector a = getPosition(m, ev[0]);
+  Vector b = getPosition(m, ev[1]);
+  getClosestPointParametricCoordinates(m, g, t, a, b, p);
+}
+
+void transferToClosestPointOnTriSplit(
+    Mesh* m,
+    Entity* face,
+    const Vector& xi,
+    Vector& param)
+{
+  Model* g = m->toModel(face);
+  int modelDimension = m->getModelType(g);
+  if (m->getDimension() == 3 && modelDimension == 3) return;
+  Entity* tv[3];
+  m->getDownward(face, 0, tv);
+  Vector x[3];
+
+  for (int i = 0; i < 3; i++) {
+    x[i] = getPosition(m, tv[i]);
+  }
+
+  Vector testPt = x[0] * xi[0] + x[1] * xi[1] + x[2] * xi[2];
+  Vector targetPt;
+
+  m->getClosestPoint(g, testPt, targetPt, param);
+  (void) targetPt;
+}
+
+void transferToClosestPointOnQuadSplit(
+    Mesh* m,
+    Entity* quad,
+    Entity* v01,
+    Entity* v32,
+    double y,
+    Vector& p)
+{
+  Model* g = m->toModel(quad);
+  int modelDimension = m->getModelType(g);
+  if (modelDimension==m->getDimension()) return;
+  Vector a = getPosition(m, v01);
+  Vector b = getPosition(m, v32);
+  Vector testPt = a * (1 - y) + b * y;
+  Vector targetPt;
+  m->getClosestPoint(g, testPt, targetPt, p);
+  (void) targetPt;
+}
+
 static void getSnapPoint(Mesh* m, Entity* v, Vector& x)
 {
   m->getPoint(v,0,x);

ma/maSnap.h

@@ -46,6 +46,31 @@ void transferParametricOnQuadSplit(
     double y,
     Vector& p);
 
+void getClosestPointParametricCoordinates(
+    apf::Mesh* m,
+    Model* g,
+    double t,
+    Vector const& a,
+    Vector const& b,
+    Vector& p);
+void transferToClosestPointOnEdgeSplit(
+    Mesh* m,
+    Entity* e,
+    double t,
+    Vector& p);
+void transferToClosestPointOnTriSplit(
+    Mesh* m,
+    Entity* face,
+    const Vector& xi,
+    Vector& param);
+void transferToClosestPointOnQuadSplit(
+    Mesh* m,
+    Entity* quad,
+    Entity* v01,
+    Entity* v32,
+    double y,
+    Vector& p);
+
 }
 
 #endif

mds/apfMDS.cc

@@ -891,7 +891,10 @@ void deriveMdlFromManifold(Mesh2* mesh, bool* isModelVert,
           "MeshTag name \"_classifn_data\" is used internally in this method\n");
   apf::MeshTag* classifnTag = mesh->createIntTag("_classifn_data", 2);
   int tagData[2], newTagData[2];
-  long minAvbl = 0;
+  long minAvbl = 1;
+  // This is set by during apf::construct, and using anything else leads to an
+  // additional region
+  long DEFAULT_REGION_ID = 0;
 
   PCU_ALWAYS_ASSERT_VERBOSE(!mesh->findTag("_vert_id"),
           "MeshTag name \"_vert_id\" is used internally in this method\n");
@@ -984,7 +987,8 @@ void deriveMdlFromManifold(Mesh2* mesh, bool* isModelVert,
           m->getIntTag(ent, classifnTag, tagData);
           mds_update_model_for_entity(m->mesh, id, tagData[0], tagData[1]);
         } else {
-          mds_update_model_for_entity(m->mesh, id, m->getDimension(), minAvbl);
+          mds_update_model_for_entity(m->mesh, id, m->getDimension(), 
+              DEFAULT_REGION_ID);
         }
       }
     }
@@ -1003,7 +1007,10 @@ void derive2DMdlFromManifold(Mesh2* mesh, bool* isModelVert,
           "MeshTag name \"_classifn_data\" is used internally in this method\n");
   apf::MeshTag* classifnTag = mesh->createIntTag("_classifn_data", 2);
   int tagData[2], newTagData[2];
-  long minAvbl = 0;
+  long minAvbl = 1;
+  // This is set by during apf::construct, and using anything else leads to an
+  // additional region
+  long DEFAULT_REGION_ID = 0;
 
   PCU_ALWAYS_ASSERT_VERBOSE(!mesh->findTag("_vert_id"),
           "MeshTag name \"_vert_id\" is used internally in this method\n");
@@ -1069,7 +1076,8 @@ void derive2DMdlFromManifold(Mesh2* mesh, bool* isModelVert,
           m->getIntTag(ent, classifnTag, tagData);
           mds_update_model_for_entity(m->mesh, id, tagData[0], tagData[1]);
         } else {
-          mds_update_model_for_entity(m->mesh, id, m->getDimension(), minAvbl);
+          mds_update_model_for_entity(m->mesh, id, m->getDimension(),
+              DEFAULT_REGION_ID);
         }
       }
     }