Coarsen and Snapping (#514)

Changes to coarsen
- improvement to independent set
- fixed bug that made it impossible to collapse edge in both directions
- fixed bug that made it slower to collapse in one direction in parallel
- added multi-step coarsen

Changes to snapping
- fixed bug with snapping coarsening the mesh too mush
- fixed bug where not all regions were included in first problem plane
- should snap vertices in most cases
- added quality settings for split collapse

More testing
- added test for common errors in coarsen, refinement and snapping

Most of the these changes to snapping and coarsen came from Li's thesis with some minor exceptions:
- In the coarsen procedure in Li's thesis it describes getting the shortest edge next to a vertex and collapsing it. We found we had better results if we try collapsing the shortest edge and then grab the next shortest until one succeeds.
- Another difference is that Li's thesis mentions collapses, but it has very little detail in how those collapses are selected, so our logic for collapses is taken from our old mesh adapt library.
- Li's thesis also describes using FaceSwap in snapping and vertex repositioning in snapping and coarsen, but we haven't created either of those operators yet.


---------

Signed-off-by: Aiden Woodruff <[email protected]>
Co-authored-by: Nightly Bot <[email protected]>
Co-authored-by: Cameron Smith <[email protected]>
Co-authored-by: Aiden Woodruff <[email protected]>
Co-authored-by: Aditya Joshi <[email protected]>

Author: 5 people

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;