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,40 @@ void CavityOp::applyToDimension(int d)
   sharing = 0;
 }
 
+void CavityOp::applyToList(std::list<MeshEntity*>& elements)
+{
+  do {
+    delete sharing;
+    sharing = apf::getSharing(mesh);
+
+    isRequesting = false;
+    for (auto iter = elements.begin(); iter != elements.end();) {
+      if (sharing->isOwned(*iter)) {
+        if (setEntity(*iter) == OK) {
+          movedByDeletion=false;
+          apply();
+          if (movedByDeletion) {iter = elements.erase(iter); 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 (auto 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,8 @@
 #include "apfMesh.h"
 #include <vector>
 #include <cstring>
+#include <functional>
+#include <list>
 
 namespace apf {
 
@@ -87,12 +89,15 @@ 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 */
+    void applyToList(std::list<MeshEntity*>& 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 */
     void preDeletion(MeshEntity* e);
     /** \brief mesh pointer for convenience */
     Mesh* mesh;
+    bool movedByDeletion;
   private:
     typedef std::vector<MeshEntity*> Requests;
     Requests requests;
@@ -102,8 +107,8 @@ class CavityOp
     bool tryToPull();
     void applyLocallyWithModification(int d);
     void applyLocallyWithoutModification(int d);
+
     bool canModify;
-    bool movedByDeletion;
     MeshIterator* iterator;
   protected:
     Sharing* sharing;

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,8 +37,8 @@ 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;

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/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

@@ -15,6 +15,7 @@
 #include <pcu_util.h>
 #include <lionPrint.h>
 #include <algorithm>
+#include <array>
 
 namespace apf {
 
@@ -221,11 +222,23 @@ void Mesh::getParamOn(ModelEntity* g, MeshEntity* e, Vector3& p)
   gmi_reparam(getModel(), from, &from_p[0], to, &p[0]);
 }
 
+std::map<std::pair<gmi_ent*, int>,
+  std::pair<std::array<double, 2>, bool>> periodicRanges;
 bool Mesh::getPeriodicRange(ModelEntity* g, int axis, double range[2])
 {
   gmi_ent* e = (gmi_ent*)g;
-  gmi_range(getModel(), e, axis, range);
-  return gmi_periodic(getModel(), e, axis);
+  auto it = periodicRanges.find({e, axis});
+  if (it != periodicRanges.end()) {
+    range[0] = it->second.first[0];
+    range[1] = it->second.first[1];
+    return it->second.second;
+  }
+  else{
+    gmi_range(getModel(), e, axis, range);
+    bool output = gmi_periodic(getModel(), e, axis);
+    periodicRanges[{e, axis}] = {{range[0], range[1]}, output};
+    return output;
+  }
 }
 
 void Mesh::getClosestPoint(ModelEntity* g, Vector3 const& from,

apf/apfVectorElement.cc

@@ -10,6 +10,11 @@
 
 namespace apf {
 
+VectorElement::VectorElement(): 
+  ElementOf<Vector3>()
+{
+}
+
 VectorElement::VectorElement(VectorField* f, MeshEntity* e):
   ElementOf<Vector3>(f,e)
 {

apf/apfVectorElement.h

@@ -18,6 +18,7 @@ class VectorField;
 class VectorElement : public ElementOf<Vector3>
 {
   public:
+    VectorElement();
     VectorElement(VectorField* f, MeshEntity* e);
     VectorElement(VectorField* f, VectorElement* p);
     virtual ~VectorElement() {}

ma/CMakeLists.txt

@@ -28,6 +28,8 @@ set(SOURCES
   maSolutionTransferHelper.cc
   maSnap.cc
   maEdgeSwap.cc
+  maFaceSwap.cc
+  maFixShape.cc
   maShape.cc
   maShapeHandler.cc
   maQuality.cc
@@ -59,6 +61,7 @@ set(HEADERS
   maMesh.h
   maSize.h
   maShape.h
+  maFixShape.h
   maTables.h
   maSolutionTransfer.h
   maExtrude.h