Coarsen and Fix Shape Mesh Adpat Improvements

CMakeLists.txt

@@ -172,6 +172,22 @@ if(ENABLE_CGNS)
   add_definitions(-DHAVE_CGNS)
 endif()
 
+if(DEBUG_FPP)
+  add_definitions(-DDEBUG_FPP)
+endif()
+
+if(ENABLE_EIGEN)
+  include(FetchContent)
+  FetchContent_Declare(
+    eigen
+    GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
+    GIT_TAG        5.0
+  )
+  FetchContent_MakeAvailable(eigen)
+  include_directories(${eigen_SOURCE_DIR})
+  add_definitions(-DEIGEN_ENABLED)
+endif()
+
 configure_file(SCOREC_config.h.in SCOREC_config.h)
 install(FILES "${CMAKE_BINARY_DIR}/SCOREC_config.h" DESTINATION include)
 include_directories(PUBLIC "$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}>")

apf/apf.cc

@@ -357,6 +357,13 @@ void getIntPoint(MeshElement* e, int order, int point, Vector3& param)
   param = p->param;
 }
 
+Vector3 getIntPoint(MeshElement* e, int order, int point)
+{
+  IntegrationPoint const* p = 
+    getIntegration(e->getType())->getAccurate(order)->getPoint(point);
+  return p->param;
+}
+
 double getIntWeight(MeshElement* e, int order, int point)
 {
   IntegrationPoint const* p = 

apf/apf.h

@@ -383,6 +383,15 @@ int countIntPoints(MeshElement* e, int order);
   */
 void getIntPoint(MeshElement* e, int order, int point, Vector3& param);
 
+
+/** \brief Get an integration point in an element.
+  *
+  * \param order The polynomial order of accuracy.
+  * \param point The integration point number.
+  * \returns The resulting local coordinate of the integration point.
+  */
+Vector3 getIntPoint(MeshElement* e, int order, int point);
+
 /** \brief Get the weight of an integration point in an element.
   *
   * \details All integration point tables in APF are scaled

apf/apfCavityOp.cc

@@ -119,6 +119,68 @@ void CavityOp::applyToDimension(int d)
   sharing = 0;
 }
 
+void CavityOp::prepareListMigration(EntityList& elements, EntityList::iterator& iter, std::function<void()> migrationPossible)
+{
+  MeshTag* migTag = mesh->createIntTag("ma_migrate_list", 1);
+  std::vector<EntityList::iterator> iterators(elements.size());
+  for (iter = elements.begin(); iter != elements.end();) {
+    int size = (int)iterators.size();
+    mesh->setIntTag(*iter, migTag, &size);
+    iterators.push_back(iter++);
+  }
+  mesh->onDestroy = [&](MeshEntity* e) {
+    if (!mesh->hasTag(e, migTag)) return;
+    int i = 0;
+    mesh->getIntTag(e, migTag, &i);
+    iter = elements.erase(iterators[i]);
+    movedByDeletion = true;
+  };
+
+  migrationPossible();
+
+  mesh->onDestroy = 0;
+  for (iter = elements.begin(); iter != elements.end();)
+    mesh->removeTag(*iter++, migTag);
+  mesh->destroyTag(migTag);
+}
+
+void CavityOp::applyToList(EntityList& elements)
+{
+  auto iter = elements.begin();
+  prepareListMigration(elements, iter, [&]() {
+  do {
+    delete sharing;
+    sharing = apf::getSharing(mesh);
+
+    isRequesting = false;
+    for (iter = elements.begin(); iter != elements.end();) {
+      if (sharing->isOwned(*iter)) {
+        if (setEntity(*iter) == OK) {
+          movedByDeletion=false;
+          apply();
+          if (movedByDeletion) continue;
+        }
+      }
+      ++iter;
+    }
+    /* request any non-local cavities.
+      note: it is critical that this loop
+      be separated from the one above for
+      mesh-modifying operators, since an apply()
+      call could change other overlapping cavities */
+    isRequesting = true;
+    for (iter = elements.begin(); iter != elements.end();) {
+      if (sharing->isOwned(*iter))
+        setEntity(*iter);
+      ++iter;
+    }
+  } while (tryToPull());
+
+  delete sharing;
+  sharing = 0;
+  });
+}
+
 bool CavityOp::requestLocality(MeshEntity** entities, int count)
 {
   bool areLocal = true;

apf/apfCavityOp.h

@@ -14,6 +14,7 @@
 #include "apfMesh.h"
 #include <vector>
 #include <cstring>
+#include <list>
 
 namespace apf {
 
@@ -87,6 +88,9 @@ class CavityOp
     virtual void apply() = 0;
     /** \brief parallel collective operation over entities of one dimension */
     void applyToDimension(int d);
+    /** \brief parallel collective operation over entities in list */
+    typedef std::list<MeshEntity*> EntityList;
+    void applyToList(EntityList& elements);
     /** \brief within setEntity, require that entities be made local */
     bool requestLocality(MeshEntity** entities, int count);
     /** \brief call before deleting a mesh entity during the operation */
@@ -102,6 +106,7 @@ class CavityOp
     bool tryToPull();
     void applyLocallyWithModification(int d);
     void applyLocallyWithoutModification(int d);
+    void prepareListMigration(EntityList& elements, EntityList::iterator& iter, std::function<void()> migrationPossible);
     bool canModify;
     bool movedByDeletion;
     MeshIterator* iterator;

apf/apfElement.cc

@@ -26,6 +26,10 @@ void Element::init(Field* f, MeshEntity* e, VectorElement* p)
   getNodeData();
 }
 
+Element::Element() 
+{
+}
+
 Element::Element(Field* f, MeshEntity* e)
 {
   init(f,e,0);

apf/apfElement.h

@@ -21,6 +21,7 @@ class VectorElement;
 class Element
 {
   public:
+    Element();
     Element(Field* f, MeshEntity* e);
     Element(Field* f, VectorElement* p);
     virtual ~Element();
@@ -36,12 +37,12 @@ class Element
     FieldShape* getFieldShape() {return field->getShape();}
     void getComponents(Vector3 const& xi, double* c);
     void getElementNodeData(NewArray<double>& d);
-  protected:
     void init(Field* f, MeshEntity* e, VectorElement* p);
+  protected:
     void getNodeData();
     Field* field;
     Mesh* mesh;
-    MeshEntity* entity;
+    MeshEntity* entity=0;
     EntityShape* shape;
     VectorElement* parent;
     int nen;

apf/apfElementOf.h

@@ -18,6 +18,10 @@ template <class T, class S = T>
 class ElementOf : public Element
 {
   public:
+    ElementOf():
+      Element()
+    {
+    }
     ElementOf(FieldOf<S>* f, MeshEntity* e):
       Element(f,e)
     {

apf/apfIntegrate.cc

@@ -5,6 +5,7 @@
  * BSD license as described in the LICENSE file in the top-level directory.
  */
 
+#include "apfVectorElement.h"
 #include "apfIntegrate.h"
 #include "apfMesh.h"
 #include "apf.h"
@@ -651,9 +652,8 @@ void Integrator::process(Mesh* m, int d)
   while ((entity = m->iterate(elements)))
   {
     if ( ! m->isOwned(entity)) continue;
-    MeshElement* e = createMeshElement(m,entity);
-    this->process(e);
-    destroyMeshElement(e);
+    MeshElement e(m->getCoordinateField(),entity);
+    this->process(&e);
   }
   m->end(elements);
   this->parallelReduce(m->getPCU());
@@ -695,9 +695,8 @@ double measure(MeshElement* e)
 
 double measure(Mesh* m, MeshEntity* e)
 {
-  MeshElement* me = createMeshElement(m,e);
-  double v = measure(me);
-  destroyMeshElement(me);
+  MeshElement me(m->getCoordinateField(),e);
+  double v = measure(&me);
   return v;
 }
 

apf/apfMatrix.cc

@@ -9,6 +9,9 @@
 #include "apf2mth.h"
 #include <mthQR.h>
 #include <pcu_util.h>
+#ifdef EIGEN_ENABLED
+#include <Eigen/Dense>
+#endif
 
 namespace apf {
 
@@ -69,6 +72,22 @@ int eigen(Matrix3x3 const& A,
           Vector<3>* eigenVectors,
           double* eigenValues)
 {
+#ifdef EIGEN_ENABLED
+  Eigen::Matrix3d matrix;
+  for (int i = 0; i < 3; ++i)
+  for (int j = 0; j < 3; ++j)
+      matrix(i, j) = A[i][j];
+
+  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> solver(matrix);
+  Eigen::Vector3d values = solver.eigenvalues();
+  Eigen::Matrix3d vectors = solver.eigenvectors();
+
+  for (unsigned i = 0; i < 3; ++i)
+    eigenValues[i] = values(i);
+  for (unsigned i = 0; i < 3; ++i)
+  for (unsigned j = 0; j < 3; ++j)
+    eigenVectors[j][i] = vectors(i,j);
+#else
   mth::Matrix<double,3,3> A2 = to_mth(A);
   mth::Matrix<double,3,3> L;
   mth::Matrix<double,3,3> Q;
@@ -79,6 +98,7 @@ int eigen(Matrix3x3 const& A,
   for (unsigned i = 0; i < 3; ++i)
   for (unsigned j = 0; j < 3; ++j)
     eigenVectors[j][i] = Q(i,j);
+#endif
   return 3;
 }
 

apf/apfMesh.cc

@@ -174,6 +174,12 @@ Mesh::~Mesh()
   delete coordinateField;
 }
 
+void Mesh::clearModelCache()
+{
+  periodicRanges.clear();
+  boundingBoxes.clear();
+}
+
 int Mesh::getModelType(ModelEntity* e)
 {
   return gmi_dim(getModel(), (gmi_ent*)e);
@@ -228,6 +234,21 @@ bool Mesh::getPeriodicRange(ModelEntity* g, int axis, double range[2])
   return gmi_periodic(getModel(), e, axis);
 }
 
+bool Mesh::getPeriodicRangeCached(ModelEntity* g, int axis, double range[2])
+{
+  auto it = periodicRanges.find({g, axis});
+  if (it != periodicRanges.end()) {
+    range[0] = it->second.first[0];
+    range[1] = it->second.first[1];
+    return it->second.second;
+  }
+  else{
+    bool output = getPeriodicRange(g, axis, range);
+    periodicRanges[{g, axis}] = {{range[0], range[1]}, output};
+    return output;
+  }
+}
+
 void Mesh::getClosestPoint(ModelEntity* g, Vector3 const& from,
     Vector3& to, Vector3& p)
 {
@@ -284,6 +305,20 @@ bool Mesh::isInClosureOf(ModelEntity* g, ModelEntity* target){
   return (res == 1) ? true : false;
 }
 
+void Mesh::boundingBoxCached(ModelEntity* g,
+    Vector3& bmin, Vector3& bmax)
+{
+  auto it = boundingBoxes.find(g);
+  if (it != boundingBoxes.end()) {
+    bmin = it->second.first;
+    bmax = it->second.second;
+  }
+  else{
+    boundingBox(g, bmin, bmax);
+    boundingBoxes.insert(it, {g, {bmin, bmax}});
+  }
+}
+
 void Mesh::boundingBox(ModelEntity* g,
     Vector3& bmin, Vector3& bmax)
 {

apf/apfMesh.h

@@ -12,11 +12,14 @@
     \brief The APF Mesh interface*/
 
 #include <vector>
+#include <array>
 #include <map>
+#include <functional>
 #include <set>
 #include <PCU.h>
 #include "apfVector.h"
 #include "apfDynamicArray.h"
+#include "gmi.h"
 
 struct gmi_model;
 
@@ -288,6 +291,8 @@ class Mesh
     virtual gmi_model* getModel() = 0;
     /** \brief set the geometric model */
     virtual void setModel(gmi_model* newModel) = 0;
+    /** \brief clear cached model info */
+    void clearModelCache();
     /** \brief return the model entity dimension */
     int getModelType(ModelEntity* e);
     /** \brief get the dimension-unique model entity identifier */
@@ -309,6 +314,11 @@ class Mesh
       \returns true if (g) is periodic along this axis */
     bool getPeriodicRange(ModelEntity* g, int axis,
         double range[2]);
+    /** \brief get the periodic properties of a model entity and chache the value
+      \param range if periodic, the parametric range
+      \returns true if (g) is periodic along this axis */
+    bool getPeriodicRangeCached(ModelEntity* g, int axis, 
+        double range[2]);
     /** \brief get closest point on geometry */
     void getClosestPoint(ModelEntity* g, Vector3 const& from,
         Vector3& to, Vector3& p);
@@ -325,6 +335,8 @@ class Mesh
     bool isInClosureOf(ModelEntity* g, ModelEntity* target);
     /** \brief get the bounding box of the model entity g */
     void boundingBox(ModelEntity* g, Vector3& bmin, Vector3& bmax);
+    /** \brief get the bounding box of the model entity g and caches the value*/
+    void boundingBoxCached(ModelEntity* g, Vector3& bmin, Vector3& bmax);
     /** \brief checks if p is on model g */
     bool isOnModel(ModelEntity* g, Vector3 p, double scale);
     /** \brief get the distribution of the mesh's coordinate field */
@@ -401,12 +413,16 @@ class Mesh
     void switchPCU(pcu::PCU *newPCU);
     /** \brief true if any associated fields use array storage */
     bool hasFrozenFields;
+    std::function<void(MeshEntity*)> onDestroy=0;
   protected:
     Field* coordinateField;
     std::vector<Field*> fields;
     std::vector<Numbering*> numberings;
     std::vector<GlobalNumbering*> globalNumberings;
     pcu::PCU* pcu_;
+    std::map<std::pair<ModelEntity*, int>,
+             std::pair<std::array<double, 2>, bool>> periodicRanges;
+    std::map<ModelEntity*, std::pair<Vector3, Vector3>> boundingBoxes;
 };
 
 /** \brief run consistency checks on an apf::Mesh structure

apf/apfMesh2.h

@@ -102,6 +102,7 @@ class Mesh2 : public Mesh
   downward adjacencies */
     void destroy(MeshEntity* e)
     {
+      if (onDestroy) onDestroy(e);
       requireUnfrozen();
       destroy_(e);
     }