Merging develop into master

Author: Nightly Bot

cdash/nightly.sh

@@ -11,7 +11,7 @@ unset MODULEPATH
 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

cmake/FindSimModSuite.cmake

@@ -84,7 +84,7 @@ string(REGEX REPLACE
   "${SIM_VERSION}")
 
 set(MIN_VALID_SIM_VERSION 15.0.191017)
-set(MAX_VALID_SIM_VERSION 2025.0.241016)
+set(MAX_VALID_SIM_VERSION 2025.0.250108)
 if( ${SKIP_SIMMETRIX_VERSION_CHECK} )
   message(STATUS "Skipping Simmetrix SimModSuite version check."
     " This may result in undefined behavior")

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);
 

ma/maCollapse.cc

@@ -55,6 +55,34 @@ bool Collapse::tryThisDirectionNoCancel(double qualityToBeat)
   return true;
 }
 
+/*Sometimes the snapping procedure will attempt to collapse a edge that was just 
+refined. This function will prevent that when the growth of the edge is above a 
+thresh hold where it might not reach the target edge length*/
+bool Collapse::edgesGoodSize() {
+  PCU_ALWAYS_ASSERT(elementsToKeep.size());
+  PCU_ALWAYS_ASSERT(elementsToKeep.size() == newElements.size());
+  size_t i=0;
+  double maxSize=0;
+  double ratioAtMaxSize=0;
+  APF_ITERATE(EntitySet,elementsToKeep,it) {
+    Entity* edgesBefore[6];
+    adapt->mesh->getDownward(*it,1,edgesBefore);
+    Entity* edgesAfter[6];
+    adapt->mesh->getDownward(newElements[i++],1,edgesAfter);
+    for (int j=0; j<6; j++) {
+      double sizeBefore = adapt->sizeField->measure(edgesBefore[j]);
+      double sizeAfter = adapt->sizeField->measure(edgesAfter[j]);
+      double ratio = sizeAfter/sizeBefore;
+      if (sizeAfter > maxSize) {
+        maxSize = sizeAfter;
+        ratioAtMaxSize = ratio;
+      }
+    }
+  }
+  if (maxSize > MAXLENGTH && ratioAtMaxSize > MAXLENGTHRATIO) return false;
+  return true;
+}
+
 bool Collapse::tryThisDirection(double qualityToBeat)
 {
   if (!tryThisDirectionNoCancel(qualityToBeat)) {
@@ -74,17 +102,20 @@ bool Collapse::tryBothDirections(double qualityToBeat)
     qualityToBeat = std::min(adapt->input->goodQuality,
         std::max(getOldQuality(),adapt->input->validQuality));
 
-  if (tryThisDirection(qualityToBeat))
+  if (tryThisDirectionNoCancel(qualityToBeat))
     return true;
-  if ( ! getFlag(adapt,vertToKeep,COLLAPSE))
-    return false;
-  std::swap(vertToKeep,vertToCollapse);
-  computeElementSets();
-  if (!adapt->input->shouldForceAdaptation)
-    qualityToBeat = std::min(adapt->input->goodQuality,
-        std::max(getOldQuality(),adapt->input->validQuality));
+  else destroyNewElements();
+
+  if (getFlag(adapt,vertToKeep,COLLAPSE)) {
+    std::swap(vertToKeep,vertToCollapse);
+    computeElementSets();
+    if (tryThisDirectionNoCancel(qualityToBeat))
+      return true;
+    else destroyNewElements();
+  }
 
-  return tryThisDirection(qualityToBeat);
+  unmark();
+  return false;
 }
 
 bool Collapse::setEdge(Entity* e)
@@ -127,8 +158,10 @@ bool checkEdgeCollapseEdgeRings(Adapt* a, Entity* edge)
   Mesh* m = a->mesh;
   Entity* v[2];
   m->getDownward(edge,0,v);
-  PCU_ALWAYS_ASSERT( ! m->isShared(v[0]));
-  PCU_ALWAYS_ASSERT( ! m->isShared(v[1]));
+  if (!getFlag(a, v[0], DONT_COLLAPSE)) //Allow collapse in one direction
+    PCU_ALWAYS_ASSERT( ! m->isShared(v[0]));
+  if (!getFlag(a, v[1], DONT_COLLAPSE)) //Allow collapse in one direction
+    PCU_ALWAYS_ASSERT( ! m->isShared(v[1]));
   apf::Up ve[2];
   m->getUp(v[0],ve[0]);
   m->getUp(v[1],ve[1]);

ma/maCollapse.h

@@ -41,6 +41,7 @@ class Collapse
     bool tryThisDirectionNoCancel(double qualityToBeat);
     bool tryBothDirections(double qualityToBeat);
     void getOldElements(EntityArray& oldElements);
+    bool edgesGoodSize();
     double getOldQuality();
     Adapt* adapt;
     Entity* edge;