Coarsen and Snapping

.github/workflows/simmetrix_enabled_pr_comment_trigger_self_hosted.yml

@@ -37,7 +37,7 @@ jobs:
           module use /opt/scorec/spack/rhel9/v0201_4/lmod/linux-rhel9-x86_64/Core/
           module load gcc/12.3.0-iil3lno
           module load mpich/4.1.1-xpoyz4t
-          module load simmetrix-simmodsuite/2024.0-240119dev-7abimo4
+          module load simmetrix-simmodsuite/2025.0-250108dev-llxq6sk
           module load zoltan/3.83-hap4ggo
           module load cmake/3.26.3-2duxfcd
           set -e
@@ -70,7 +70,7 @@ jobs:
           module use /opt/scorec/spack/rhel9/v0201_4/lmod/linux-rhel9-x86_64/Core/
           module load gcc/12.3.0-iil3lno
           module load mpich/4.1.1-xpoyz4t
-          module load simmetrix-simmodsuite/2024.0-240119dev-7abimo4
+          module load simmetrix-simmodsuite/2025.0-250108dev-llxq6sk
           module load zoltan/3.83-hap4ggo
           module load cmake/3.26.3-2duxfcd
           module load cgns/develop-cc4dfwp

CMakeLists.txt

@@ -6,7 +6,7 @@ endif()
 # This is the top level CMake file for the SCOREC build
 cmake_minimum_required(VERSION 3.12)
 
-project(SCOREC VERSION 4.0.0 LANGUAGES CXX C)
+project(SCOREC VERSION 4.1.0 LANGUAGES CXX C)
 
 include(cmake/bob.cmake)
 include(cmake/xsdk.cmake)

gmi_cap/gmi_cap.cc

@@ -541,3 +541,7 @@ GDBI* gmi_export_cap(gmi_model* m)
   cap_model* cm = (cap_model*)m;
   return cm->geomInterface;
 }
+
+int gmi_cap_test(struct gmi_model* model) {
+  return model->ops == &ops ? 1 : 0;
+}

gmi_cap/gmi_cap.h

@@ -57,6 +57,7 @@ struct gmi_ent;
  * \note This call is required before calling any other gmi_cap functions.
  */
 void gmi_cap_start(void);
+
 /**
  * \brief Finalize gmi_cap library.
  *
@@ -94,6 +95,15 @@ struct gmi_model* gmi_cap_load(const char* creFileName);
  */
 void gmi_cap_write(struct gmi_model* model, const char* creFileName);
 
+/**
+ * \brief Test if model is a Capstone gmi_model.
+ *
+ * \param model An abstract gmi_model which may have an underlying Capstone
+ * database.
+ * \return 1 if model is a Capstone gmi_model and 0 otherwise.
+ */
+int gmi_cap_test(struct gmi_model* model);
+
 #ifdef __cplusplus
 
 /**

ma/CMakeLists.txt

@@ -85,6 +85,10 @@ target_link_libraries(ma
      pcu
    )
 
+if (ENABLE_CAPSTONE)
+target_link_libraries(ma PRIVATE apf_cap)
+endif()
+
 scorec_export_library(ma)
 
 bob_end_subdir()

ma/ma.cc

@@ -100,7 +100,7 @@ void adaptVerbose(Input* in, bool verbose)
   int count = 0;
   double lMax = ma::getMaximumEdgeLength(a->mesh, a->sizeField);
   print(a->mesh->getPCU(), "Maximum (metric) edge length in the mesh is %f", lMax);
-  while (lMax > 1.5) {
+  while (lMax > MAXLENGTH) {
     print(a->mesh->getPCU(), "%dth additional refine-snap call", count);
     refine(a);
     snap(a);

ma/maCoarsen.cc

@@ -7,15 +7,81 @@
   of the SCOREC Non-Commercial License this program is distributed under.
 
 *******************************************************************************/
+/*
+This file contains two coarsening alogrithms. 
+  1. coarsenMultiple: will repeatedly create independent sets and collapse the mesh 
+      until there are no more short edges left. This code is currently only used for testing 
+      since more work needs to be done to achieve better performance and have the code work in 
+      parallel, at which point the other coarsen algorithm can be deleted.
+  2. coarsen: The first will create one independent set and then collapse all of the edges 
+      in that independent set. 
+*/
 #include "maCoarsen.h"
 #include "maAdapt.h"
 #include "maCollapse.h"
 #include "maMatchedCollapse.h"
 #include "maOperator.h"
+#include "maDBG.h"
 #include <pcu_util.h>
+#include "apfShape.h"
+#include <vector>
+#include <list>
+#include <algorithm>
 
 namespace ma {
 
+namespace {
+
+//Measures edge length and stores the result so it doesn't have to be calculated again
+double getLength(Adapt* a, Tag* lengthTag, Entity* edge)
+{
+  double length = 0;
+  if (a->mesh->hasTag(edge, lengthTag))
+    a->mesh->getDoubleTag(edge, lengthTag, &length);
+  else {
+    length = a->sizeField->measure(edge);
+    a->mesh->setDoubleTag(edge, lengthTag, &length);
+  }
+  return length;
+}
+
+//Make sure that a collapse will not create an edge longer than the max
+bool collapseSizeCheck(Adapt* a, Entity* vertex, Entity* edge, apf::Up& adjacent)
+{
+  Entity* vCollapse = getEdgeVertOppositeVert(a->mesh, edge, vertex);
+  for (int i=0; i<adjacent.n; i++) {
+    Entity* newEdgeVerts[2]{vertex, getEdgeVertOppositeVert(a->mesh, adjacent.e[i], vCollapse)};
+    Entity* newEdge = a->mesh->createEntity(apf::Mesh::EDGE, 0, newEdgeVerts);
+    double length = a->sizeField->measure(newEdge);
+    destroyElement(a, newEdge);
+    if (length > MAXLENGTH) return false;
+  }
+  return true;
+}
+
+bool tryCollapseEdge(Adapt* a, Entity* edge, Entity* keep, Collapse& collapse, apf::Up& adjacent)
+{
+  PCU_ALWAYS_ASSERT(a->mesh->getType(edge) == apf::Mesh::EDGE);
+  bool alreadyFlagged = true;
+  if (keep) alreadyFlagged = getFlag(a, keep, DONT_COLLAPSE);
+  if (!alreadyFlagged) setFlag(a, keep, DONT_COLLAPSE);
+
+  double quality = a->input->shouldForceAdaptation ? a->input->validQuality 
+                                                  : a->input->goodQuality;
+
+  bool result = false;
+  if (collapse.setEdge(edge) && 
+      collapse.checkClass() &&
+      collapse.checkTopo() &&
+      collapseSizeCheck(a, keep, edge, adjacent) &&
+      collapse.tryBothDirections(quality)) {
+    result = true;
+  }  
+  if (!alreadyFlagged) clearFlag(a, keep, DONT_COLLAPSE);
+  return result;
+}
+}
+
 class CollapseChecker : public apf::CavityOp
 {
   public:
@@ -253,4 +319,153 @@ bool coarsen(Adapt* a)
   return true;
 }
 
+/*
+  To be used after collapsing a vertex to flag adjacent vertices as NEED_NOT_COLLAPSE. This will create a
+  set of collapses where no two adjacent are vertices collapsed. Will also clear adjacent vertices for 
+  collapse in next independent set since they might succeed after a collapse.
+*/
+void flagIndependentSet(Adapt* a, apf::Up& adjacent, size_t& checked)
+{
+  for (int adj=0; adj < adjacent.n; adj++) {
+    Entity* vertices[2];
+    a->mesh->getDownward(adjacent.e[adj],0, vertices);
+    for (int v = 0; v < 2; v++) {
+      setFlag(a, vertices[v], NEED_NOT_COLLAPSE);
+      if (getFlag(a, vertices[v], CHECKED)){
+        clearFlag(a, vertices[v], CHECKED); //needs to be checked again in next independent set
+        checked--;
+      }
+    }
+  }
+}
+
+struct EdgeLength
+{
+  Entity* edge;
+  double length;
+  bool operator<(const EdgeLength& other) const {
+    return length < other.length;
+  }
+};
+
+/*
+  Given an iterator pointing to a vertex we will collapse the shortest adjacent edge and try the next
+  shorted until one succeeds and then it will expand independent set. In Li's thesis it only attempts
+  to collapse the shortest edge, but this gave us better results.
+*/
+bool collapseShortest(Adapt* a, Collapse& collapse, std::list<Entity*>& shortEdgeVerts, std::list<Entity*>::iterator& itr, size_t& checked, apf::Up& adjacent, Tag* lengthTag)
+{
+  Entity* vertex = *itr;
+  std::vector<EdgeLength> sorted;
+  for (int i=0; i < adjacent.n; i++) {
+    double length = getLength(a, lengthTag, adjacent.e[i]);
+    EdgeLength measured{adjacent.e[i], length};
+    if (measured.length > MINLENGTH) continue;
+    sorted.push_back(measured);
+  }
+  if (sorted.size() == 0) { //performance optimization, will rarely result in a missed edge
+    itr = shortEdgeVerts.erase(itr);
+    return false;
+  }
+  std::sort(sorted.begin(), sorted.end());
+  for (size_t i=0; i < sorted.size(); i++) {
+    Entity* keepVertex = getEdgeVertOppositeVert(a->mesh, sorted[i].edge, vertex);
+    if (!tryCollapseEdge(a, sorted[i].edge, keepVertex, collapse, adjacent)) continue;
+    flagIndependentSet(a, adjacent, checked);
+    itr = shortEdgeVerts.erase(itr);
+    collapse.destroyOldElements();
+    return true;
+  }
+  setFlag(a, vertex, CHECKED);
+  checked++;
+  return false;
+}
+
+/*
+  Iterates through shortEdgeVerts until it finds a vertex that is adjacent to an 
+  independent set. We want our collapses to touch the independent set in order to
+  reduce adjacent collapses, since collapsing adjacent vertices will result in a
+  lower quality mesh.
+*/
+bool getAdjIndependentSet(Adapt* a, std::list<Entity*>& shortEdgeVerts, std::list<Entity*>::iterator& itr, bool& independentSetStarted, apf::Up& adjacent)
+{
+  size_t numItr=0;
+  do {
+    numItr++;
+    if (itr == shortEdgeVerts.end()) itr = shortEdgeVerts.begin();
+    Entity* vertex = *itr;
+    if (getFlag(a, vertex, CHECKED)) {itr++; continue;} //Already tried to collapse
+    a->mesh->getUp(vertex, adjacent);
+    if (!independentSetStarted) return true;
+    if (getFlag(a, vertex, NEED_NOT_COLLAPSE)) {itr++; continue;} //Too close to last collapse
+    for (int i=0; i < adjacent.n; i++)
+    {
+      Entity* opposite = getEdgeVertOppositeVert(a->mesh, adjacent.e[i], vertex);
+      if (getFlag(a, opposite, NEED_NOT_COLLAPSE)) return true; //Touching independent set
+    }
+    itr++;
+  } while (numItr < shortEdgeVerts.size());
+  for (auto v : shortEdgeVerts) clearFlag(a, v, NEED_NOT_COLLAPSE);
+  independentSetStarted = false;
+  return false;
+}
+
+//returns a list of vertices that bound a short edge and flags them
+std::list<Entity*> getShortEdgeVerts(Adapt* a, Tag* lengthTag)
+{
+  std::list<Entity*> shortEdgeVerts;
+  Iterator* it = a->mesh->begin(1);
+  Entity* edge;
+  while ((edge = a->mesh->iterate(it))) 
+  {
+    double length = getLength(a, lengthTag, edge);
+    if (length > MINLENGTH) continue;
+    Entity* vertices[2];
+    a->mesh->getDownward(edge,0,vertices);
+    for (int i = 0; i < 2; i++) {
+      if (getFlag(a, vertices[i], CHECKED)) continue;
+      setFlag(a, vertices[i], CHECKED);
+      shortEdgeVerts.push_back(vertices[i]);
+    }
+  }
+  for (auto v : shortEdgeVerts) clearFlag(a, v, CHECKED);
+  a->mesh->end(it);
+  return shortEdgeVerts;
+}
+
+/*
+  Follows the alogritm in Li's thesis in order to coarsen all short edges 
+  in a mesh while maintaining a decent quality mesh.
+*/
+bool coarsenMultiple(Adapt* a)
+{
+  if (!a->input->shouldCoarsen)
+    return false;
+  double t0 = pcu::Time();
+  Tag* lengthTag = a->mesh->createDoubleTag("edge_length", 1);
+  std::list<Entity*> shortEdgeVerts = getShortEdgeVerts(a, lengthTag);
+
+  Collapse collapse;
+  collapse.Init(a);
+  int success = 0;
+  size_t checked = 0;
+  bool independentSetStarted = false;
+  std::list<Entity*>::iterator itr = shortEdgeVerts.begin();
+  while (checked < shortEdgeVerts.size())
+  {
+    apf::Up adjacent;
+    if (!getAdjIndependentSet(a, shortEdgeVerts, itr, independentSetStarted, adjacent)) continue;
+    if (collapseShortest(a, collapse, shortEdgeVerts, itr, checked, adjacent, lengthTag)) {
+      independentSetStarted=true;
+      success++;
+    }
+  }
+  ma::clearFlagFromDimension(a, NEED_NOT_COLLAPSE | CHECKED, 0);
+  a->mesh->destroyTag(lengthTag);
+  double t1 = pcu::Time();
+  print(a->mesh->getPCU(), "coarsened %d edges in %f seconds", success, t1-t0);
+  return true;
+}
+
+
 }

ma/maCoarsen.h

@@ -14,6 +14,7 @@ namespace ma {
 
 class Adapt;
 
+bool coarsenMultiple(Adapt* a);
 bool coarsen(Adapt* a);
 bool coarsenLayer(Adapt* a);