diff --git a/.github/workflows/simmetrix_enabled_pr_comment_trigger_self_hosted.yml b/.github/workflows/simmetrix_enabled_pr_comment_trigger_self_hosted.yml index 02999fdbf..1af4cd3b2 100644 --- a/.github/workflows/simmetrix_enabled_pr_comment_trigger_self_hosted.yml +++ b/.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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 621afc77b..7a2ad4c00 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/gmi_cap/gmi_cap.cc b/gmi_cap/gmi_cap.cc index 9996caed6..ab1b697e0 100644 --- a/gmi_cap/gmi_cap.cc +++ b/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; +} diff --git a/gmi_cap/gmi_cap.h b/gmi_cap/gmi_cap.h index 6ee457198..fcf9e31e2 100644 --- a/gmi_cap/gmi_cap.h +++ b/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 /** diff --git a/ma/CMakeLists.txt b/ma/CMakeLists.txt index c8af8e6fe..e54b59b38 100644 --- a/ma/CMakeLists.txt +++ b/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() diff --git a/ma/ma.cc b/ma/ma.cc index 24c5100ce..b1aa8c3f9 100644 --- a/ma/ma.cc +++ b/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); diff --git a/ma/maCoarsen.cc b/ma/maCoarsen.cc index 35089e344..9baa67130 100644 --- a/ma/maCoarsen.cc +++ b/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 +#include "apfShape.h" +#include +#include +#include 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; imesh, 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& shortEdgeVerts, std::list::iterator& itr, size_t& checked, apf::Up& adjacent, Tag* lengthTag) +{ + Entity* vertex = *itr; + std::vector 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& shortEdgeVerts, std::list::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 getShortEdgeVerts(Adapt* a, Tag* lengthTag) +{ + std::list 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 shortEdgeVerts = getShortEdgeVerts(a, lengthTag); + + Collapse collapse; + collapse.Init(a); + int success = 0; + size_t checked = 0; + bool independentSetStarted = false; + std::list::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; +} + + } diff --git a/ma/maCoarsen.h b/ma/maCoarsen.h index bf83b1e54..2e4283c77 100644 --- a/ma/maCoarsen.h +++ b/ma/maCoarsen.h @@ -14,6 +14,7 @@ namespace ma { class Adapt; +bool coarsenMultiple(Adapt* a); bool coarsen(Adapt* a); bool coarsenLayer(Adapt* a); diff --git a/ma/maCollapse.cc b/ma/maCollapse.cc index 184e7e2be..af13ba324 100644 --- a/ma/maCollapse.cc +++ b/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]); diff --git a/ma/maCollapse.h b/ma/maCollapse.h index bf70c50a3..d9a449adf 100644 --- a/ma/maCollapse.h +++ b/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; diff --git a/ma/maDBG.cc b/ma/maDBG.cc index 2460ec809..ed78f92df 100644 --- a/ma/maDBG.cc +++ b/ma/maDBG.cc @@ -30,10 +30,11 @@ static double PI = 3.14159265359; namespace ma_dbg { - +// If no dimension specified then it will write the greatest dimension to the file. void writeMesh(ma::Mesh* m, const char* prefix, - const char* suffix) + const char* suffix, + int dim=-1) { std::stringstream ss; if (std::string(suffix) != "") @@ -42,7 +43,28 @@ void writeMesh(ma::Mesh* m, ss << prefix; std::string tmp = ss.str(); const char* fileName = tmp.c_str(); - apf::writeVtkFiles(fileName, m); + apf::writeVtkFiles(fileName, m, dim); +} + +void addClassification(ma::Adapt* a, + const char* fieldName) +{ + ma::Mesh* m = a->mesh; + apf::Field* field; + field = m->findField(fieldName); + if (field) + apf::destroyField(field); + + field = apf::createFieldOn(m, fieldName, apf::SCALAR); + ma::Entity* ent; + ma::Iterator* it; + it = m->begin(0); + while ( (ent = m->iterate(it)) ){ + ma::Model* g = m->toModel(ent); + double modelDimension = (double)m->getModelType(g); + apf::setComponents(field, ent, 0, &modelDimension); + } + m->end(it); } void addTargetLocation(ma::Adapt* a, @@ -168,7 +190,7 @@ void dumpMeshWithQualities(ma::Adapt* a, ss << std::setfill('0') << std::setw(3) << iter << "_"; ss << prefix; - writeMesh(a->mesh, ss.str().c_str(), ""); + writeMesh(a->mesh, ss.str().c_str(), "", -1); apf::Field* colorField; colorField = a->mesh->findField("qual_metric"); @@ -208,7 +230,7 @@ void dumpMeshWithFlag(ma::Adapt* a, } ss << prefix << "_" << std::setfill('0') << std::setw(3) << iter; - writeMesh(a->mesh, ss.str().c_str(), ""); + writeMesh(a->mesh, ss.str().c_str(), "", dim); apf::Field* colorField; colorField = a->mesh->findField(flagName); @@ -240,7 +262,10 @@ void createCavityMesh(ma::Adapt* a, cavityMesh->acceptChanges(); std::stringstream ss; - ss << a->input->debugFolder << "/"; + //Allows the user to print to current directory or to folder + if (a->input->debugFolder) { + ss << a->input->debugFolder << "/"; + } ss << prefix; diff --git a/ma/maDBG.h b/ma/maDBG.h index 027a245fd..16a65d882 100644 --- a/ma/maDBG.h +++ b/ma/maDBG.h @@ -25,6 +25,11 @@ void writeMesh(ma::Mesh* m, const char* prefix, const char* suffix); +/* Creates a field to contain the model classification for each vertex. Which can be + printed to a file using writeMesh() with dim=0. */ +void addClassification(ma::Adapt* a, + const char* fieldName); + void addTargetLocation(ma::Adapt* a, const char* fieldName); diff --git a/ma/maDoubleSplitCollapse.cc b/ma/maDoubleSplitCollapse.cc index d0457e54c..dfd69f9c1 100644 --- a/ma/maDoubleSplitCollapse.cc +++ b/ma/maDoubleSplitCollapse.cc @@ -55,13 +55,13 @@ void DoubleSplitCollapse::accept() collapse.destroyOldElements(); } -bool DoubleSplitCollapse::run(Entity** edges) +bool DoubleSplitCollapse::run(Entity** edges, double quality) { Adapt* a = getAdapt(); Mesh* m = a->mesh; if ( ! splits.setEdges(edges,2)) return false; - oldQuality = getWorstQuality(a,splits.getTets()); + oldQuality = (quality >= 0) ? quality : getWorstQuality(a,splits.getTets()); splits.makeNewElements(); splits.transfer(); Entity* splitVerts[2]; diff --git a/ma/maDoubleSplitCollapse.h b/ma/maDoubleSplitCollapse.h index 52124378b..048faa0aa 100644 --- a/ma/maDoubleSplitCollapse.h +++ b/ma/maDoubleSplitCollapse.h @@ -24,7 +24,10 @@ class DoubleSplitCollapse bool tryThisCollapse(); bool tryBothCollapses(Entity* e); void accept(); - bool run(Entity** edges); + /* Quality is optional since during fix shape we don't need a target quality as + the quality improves. However we need a target quality during snapping since + we don't care if the quality gets worse. */ + bool run(Entity** edges, double quality = -1); Adapt* getAdapt(); private: Splits splits; diff --git a/ma/maMatch.cc b/ma/maMatch.cc index c17ecdda5..1a0582c1b 100644 --- a/ma/maMatch.cc +++ b/ma/maMatch.cc @@ -75,6 +75,9 @@ void matchNewElements(Refine* r) print(m->getPCU(), "updated matching for %li faces", face_count); } +/* we are starting to support a few operations on matched +meshes, including snapping+UR. this should prevent snapping +from modifying any matched entities */ void preventMatchedCavityMods(Adapt* a) { Mesh* m = a->mesh; diff --git a/ma/maMatchedSnapper.cc b/ma/maMatchedSnapper.cc index c5e87feb1..311c571eb 100644 --- a/ma/maMatchedSnapper.cc +++ b/ma/maMatchedSnapper.cc @@ -16,11 +16,10 @@ namespace ma { -MatchedSnapper::MatchedSnapper(Adapt* a, Tag* st, bool is) +MatchedSnapper::MatchedSnapper(Adapt* a, Tag* st) { adapter = a; snapTag = st; - isSimple = is; sharing = apf::getSharing(a->mesh); vert = 0; } @@ -38,7 +37,7 @@ void MatchedSnapper::setVert(Entity* v) { snappers.setSize(0); snappers.setSize(1); - snappers[0] = new Snapper(adapter, snapTag, isSimple); + snappers[0] = new Snapper(adapter, snapTag); snappers[0]->setVert(v); vert = v; } @@ -63,7 +62,7 @@ void MatchedSnapper::setVerts() snappers.setSize(copies.getSize() + 1); for (unsigned i = 0; i < snappers.getSize(); i++) - snappers[i] = new Snapper(adapter, snapTag, isSimple); + snappers[i] = new Snapper(adapter, snapTag); snappers[0]->setVert(v); for (unsigned i = 0; i < copies.getSize(); i++) snappers[i+1]->setVert(copies[i].entity); diff --git a/ma/maMatchedSnapper.h b/ma/maMatchedSnapper.h index c8fdbd487..42a7c12ce 100644 --- a/ma/maMatchedSnapper.h +++ b/ma/maMatchedSnapper.h @@ -21,7 +21,7 @@ namespace ma { class MatchedSnapper { public: - MatchedSnapper(Adapt* a, Tag* st, bool is); + MatchedSnapper(Adapt* a, Tag* st); ~MatchedSnapper(); void setVert(Entity* v); bool requestLocality(apf::CavityOp* o); @@ -34,7 +34,6 @@ class MatchedSnapper Entity* vert; apf::DynamicArray snappers; apf::DynamicArray locations; - bool isSimple; apf::Sharing* sharing; }; diff --git a/ma/maShape.cc b/ma/maShape.cc index f34ebcc29..dcd0eec17 100644 --- a/ma/maShape.cc +++ b/ma/maShape.cc @@ -881,4 +881,4 @@ void printQuality(Adapt* a) print(a->mesh->getPCU(), "worst element quality is %e", minqual); } -} +} \ No newline at end of file diff --git a/ma/maSingleSplitCollapse.cc b/ma/maSingleSplitCollapse.cc index 824816d42..1c326db85 100644 --- a/ma/maSingleSplitCollapse.cc +++ b/ma/maSingleSplitCollapse.cc @@ -59,7 +59,7 @@ void SingleSplitCollapse::accept() collapse.destroyOldElements(); } -bool SingleSplitCollapse::run(Entity* edge, Entity* vert) +bool SingleSplitCollapse::run(Entity* edge, Entity* vert, double quality) { oldEdge = edge; oldVert = vert; @@ -67,7 +67,7 @@ bool SingleSplitCollapse::run(Entity* edge, Entity* vert) Mesh* m = a->mesh; if ( ! splits.setEdges(&edge,1)) return false; - oldQuality = getWorstQuality(a,splits.getTets()); + oldQuality = (quality >= 0 ) ? quality : getWorstQuality(a,splits.getTets()); splits.makeNewElements(); splits.transfer(); Entity* collVerts[2]; diff --git a/ma/maSingleSplitCollapse.h b/ma/maSingleSplitCollapse.h index a0dd30620..19cefb6f2 100644 --- a/ma/maSingleSplitCollapse.h +++ b/ma/maSingleSplitCollapse.h @@ -32,7 +32,10 @@ class SingleSplitCollapse bool tryThisCollapse(); bool tryBothCollapses(Entity* e); void accept(); - bool run(Entity* edge, Entity* vert); + /* Quality is optional since during fix shape we don't need a target quality as + the quality improves. However we need a target quality during snapping since + we don't care if the quality gets worse. */ + bool run(Entity* edge, Entity* vert, double quality=-1); Adapt* getAdapt(); private: Entity *oldVert, *oldEdge; diff --git a/ma/maSize.cc b/ma/maSize.cc index 7f097d2eb..e596d0788 100644 --- a/ma/maSize.cc +++ b/ma/maSize.cc @@ -216,11 +216,11 @@ struct MetricSizeField : public SizeField } bool shouldSplit(Entity* edge) { - return this->measure(edge) > 1.5; + return this->measure(edge) > MAXLENGTH; } bool shouldCollapse(Entity* edge) { - return this->measure(edge) < 0.5; + return this->measure(edge) < MINLENGTH; } double getWeight(Entity* e) { diff --git a/ma/maSize.h b/ma/maSize.h index dd88dc643..e8f9b6058 100644 --- a/ma/maSize.h +++ b/ma/maSize.h @@ -19,6 +19,13 @@ namespace ma { typedef apf::Matrix3x3 Matrix; +/* Desired length bounds in metric space to replace hard coded values. Right now +these values are const, since we are not sure it is worth while to have the user, +modify these in addition to desired length. Future testing in need because Li's +thesis describes there might be utility to having the user modify them */ +const double MAXLENGTH = 1.5; +const double MINLENGTH = .5; +const double MAXLENGTHRATIO = 1.2; class SizeField { diff --git a/ma/maSnap.cc b/ma/maSnap.cc index 4fea96021..3c4d92b76 100644 --- a/ma/maSnap.cc +++ b/ma/maSnap.cc @@ -11,6 +11,7 @@ #include "maAdapt.h" #include "maOperator.h" #include "maSnapper.h" +#include "maRefine.h" #include "maMatchedSnapper.h" #include "maLayer.h" #include "maMatch.h" @@ -20,9 +21,24 @@ #include #include #include +#include +#include + +#ifdef PUMI_HAS_CAPSTONE +#include "gmi_cap.h" +#endif namespace ma { +static bool isCapstone(apf::Mesh* m) { + #ifdef PUMI_HAS_CAPSTONE + return gmi_cap_test(m->getModel()); + #else + (void) m; + return false; + #endif +} + static size_t isSurfUnderlyingFaceDegenerate( apf::Mesh* m, Model* g, // this the model entity in question @@ -68,10 +84,10 @@ static size_t isSurfUnderlyingFaceDegenerate( m->getFirstDerivative(g, p, uTan, vTan); double uTanSize = uTan.getLength(); double vTanSize = vTan.getLength(); -#ifdef PUMI_HAS_CAPSTONE - uTanSize = uTan * uTan; - vTanSize = vTan * vTan; -#endif + if (isCapstone(m)) { + uTanSize = uTan * uTan; + vTanSize = vTan * vTan; + } if (uTanSize < tol || vTanSize < tol) { axis = degenAxes; values.push_back(candidateDegenParam); @@ -119,6 +135,8 @@ static double interpolateParametricCoordinate( } double period = range[1]-range[0]; double span = b-a; + if (period < 0.5) //meshes in capstone can have partial ranges that aren't periodic + return (1-t)*a + t*b; if (!mode) { if (span < (period/2)) return (1-t)*a + t*b; @@ -150,46 +168,46 @@ static void interpolateParametricCoordinateOnEdge( p[1] = 0.0; p[2] = 0.0; -#ifdef PUMI_HAS_CAPSTONE - // account for non-uniform parameterization of model-edge - Vector X[3]; - Vector para[2] = {a, b}; - for (int i = 0; i < 2; i++) - m->snapToModel(g, para[i], X[i]); + if (isCapstone(m)) { + // account for non-uniform parameterization of model-edge + Vector X[3]; + Vector para[2] = {a, b}; + for (int i = 0; i < 2; i++) + m->snapToModel(g, para[i], X[i]); - double tMax = 1., tMin = 0.; - m->snapToModel(g, p, X[2]); + double tMax = 1., tMin = 0.; + m->snapToModel(g, p, X[2]); - // check if the snap point on the model edge is - // approximately at same length from either vertices + // check if the snap point on the model edge is + // approximately at same length from either vertices - double r = (X[0]-X[2]).getLength(); - double s = (X[1]-X[2]).getLength(); + double r = (X[0]-X[2]).getLength(); + double s = (X[1]-X[2]).getLength(); - double alpha = t/(1. - t); //parametric ratio + double alpha = t/(1. - t); //parametric ratio - int num_it = 0; - while (r/s < 0.95 * alpha || r/s > 1.05 * alpha) { - if ( r/s > alpha) { - tMax = t; - t = (tMin + t) / 2.; - } - else { - tMin = t; - t = (t + tMax) / 2.; - } + int num_it = 0; + while (r/s < 0.95 * alpha || r/s > 1.05 * alpha) { + if ( r/s > alpha) { + tMax = t; + t = (tMin + t) / 2.; + } + else { + tMin = t; + t = (t + tMax) / 2.; + } - p[0] = interpolateParametricCoordinate(t, a[0], b[0], range, isPeriodic, 0); + p[0] = interpolateParametricCoordinate(t, a[0], b[0], range, isPeriodic, 0); - m->snapToModel(g, p, X[2]); + m->snapToModel(g, p, X[2]); - r = (X[0]-X[2]).getLength(); - s = (X[1]-X[2]).getLength(); + r = (X[0]-X[2]).getLength(); + s = (X[1]-X[2]).getLength(); - if ( num_it > 20) break; - num_it++; + if ( num_it > 20) break; + num_it++; + } } -#endif } // convert (phi,theta) on unit sphere to (x,y,z) @@ -490,9 +508,9 @@ static void interpolateParametricCoordinatesOnRegularFace( * 1) we are assuming manifold surfaces * 2) we only check for faces that are periodic */ - -#ifndef PUMI_HAS_CAPSTONE // this need to be done for faces, only + if (isCapstone(m)) + return; if (dim != 2) return; if (!gface_isPeriodic) @@ -508,9 +526,6 @@ static void interpolateParametricCoordinatesOnRegularFace( bool isPeriodic = m->getPeriodicRange(g,d,range); p[d] = interpolateParametricCoordinate(t,a[d],b[d],range,isPeriodic, 1); } -#else - (void) gface_isPeriodic; -#endif } static void interpolateParametricCoordinatesOnFace( @@ -526,50 +541,51 @@ static void interpolateParametricCoordinatesOnFace( size_t num = isSurfUnderlyingFaceDegenerate(m, g, axes, vals); if (num > 0) { // the underlying surface is degenerate -#ifndef PUMI_HAS_CAPSTONE - interpolateParametricCoordinatesOnDegenerateFace(m, g, t, a, b, axes, vals, p); -#else - // account for non-uniform parameterization of model-edge - Vector X[3]; - Vector para[2] = {a, b}; - for (int i = 0; i < 2; i++) - m->snapToModel(g, para[i], X[i]); + if (!isCapstone(m)) { + interpolateParametricCoordinatesOnDegenerateFace(m, g, t, a, b, axes, vals, p); + } + else { + // account for non-uniform parameterization of model-edge + Vector X[3]; + Vector para[2] = {a, b}; + for (int i = 0; i < 2; i++) + m->snapToModel(g, para[i], X[i]); - interpolateParametricCoordinatesOnDegenerateFace(m, g, t, para[0], para[1], axes, vals, p); + interpolateParametricCoordinatesOnDegenerateFace(m, g, t, para[0], para[1], axes, vals, p); - double tMax = 1., tMin = 0.; - m->snapToModel(g, p, X[2]); + double tMax = 1., tMin = 0.; + m->snapToModel(g, p, X[2]); - // check if the snap point on the model edge is - // approximately at same length from either vertices + // check if the snap point on the model edge is + // approximately at same length from either vertices - double r = (X[0]-X[2]).getLength(); - double s = (X[1]-X[2]).getLength(); + double r = (X[0]-X[2]).getLength(); + double s = (X[1]-X[2]).getLength(); - double alpha = t/(1. - t); //parametric ratio + double alpha = t/(1. - t); //parametric ratio - int num_it = 0; - while (r/s < 0.95 * alpha || r/s > 1.05 * alpha) { - if ( r/s > alpha) { - tMax = t; - t = (tMin + t) / 2.; - } - else { - tMin = t; - t = (t + tMax) / 2.; - } + int num_it = 0; + while (r/s < 0.95 * alpha || r/s > 1.05 * alpha) { + if ( r/s > alpha) { + tMax = t; + t = (tMin + t) / 2.; + } + else { + tMin = t; + t = (t + tMax) / 2.; + } - interpolateParametricCoordinatesOnDegenerateFace(m, g, t, para[0], para[1], axes, vals, p); + interpolateParametricCoordinatesOnDegenerateFace(m, g, t, para[0], para[1], axes, vals, p); - m->snapToModel(g, p, X[2]); + m->snapToModel(g, p, X[2]); - r = (X[0]-X[2]).getLength(); - s = (X[1]-X[2]).getLength(); + r = (X[0]-X[2]).getLength(); + s = (X[1]-X[2]).getLength(); - if ( num_it > 20) break; - num_it++; + if ( num_it > 20) break; + num_it++; + } } -#endif } else interpolateParametricCoordinatesOnRegularFace(m, g, t, a, b, p); @@ -742,12 +758,11 @@ static void getSnapPoint(Mesh* m, Entity* v, Vector& x) class SnapAll : public Operator { public: - SnapAll(Adapt* a, Tag* t, bool simple): - snapper(a, t, simple) + SnapAll(Adapt* a, Tag* t, Snapper& snapperIn) { adapter = a; tag = t; - successCount = 0; + snapper = &snapperIn; didAnything = false; vert = 0; } @@ -756,21 +771,20 @@ class SnapAll : public Operator { if ( ! getFlag(adapter, e, SNAP)) return false; + if (getFlag(adapter, e, LAYER)) + return false; vert = e; - snapper.setVert(e); + snapper->setVert(e); return true; } bool requestLocality(apf::CavityOp* o) { - return snapper.requestLocality(o); + return snapper->requestLocality(o); } void apply() { - bool snapped = snapper.run(); - didAnything = didAnything || snapped || snapper.dug; - if (snapped) - ++successCount; - clearFlag(adapter, vert, SNAP); + bool snapped = snapper->run(); + didAnything = didAnything || snapped; } int successCount; bool didAnything; @@ -778,22 +792,21 @@ class SnapAll : public Operator Adapt* adapter; Tag* tag; Entity* vert; - Snapper snapper; + Snapper* snapper; }; -bool snapAllVerts(Adapt* a, Tag* t, bool isSimple, long& successCount) +bool snapAllVerts(Adapt* a, Tag* t, Snapper& snapper) { - SnapAll op(a, t, isSimple); + SnapAll op(a, t, snapper); applyOperator(a, &op); - successCount += a->mesh->getPCU()->Add(op.successCount); return a->mesh->getPCU()->Or(op.didAnything); } class SnapMatched : public Operator { public: - SnapMatched(Adapt* a, Tag* t, bool simple): - snapper(a, t, simple) + SnapMatched(Adapt* a, Tag* t): + snapper(a, t) { adapter = a; tag = t; @@ -832,19 +845,18 @@ class SnapMatched : public Operator MatchedSnapper snapper; }; -bool snapMatchedVerts(Adapt* a, Tag* t, bool isSimple, long& successCount) +bool snapMatchedVerts(Adapt* a, Tag* t, long& successCount) { - SnapMatched op(a, t, isSimple); + SnapMatched op(a, t); applyOperator(a, &op); successCount += a->mesh->getPCU()->Add(op.successCount); return a->mesh->getPCU()->Or(op.didAnything); } -long tagVertsToSnap(Adapt* a, Tag*& t) +long tagVertsToSnap(Adapt* a, Tag*& tag) { Mesh* m = a->mesh; int dim = m->getDimension(); - t = m->createDoubleTag("ma_snap", 3); Entity* v; long n = 0; Iterator* it = m->begin(0); @@ -857,7 +869,8 @@ long tagVertsToSnap(Adapt* a, Tag*& t) Vector x = getPosition(m, v); if (apf::areClose(s, x, 1e-12)) continue; - m->setDoubleTag(v, t, &s[0]); + m->setDoubleTag(v, tag, &s[0]); + setFlag(a, v, SNAP); if (m->isOwned(v)) ++n; } @@ -865,108 +878,62 @@ long tagVertsToSnap(Adapt* a, Tag*& t) return m->getPCU()->Add(n); } -static void markVertsToSnap(Adapt* a, Tag* t) +long snapTaggedVerts(Adapt* a, Tag* tag, Snapper& snapper) { - HasTag p(a->mesh, t); - markEntities(a, 0, p, SNAP, DONT_SNAP); -} + bool shouldForce = a->input->shouldForceAdaptation; + a->input->shouldForceAdaptation = true; //Allows quality to decrease from snapping -bool snapOneRound(Adapt* a, Tag* t, bool isSimple, long& successCount) -{ - markVertsToSnap(a, t); - if (a->mesh->hasMatching()) - return snapMatchedVerts(a, t, isSimple, successCount); - else - return snapAllVerts(a, t, isSimple, successCount); -} - -long snapTaggedVerts(Adapt* a, Tag* tag) -{ long successCount = 0; - /* there are two approaches possible here: - * 1- first snap all the vertices we can without any additional - * operation such as digging (simple snap). And then try snapping - * the remaining vertices that will need extra modifications (non- - * simple snap). - * 2- first do the non-simple snaps and then the simple snaps. - * - * Here we choose approach 2 for the following reasons - * (a) approach 2 is approximately as fast as approach 1 - * (b) the problematic snaps will be attempted as soon as possible. - * This is extremely helpful because if we wait until later on - * bringing the vert to-be-snapped to the boundary might become more - * difficult due to the change in location of neighboring verticies - * that will be snapped before the problematic vert to-be-snapped. - */ - while (snapOneRound(a, tag, false, successCount)); - while (snapOneRound(a, tag, true, successCount)); + bool snapped = true; + int prevTagged = 0; + while (snapped) { + int tagged = tagVertsToSnap(a, tag); + if (tagged == prevTagged) break; + prevTagged = tagged; + if (a->mesh->hasMatching()) + snapped = snapMatchedVerts(a, tag, successCount); + else + snapped = snapAllVerts(a, tag, snapper); + } + + a->input->shouldForceAdaptation = shouldForce; return successCount; } +int collect(Adapt* a, int val) { + return a->mesh->getPCU()->Add(val); +} + +/* + This function follows the snapping procedure specified in Li's thesis, in which verticies + will be moved directly to the model boundary, however if the movement fails due to making + an element invalid, then several operations will be attempted in order to resolve invalid + elements. +*/ void snap(Adapt* a) { - if ( ! a->input->shouldSnap) + if (!a->input->shouldSnap) return; double t0 = pcu::Time(); - Tag* tag; - /* we are starting to support a few operations on matched - meshes, including snapping+UR. this should prevent snapping - from modifying any matched entities */ + Tag* snapTag = a->mesh->createDoubleTag("ma_snap", 3); preventMatchedCavityMods(a); - long targets = tagVertsToSnap(a, tag); - long success = snapTaggedVerts(a, tag); - snapLayer(a, tag); - apf::removeTagFromDimension(a->mesh, tag, 0); - a->mesh->destroyTag(tag); - double t1 = pcu::Time(); - print(a->mesh->getPCU(), "snapped in %f seconds: %ld targets, %ld non-layer snaps", - t1 - t0, targets, success); - if (a->hasLayer) - checkLayerShape(a->mesh, "after snapping"); -} - -void visualizeGeometricInfo(Mesh* m, const char* name) -{ - Tag* dimensionTag = m->createIntTag("ma_geom_dim",1); - Tag* idTag = m->createIntTag("ma_geom_id",1); - apf::Field* field = apf::createLagrangeField(m,"ma_param",apf::VECTOR,1); - apf::Field* targetField = apf::createLagrangeField(m,"ma_target",apf::VECTOR,1); - Iterator* it = m->begin(0); - Entity* v; - while ((v = m->iterate(it))) - { - Model* c = m->toModel(v); - int dimension = m->getModelType(c); - m->setIntTag(v,dimensionTag,&dimension); - int id = m->getModelTag(c); - m->setIntTag(v,idTag,&id); - Vector p; - Vector xp = getPosition(m, v); - m->getParam(v,p); - if (dimension == 2 || dimension == 1) { - Vector x; - m->isParamPointInsideModel(c, &p[0], x); - apf::setVector(targetField, v, 0, x - xp); - } - else { - Vector x(0., 0., 0.); - apf::setVector(targetField, v, 0, x); - } - apf::setVector(field,v,0,p); - } - m->end(it); - apf::writeVtkFiles(name,m); - it = m->begin(0); - while ((v = m->iterate(it))) + int toSnap = tagVertsToSnap(a, snapTag); + if (toSnap) { - m->removeTag(v,dimensionTag); - m->removeTag(v,idTag); + Snapper snapper(a, snapTag); + snapTaggedVerts(a, snapTag, snapper); + snapLayer(a, snapTag); + + double t1 = pcu::Time(); + print(a->mesh->getPCU(), "ToSnap %d - Moved %d - Failed %d - CollapseToVtx %d" + " - Collapse %d - Swap %d - SplitCollapse %d" + " - completed in %f seconds", + toSnap, collect(a,snapper.numSnapped), collect(a,snapper.numFailed), + collect(a,snapper.numCollapseToVtx), collect(a,snapper.numCollapse), + collect(a,snapper.numSwap), collect(a,snapper.numSplitCollapse), t1 - t0); + if (a->hasLayer) + checkLayerShape(a->mesh, "after snapping"); } - m->end(it); - m->destroyTag(dimensionTag); - m->destroyTag(idTag); - apf::destroyField(field); - apf::destroyField(targetField); + a->mesh->destroyTag(snapTag); } - } diff --git a/ma/maSnapper.cc b/ma/maSnapper.cc index 14c9a1651..2f3a426d2 100644 --- a/ma/maSnapper.cc +++ b/ma/maSnapper.cc @@ -7,10 +7,19 @@ of the SCOREC Non-Commercial License this program is distributed under. *******************************************************************************/ +/** + * \file maSnapper.cc + * \brief Definition of maSnapper.h file. + * This file contains functions to move a point to the model surface. As described + * in Li's thesis it will first try to collapse in the target direction. Otherwise + * it will collapse to simplify the region and attempt other operators such as + * swap, split collapse, double split collapse. +*/ #include "maSnapper.h" #include "maAdapt.h" #include "maShapeHandler.h" #include "maSnap.h" +#include "maDBG.h" #include #include #include @@ -18,16 +27,20 @@ namespace ma { -Snapper::Snapper(Adapt* a, Tag* st, bool is) +Snapper::Snapper(Adapt* a, Tag* st) : mesh(a->mesh), splitCollapse(a), doubleSplitCollapse(a) { - adapter = a; + adapt = a; snapTag = st; collapse.Init(a); - isSimple = is; - dug = false; + edgeSwap = makeEdgeSwap(a); vert = 0; } +Snapper::~Snapper() +{ + delete edgeSwap; +} + void Snapper::setVert(Entity* v) { vert = v; @@ -42,241 +55,589 @@ bool Snapper::requestLocality(apf::CavityOp* o) { if (!o->requestLocality(&vert, 1)) return false; - if (isSimple) - return true; /* in order to try an edge collapse (we don't yet know which edge), bring in a cavity such that all adjacent edges have both vertices local. This is basically two layers of elements around the vertex */ apf::Up edges; - adapter->mesh->getUp(vert,edges); + mesh->getUp(vert,edges); apf::Up ovs; ovs.n = edges.n; for (int i = 0; i < edges.n; ++i) - ovs.e[i] = apf::getEdgeVertOppositeVert(adapter->mesh, edges.e[i], vert); + ovs.e[i] = apf::getEdgeVertOppositeVert(mesh, edges.e[i], vert); return o->requestLocality(&ovs.e[0], ovs.n); } -static void computeNormals(Mesh* m, Upward& es, apf::NewArray& normals) +//Write snapping data to files for debugging purposes +//In order to view relevant information it is neccessary to hide entities with relevent flag in vtk viewer +#if defined(DEBUG_FPP) +static void flagAndPrint(Adapt* a, Entity* ent, int dim, const char* name) { - if (m->getDimension() != 2) - return; - normals.allocate(es.getSize()); - for (size_t i = 0; i < es.getSize(); ++i) - normals[i] = getTriNormal(m, es[i]); + setFlag(a, ent, CHECKED); + ma_dbg::dumpMeshWithFlag(a, 0, dim, CHECKED, name, name); + clearFlag(a, ent, CHECKED); } +#endif -static bool didInvert(Mesh* m, Vector& oldNormal, Entity* tri) +//Write snapping data to files for debugging purposes +#if defined(DEBUG_FPP) +static void printFPP(Adapt* a, FirstProblemPlane* FPP) { - return (oldNormal * getTriNormal(m, tri)) < 0; -} + ma_dbg::addTargetLocation(a, "snap_target"); + ma_dbg::addClassification(a, "classification"); -static void collectBadElements(Adapt* a, Upward& es, - apf::NewArray& normals, apf::Up& bad) -{ - Mesh* m = a->mesh; - bad.n = 0; - for (size_t i = 0; i < es.getSize(); ++i) { -/* for now, when snapping a vertex on the boundary - layer, ignore the quality of layer elements. - not only do we not have metrics for this, but the - algorithm that moves curves would need to change */ - if (getFlag(a, es[i], LAYER)) - continue; - double quality = a->shape->getQuality(es[i]); - if (quality < a->input->validQuality) - bad.e[bad.n++] = es[i]; -/* check for triangles whose normals have changed by - more than 90 degrees when the vertex was snapped */ - else if ((m->getDimension() == 2) && - didInvert(m, normals[i], es[i])) - bad.e[bad.n++] = es[i]; + apf::writeVtkFiles("FPP_Mesh", a->mesh); + EntityArray invalid; + for (int i=0; iproblemRegions.n; i++){ + invalid.append(FPP->problemRegions.e[i]); } - PCU_ALWAYS_ASSERT(bad.n < (int)(sizeof(bad.e) / sizeof(Entity*))); + ma_dbg::createCavityMesh(a, invalid, "FPP_Invalid"); + + for (int i=0; icommEdges.n; i++) setFlag(a, FPP->commEdges.e[i], CHECKED); + ma_dbg::dumpMeshWithFlag(a, 0, 1, CHECKED, "FPP_CommEdges", "FPP_CommEdges"); + for (int i=0; icommEdges.n; i++) clearFlag(a, FPP->commEdges.e[i], CHECKED); + + flagAndPrint(a, FPP->vert, 0, "FPP_Vertex"); + flagAndPrint(a, FPP->problemFace, 2, "FPP_Face"); + flagAndPrint(a, FPP->problemRegion, 3, "FPP_Region"); } +#endif -static bool trySnapping(Adapt* adapter, Tag* tag, Entity* vert, - apf::Up& badElements) +//returns the greater index in the case of equality +static int indexOfMin(double a0, double a1, double a2) { - Mesh* mesh = adapter->mesh; - Vector x = getPosition(mesh, vert); - Vector s; - mesh->getDoubleTag(vert, tag, &s[0]); -/* gather the adjacent elements */ - Upward elements; - mesh->getAdjacent(vert, mesh->getDimension(), elements); -/* in 2D, get the old triangle normals */ - apf::NewArray normals; - computeNormals(mesh, elements, normals); -/* move the vertex to the desired point */ - mesh->setPoint(vert, 0, s); -/* check resulting cavity */ - collectBadElements(adapter, elements, normals, badElements); - if (badElements.n) { - /* not ok, put the vertex back where it was */ - mesh->setPoint(vert, 0, x); - return false; + if (a1 < a0) { + if (a2 < a1) return 2; + else return 1; } else { - /* ok, take off the snap tag */ - mesh->removeTag(vert, tag); - return true; + if (a2 < a0) return 2; + else return 0; } } -static bool tryDiggingEdge(Adapt* adapter, Collapse& collapse, Entity* e) +static Vector projOnTriPlane(Adapt* a, Entity* vert, Vector normal, Vector v0) { - Mesh* mesh = adapter->mesh; - PCU_ALWAYS_ASSERT(mesh->getType(e) == apf::Mesh::EDGE); - if ( ! collapse.setEdge(e)) - return false; - if ( ! collapse.checkClass()) - return false; - if ( ! collapse.checkTopo()) - return false; - double q = adapter->input->validQuality; - bool oldShouldForce = adapter->input->shouldForceAdaptation; - adapter->input->shouldForceAdaptation = true; - if ( ! collapse.tryBothDirections(q)) - return false; - adapter->input->shouldForceAdaptation = oldShouldForce; - collapse.destroyOldElements(); - return true; + double magN = normal*normal; + Vector vertPos = getPosition(a->mesh, vert); + double magCP = (vertPos-v0) * normal; + double ratio=magCP/magN; + + Vector result; + for (int i=0; i<3; ++i) + result[i]=vertPos[i]-ratio*normal[i]; + + return result; } -static bool tryDigging2(Adapt* a, Collapse& c, apf::Up& badElements, - FirstProblemPlane* FPP) +/* + Given a poorly-shaped tetrahedron, a base triangle and the opposite vertex of the base, + determine the following information: + 1. the key mesh entities to apply local mesh modification + 2. area of the four face + + return 0 : if an edge is degenerated + 3,5,6 : the tetrahedron has two large dihedral angles. The opposite edges will be stored in ents[0], ents[1]. + 1,2,4,7 : the tetrahedron has three large angles. The largeest face is stored in ents[0]. +*/ +int getTetStats(Adapt* a, Entity* vert, Entity* face, Entity* region, Entity* ents[4], double area[4]) { + Entity* faceEdges[3]; + a->mesh->getDownward(face, 1, faceEdges); - Mesh* m = a->mesh; - int dim = m->getDimension(); + Entity* verts[3]; + a->mesh->getDownward(face, 0, verts); - // first go through the candidate edges found by the first problem plane - // (if any) - std::vector edgesFromFPP; - edgesFromFPP.clear(); - if (FPP && dim == 3) { - FPP->setBadElements(badElements); - FPP->getCandidateEdges(edgesFromFPP); + Vector facePos[3]; + Entity* edges[6]; + for (int i=0; i<3; i++) { + edges[i]=faceEdges[i]; + facePos[i]=getPosition(a->mesh, verts[i]); } - for (size_t i = 0; i < edgesFromFPP.size(); ++i) { - if (tryDiggingEdge(a, c, edgesFromFPP[i])) + + Entity* faces[4]; + faces[0]=face; + + Entity* problemFaces[4]; + a->mesh->getDownward(region, 2, problemFaces); + for (int i=0; i<4; i++) { + if (problemFaces[i] == face ) continue; + else if (isLowInHigh(a->mesh, problemFaces[i], edges[0])) faces[1] = problemFaces[i]; + else if (isLowInHigh(a->mesh, problemFaces[i], edges[1])) faces[2] = problemFaces[i]; + else if (isLowInHigh(a->mesh, problemFaces[i], edges[2])) faces[3] = problemFaces[i]; + } + + for (int i=1; i<3; i++) { + Entity* problemEdges[3]; + a->mesh->getDownward(faces[i], 1, problemEdges); + for (int j=0; j<3; j++) { + if (problemEdges[j]==edges[i-1] ) continue; + else if (isLowInHigh(a->mesh, problemEdges[j], verts[0])) edges[3] = problemEdges[j]; + else if (isLowInHigh(a->mesh, problemEdges[j], verts[1])) edges[4] = problemEdges[j]; + else if (isLowInHigh(a->mesh, problemEdges[j], verts[2])) edges[5] = problemEdges[j]; + } + } + + /* find normal to the plane */ + Vector v01 = facePos[1] - facePos[0]; + Vector v02 = facePos[2] - facePos[0]; + Vector norm = apf::cross(v01, v02); + + Vector projection = projOnTriPlane(a, vert, norm, facePos[0]); + Vector ri = projection - facePos[0]; + Vector rj = projection - facePos[1]; + Vector rk = projection - facePos[2]; + + /* determine which side of the edges does the point R lie. + First get normal vectors */ + Vector normi = apf::cross(v01, ri); + Vector normj = apf::cross(facePos[2]-facePos[1], rj); + Vector normk = apf::cross(facePos[0]-facePos[2], rk); + + Vector mag; + mag[0]=normi*norm; + mag[1]=normj*norm; + mag[2]=normk*norm; + + area[0]=norm*norm; + area[1]=normi*normi; + area[2]=normj*normj; + area[3]=normk*normk; + + int filter[]={1,2,4}; + int bit=0; + /* examine signs of mag[0], mag[1] and mag[2] */ + for(int i=0; i<3; i++) + if(mag[i]>0.0) + bit = bit | filter[i]; + + /* + 010=2 | 011=3 / 001=1 + | / + ------------+--e2--+----------- + v0| /v2 + | 7 / + 110=6 e0 e1 + | / + | / 101=5 + |/ + v1+ + /| + / | + 4 + */ + + switch( bit ) { + case 1:{ + int Emap[]={0,4,3}; + int Fmap[]={0,2,3}; + ents[0]=faces[1]; + int i=indexOfMin(area[0],area[2],area[3]); + ents[1]=edges[Emap[i]]; + ents[2]=faces[Fmap[i]]; + break; + } + case 2: { + int Emap[]={1,4,5}; + int Fmap[]={0,1,3}; + ents[0]=faces[2]; + int i=indexOfMin(area[0],area[1],area[3]); + ents[1]=edges[Emap[i]]; + ents[2]=faces[Fmap[i]]; + break; + } + case 3: { + ents[0]=edges[2]; + ents[1]=edges[4]; + ents[2]=((area[0]mesh->getPCU(), "Swap warning: This swap/splt may not work consider more collapses"); + } + return bit; +} + +/* + We perform this last to make sure that we have a simple region where we can determine the + best operation to perform and because we want to avoid creating more vertices to snap since + we could get stuck in an infinite loop of creating and snapping those vertices. +*/ +bool Snapper::trySwapOrSplit(FirstProblemPlane* FPP) +{ + Entity* ents[4] = {0}; + double area[4]; + int bit = getTetStats(adapt, FPP->vert, FPP->problemFace, FPP->problemRegion, ents, area); + + double min=area[0]; + for(int i=1; i<4; i++) + if( area[i]getDownward(FPP->problemFace, 1, edges); + Entity* longest = edges[0]; + for (int i=1; i<3; i++) + if (adapt->sizeField->measure(edges[i]) > adapt->sizeField->measure(longest)) + longest = edges[i]; + + if (edgeSwap->run(longest)) { + numSwap++; + return true; + } + if (splitCollapse.run(longest, FPP->vert, adapt->input->validQuality)) { + numSplitCollapse++; return true; + } } - // next try all the edges - for (int i = 0; i < badElements.n; ++i) { - Entity* elem = badElements.e[i]; - Downward edges; - int nedges; - nedges = m->getDownward(elem, 1, edges); - for (int j = 0; j < nedges; ++j){ - if (tryDiggingEdge(a, c, edges[j])) - return true; + if (ents[0] == 0) + return false; + + // two large dihedral angles -> key problem: two mesh edges + if (bit==3 || bit==5 || bit==6) { + for (int i=0; i<2; i++) + if (edgeSwap->run(ents[i])) { + numSwap++; + return true; + } + for (int i=0; i<2; i++) + if (splitCollapse.run(ents[i], FPP->vert, adapt->input->validQuality)) { + numSplitCollapse++; + return true; + } + if (doubleSplitCollapse.run(ents, adapt->input->validQuality)) { + numSplitCollapse++; + return true; + } + print(mesh->getPCU(), "Swap failed: Consider more collapses"); + } + // three large dihedral angles -> key entity: a mesh face + else { + Entity* edges[3]; + mesh->getDownward(ents[0], 1, edges); + for (int i=0; i<3; i++) { + if (edgeSwap->run(edges[i])) { + numSwap++; + return true; + } + } + //TODO: RUN FACE SWAP HERE + if (splitCollapse.run(ents[1], FPP->vert, adapt->input->validQuality)) { + numSplitCollapse++; + return true; } + print(mesh->getPCU(), "Swap failed: face swap not implemented"); } return false; } -static void updateVertexParametricCoords( - Mesh* m, - Entity* vert, - Vector& newTarget) +static bool tryCollapseEdge(Adapt* a, Entity* edge, Entity* keep, Collapse& collapse) { - PCU_ALWAYS_ASSERT_VERBOSE(m->getType(vert) == apf::Mesh::VERTEX, - "expecting a vertex!"); + 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); + + bool result = false; + if (collapse.setEdge(edge) && + collapse.checkClass() && + collapse.checkTopo() && + collapse.tryBothDirections(a->input->validQuality)) { + collapse.destroyOldElements(); + result = true; + } + if (!alreadyFlagged) clearFlag(a, keep, DONT_COLLAPSE); + return result; +} - // if vert is classified on a model vert or edge return - Model* g = m->toModel(vert); - if (m->getModelType(g) != 2) - return; +struct BestCollapse +{ + double quality=-1; + Entity* edge; + Entity* keep; +}; + +/* + Perfomes a collapse operation and stores the operation in best if the quality is better, + then cancels the collapse. We want to pick the highest quality after collapsing to the + first problem plane so future operations are more likely to succeed. +*/ +static void getBestQualityCollapse(Adapt* a, Entity* edge, Entity* keep, Collapse& collapse, BestCollapse& best) +{ + 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); + if (collapse.setEdge(edge) && collapse.checkClass() && collapse.checkTopo()) { + collapse.computeElementSets(); + if (collapse.tryThisDirectionNoCancel(a->input->validQuality) && collapse.edgesGoodSize()) { + double quality = getWorstQuality(a, collapse.newElements); + if (quality > best.quality) { + best.quality = quality; + best.edge = edge; + best.keep = keep; + } + } + collapse.cancel(); + } + if (!alreadyFlagged) clearFlag(a, keep, DONT_COLLAPSE); +} - // get the list of upward adj edges that are - // classified on the same model face as vert - apf::Up edges; - m->getUp(vert,edges); - apf::Up oes; - oes.n = edges.n; - int counter = 0; - for (int i = 0; i < edges.n; ++i) { - Model* h = m->toModel(edges.e[i]); - if (m->getModelType(h) == 3) - continue; - PCU_ALWAYS_ASSERT_VERBOSE(g == h, - "expecting the model to be the same for current edge and vert"); - oes.e[counter] = edges.e[i]; - counter++; +//returns if testVert and refVert are on the same side of the face +static bool sameSide(Adapt* a, Entity* testVert, Entity* refVert, Entity* face) +{ + Entity* faceVert[3]; + a->mesh->getDownward(face, 0, faceVert); + Vector facePos[3]; + for (int i=0; i < 3; ++i) + facePos[i] = getPosition(a->mesh,faceVert[i]); + + Vector normal = apf::cross((facePos[1]-facePos[0]),(facePos[2]-facePos[0])); + Vector testPos = getPosition(a->mesh, testVert); + Vector refPos = getPosition(a->mesh, refVert); + const double tol=1e-12; + + double dr = (testPos - facePos[0]) * normal; + if (dr*dr < tol) return false; //testVert is on the face + double ds = (refPos - facePos[0]) * normal; + if (dr*ds < 0.0) return false; //different sides of face + return true; //same side of face +} + +/* + If collapsing the common edges failed we want to try collapsing any edge that will + move us towards the first problem plane. We try collapses first in order to simplify + the region until we can perform smarter operations. +*/ +bool Snapper::tryCollapseTetEdges(FirstProblemPlane* FPP) +{ + std::vector& commEdges = FPP->commEdges; + BestCollapse best; + + for (size_t i=0; igetDownward(commEdges[i], 0, vertex); + for (int j=0; j<2; j++) + getBestQualityCollapse(adapt, commEdges[i], vertex[j], collapse, best); } - Vector pBar(0., 0., 0.); - for (int i = 0; i < counter; i++) { - Vector pTmp; - transferParametricOnEdgeSplit(m, oes.e[i], 0.5, pTmp); - pBar += pTmp; + for (size_t i=0; igetUp(vertexFPP, adjEdges); + for (int j=0; jproblemFace, edgeDel)) continue; + Entity* vertKeep = getEdgeVertOppositeVert(mesh, edgeDel, vertexFPP); + if (sameSide(adapt, vertKeep, vert, FPP->problemFace)) continue; + getBestQualityCollapse(adapt, edgeDel, vertKeep, collapse, best); + } } - pBar = pBar / oes.n; - m->snapToModel(m->toModel(vert), pBar, newTarget); - m->setParam(vert, pBar); + if (best.quality > 0) { + numCollapse++; + return tryCollapseEdge(adapt, best.edge, best.keep, collapse); + } + else return false; } -static bool tryMoving(Adapt* adapter, Entity* v, Tag* tag) +/* + If collapsing to the first problem plane has failed then we want + to collapse edges on the first problem plane in order to simplify + the region future operations are more likely to succeed. +*/ +bool Snapper::tryReduceCommonEdges(FirstProblemPlane* FPP) { - Mesh* m = adapter->mesh; - PCU_ALWAYS_ASSERT_VERBOSE(m->hasTag(v, tag), - "expecting the vertex to have a tag!"); - bool hadItBefore = getFlag(adapter, v, DONT_MOVE); - setFlag(adapter, v, DONT_MOVE); - Vector newTarget; - m->getDoubleTag(v, tag, &newTarget[0]); // default - updateVertexParametricCoords(m, v, newTarget); - m->setDoubleTag(v, tag, &newTarget[0]); - if (!hadItBefore) - clearFlag(adapter, v, DONT_MOVE); - return true; + std::vector& commEdges = FPP->commEdges; + BestCollapse best; + + Entity* pbEdges[3]; + mesh->getDownward(FPP->problemFace, 1, pbEdges); + switch(commEdges.size()) { + case 2: { + Entity* v1 = getEdgeVertOppositeVert(mesh, commEdges[0], vert); + Entity* v2 = getEdgeVertOppositeVert(mesh, commEdges[1], vert); + + for (int i=0; i<3; i++) { + Entity* pbVert[2]; + mesh->getDownward(pbEdges[i], 0, pbVert); + if (pbVert[0] == v1 && pbVert[1] == v2) continue; + if (pbVert[1] == v1 && pbVert[0] == v2) continue; + for (int j=0; j<2; j++) + getBestQualityCollapse(adapt, pbEdges[i], pbVert[j], collapse, best); + } + break; + } + case 3: { + for (int i=0; i<3; i++) { + Entity* pbVert[2]; + mesh->getDownward(pbEdges[i], 0, pbVert); + for (int j=0; j<2; j++) + getBestQualityCollapse(adapt, pbEdges[i], pbVert[j], collapse, best); + } + break; + } + } + if (best.quality > 0) { + numCollapse++; + return tryCollapseEdge(adapt, best.edge, best.keep, collapse); + } + else return false; } -static bool tryDigging(Adapt* a, Collapse& c, Entity* v, - apf::Up& badElements, FirstProblemPlane* FPP = 0) +/* + First we try collapsing to a vertex on the first problem plane because this is the most likely + operation to succeed. We first try collapsing to the common edges among the invalid regions + since those are more likely to succeed. +*/ +bool Snapper::tryCollapseToVertex(FirstProblemPlane* FPP) { - bool hadItBefore = getFlag(a, v, DONT_COLLAPSE); - setFlag(a, v, DONT_COLLAPSE); - bool ok = tryDigging2(a, c, badElements, FPP); - if (!hadItBefore) - clearFlag(a, v, DONT_COLLAPSE); - return ok; + Vector position = getPosition(mesh, vert); + Vector target; + mesh->getDoubleTag(vert, snapTag, &target[0]); + double distTarget = (position - target).getLength(); + + BestCollapse best; + + for (size_t i = 0; i < FPP->commEdges.size(); ++i) { + Entity* edge = FPP->commEdges[i]; + Entity* vertexOnFPP = getEdgeVertOppositeVert(mesh, edge, vert); + Vector vFPPCoord = getPosition(mesh, vertexOnFPP); + double distToFPPVert = (vFPPCoord - target).getLength(); + if (distToFPPVert > distTarget) continue; + getBestQualityCollapse(adapt, edge, vert, collapse, best); + } + + if (best.quality > 0) { + numCollapseToVtx++; + return tryCollapseEdge(adapt, best.edge, best.keep, collapse); + } + else return false; +} + +static FirstProblemPlane* getFPP(Adapt* a, Entity* vertex, Tag* snapTag, apf::Up& invalid) +{ + FirstProblemPlane* FPP = new FirstProblemPlane(a, snapTag); + FPP->setVertex(vertex); + FPP->setBadElements(invalid); + std::vector commEdges; + FPP->getCandidateEdges(commEdges); + return FPP; +} + +static void getInvalid(Adapt* a, Upward& adjacentElements, apf::Up& invalid) +{ + invalid.n = 0; + for (size_t i = 0; i < adjacentElements.getSize(); ++i) { + /* for now, when snapping a vertex on the boundary + layer, ignore the quality of layer elements. + not only do we not have metrics for this, but the + algorithm that moves curves would need to change */ + if (getFlag(a, adjacentElements[i], LAYER)) + continue; + if (a->shape->getQuality(adjacentElements[i]) < a->input->validQuality) + invalid.e[invalid.n++] = adjacentElements[i]; + } +} + +//Moved vertex to model surface or returns invalid elements if not possible +static bool tryReposition(Adapt* adapt, Entity* vertex, Tag* snapTag, apf::Up& invalid) +{ + Mesh* mesh = adapt->mesh; + if (!mesh->hasTag(vertex, snapTag)) return true; + Vector prev = getPosition(mesh, vertex); + Vector target; + mesh->getDoubleTag(vertex, snapTag, &target[0]); + Upward adjacentElements; + mesh->getAdjacent(vertex, mesh->getDimension(), adjacentElements); + + mesh->setPoint(vertex, 0, target); + getInvalid(adapt, adjacentElements, invalid); + if (invalid.n == 0) return true; + mesh->setPoint(vertex, 0, prev); + return false; } bool Snapper::trySimpleSnap() { - apf::Up badElements; - return trySnapping(adapter, snapTag, vert, badElements); + apf::Up invalid; + return tryReposition(adapt, vert, snapTag, invalid); } +#if defined(DEBUG_FPP) +static int DEBUGFAILED=0; +#endif +/* +This function will attempt to move vert to the model surface, if it can not do so then it +will atleast move to the first problem plane as described in Li's thesis. It might take multiple +iterations for vert to reach the model surface. Li's thesis was missing some details on how +to apply certain opperators so other algoritms were adapted from old scorec libraries. +*/ bool Snapper::run() { - dug = false; - moved = false; - apf::Up badElements; - bool ok = trySnapping(adapter, snapTag, vert, badElements); - if (isSimple) - return ok; - // there is no need for the following if there exists no bad elements - if (badElements.n == 0) return true; - FirstProblemPlane* FPP; -#ifdef DO_FPP - FPP = new FirstProblemPlane(adapter, snapTag); - FPP->setVertex(vert); -#else - FPP = 0; -#endif - if (adapter->mesh->getDimension() == 2) - moved = tryMoving(adapter, vert, snapTag); - else - dug = tryDigging(adapter, collapse, vert, badElements, FPP); - delete FPP; - if (!dug && !moved) - return false; - return trySnapping(adapter, snapTag, vert, badElements); + apf::Up invalid; + bool success = tryReposition(adapt, vert, snapTag, invalid); + if (success) { + numSnapped++; + mesh->removeTag(vert,snapTag); + clearFlag(adapt, vert, SNAP); + return true; + } + + FirstProblemPlane* FPP=0; + if (mesh->getDimension() == 3) { + if (!success) FPP = getFPP(adapt, vert, snapTag, invalid); + if (!success) success = tryCollapseToVertex(FPP); + if (!success) success = tryReduceCommonEdges(FPP); + if (!success) success = tryCollapseTetEdges(FPP); + if (!success) success = trySwapOrSplit(FPP); + } + + if (!success) { + numFailed++; + mesh->removeTag(vert,snapTag); + clearFlag(adapt, vert, SNAP); + } + #if defined(DEBUG_FPP) + if (!success && ++DEBUGFAILED == 1) printFPP(adapt, FPP); + #endif + if (FPP) delete FPP; + return success; } FirstProblemPlane::FirstProblemPlane(Adapt* a, Tag* st) @@ -285,7 +646,7 @@ FirstProblemPlane::FirstProblemPlane(Adapt* a, Tag* st) snapTag = st; problemFace = 0; problemRegion = 0; - commEdges.n = 0; + commEdges.clear(); tol = 1.0e-14; } @@ -296,9 +657,8 @@ void FirstProblemPlane::setVertex(Entity* v) void FirstProblemPlane::setBadElements(apf::Up& badElements) { - problemRegions.n = badElements.n; for (int i = 0; i < badElements.n; i++) { - problemRegions.e[i] = badElements.e[i]; + problemRegions.push_back(badElements.e[i]); } } @@ -318,8 +678,7 @@ bool FirstProblemPlane::find() // determine distances to all possible problem faces, the shortest // distance and its intersection on first problem plane (FPP) - int n; - n = problemRegions.n; + size_t n = problemRegions.size(); Entity* elem; Entity* face; Ray ray; @@ -330,8 +689,8 @@ bool FirstProblemPlane::find() ray.dir = target - ray.start; dists.clear(); - for (int i = 0; i < n; i++) { - elem = problemRegions.e[i]; + for (size_t i = 0; i < n; i++) { + elem = problemRegions[i]; face = getTetFaceOppositeVert(mesh, elem, vert); std::vector coords; getFaceCoords(mesh, face, coords); @@ -345,16 +704,16 @@ bool FirstProblemPlane::find() lion_oprint(1, "Info: Found Infinitely Many Intersection Points!\n"); Vector newDirection = intersect - ray.start; if (newDirection.getLength() < minDist) { - dists.push_back(newDirection.getLength()); - minDist = dists.back(); - problemFace = face; - problemRegion = elem; - intersection = intersect; - // do not need to check whether the move is valid since the valid - // ones should have been taken care of by this point + dists.push_back(newDirection.getLength()); + minDist = dists.back(); + problemFace = face; + problemRegion = elem; + intersection = intersect; + // do not need to check whether the move is valid since the valid + // ones should have been taken care of by this point } else - dists.push_back(minDist + 1.0 + tol); + dists.push_back(newDirection.getLength()); } } @@ -363,19 +722,19 @@ bool FirstProblemPlane::find() coplanarProblemRegions.n = 0; if (!problemRegion) { - problemRegion = problemRegions.e[0]; + problemRegion = problemRegions[0]; problemFace = getTetFaceOppositeVert(mesh, problemRegion, vert); coplanarProblemRegions.n = n; - for (int i = 0; i < n; i++) { - coplanarProblemRegions.e[i] = problemRegions.e[i]; + for (size_t i = 0; i < n; i++) { + coplanarProblemRegions.e[i] = problemRegions[i]; } } else { minDist += tol; - for (int i = 0; i < n; i++) { + for (size_t i = 0; i < n; i++) { if (dists[i] < minDist) { - coplanarProblemRegions.e[coplanarProblemRegions.n] = problemRegions.e[i]; - coplanarProblemRegions.n++; + coplanarProblemRegions.e[coplanarProblemRegions.n] = problemRegions[i]; + coplanarProblemRegions.n++; } } } @@ -401,8 +760,8 @@ void FirstProblemPlane::findCandidateEdges(std::vector &edges) Entity* edge; Entity* v; - for (int i = 0; i < commEdges.n; i++) { - edge = commEdges.e[i]; + for (size_t i = 0; i < commEdges.size(); i++) { + edge = commEdges[i]; Downward dv; mesh->getDownward(edge, 0, dv); (dv[0] == vert) ? v = dv[1] : v = dv[0]; @@ -476,10 +835,8 @@ void FirstProblemPlane::findCommonEdges(apf::Up& cpRegions) Downward edges; int nDownEdges = mesh->getDownward(cpRegions.e[0], 1, edges); for (int i = 0; i < nDownEdges; i++) { - if (isLowInHigh(mesh, edges[i], vert)) { - commEdges.e[commEdges.n] = edges[i]; - commEdges.n++; - } + if (isLowInHigh(mesh, edges[i], vert)) + commEdges.push_back(edges[i]); } return; } @@ -517,12 +874,9 @@ void FirstProblemPlane::findCommonEdges(apf::Up& cpRegions) if (!isLowInHigh(mesh, region, edges[i])) { flag = 0; break; - } - if (flag) { - commEdges.e[commEdges.n] = edges[i]; - commEdges.n++; - } + } } + if (flag) commEdges.push_back(edges[i]); } } } diff --git a/ma/maSnapper.h b/ma/maSnapper.h index 9711d5e65..8aa17f5da 100644 --- a/ma/maSnapper.h +++ b/ma/maSnapper.h @@ -11,6 +11,9 @@ #define MA_SNAPPER_H #include "maCollapse.h" +#include "maSingleSplitCollapse.h" +#include "maDoubleSplitCollapse.h" +#include "maEdgeSwap.h" namespace apf { class CavityOp; @@ -27,20 +30,34 @@ class FirstProblemPlane; class Snapper { public: - Snapper(Adapt* a, Tag* st, bool is); + int numFailed = 0; + int numSnapped = 0; + int numCollapseToVtx = 0; + int numCollapse = 0; + int numSwap = 0; + int numSplitCollapse = 0; + + Snapper(Adapt* a, Tag* st); + ~Snapper(); void setVert(Entity* v); Entity* getVert(); bool requestLocality(apf::CavityOp* o); bool trySimpleSnap(); bool run(); - bool dug; - bool moved; private: - Adapt* adapter; + Adapt* adapt; + Mesh* mesh; Tag* snapTag; Entity* vert; Collapse collapse; - bool isSimple; + SingleSplitCollapse splitCollapse; + DoubleSplitCollapse doubleSplitCollapse; + EdgeSwap* edgeSwap; + + bool tryCollapseToVertex(FirstProblemPlane* FPP); + bool tryCollapseTetEdges(FirstProblemPlane* FPP); + bool tryReduceCommonEdges(FirstProblemPlane* FPP); + bool trySwapOrSplit(FirstProblemPlane* FPP); }; struct Ray{ @@ -55,15 +72,15 @@ class FirstProblemPlane void setVertex(Entity* v); void setBadElements(apf::Up& badElements); void getCandidateEdges(std::vector &edges); - private: - Adapt* adapter; - Tag* snapTag; Entity* vert; - apf::Up problemRegions; Entity* problemFace; + std::vector commEdges; + std::vector problemRegions; Entity* problemRegion; + Tag* snapTag; + private: + Adapt* adapter; Vector intersection; - apf::Up commEdges; double tol; bool find(); void findCandidateEdges(std::vector &edges); @@ -72,6 +89,7 @@ class FirstProblemPlane void findCommonEdges(apf::Up& cpRegions); }; +int getTetStats(Adapt* a, Entity* vert, Entity* face, Entity* region, Entity* ents[4], double area[4]); Entity* getTetFaceOppositeVert(Mesh* m, Entity* e, Entity* v); void getFaceCoords(Mesh* m, Entity* face, std::vector& coords); Vector getCenter(Mesh* mesh, Entity* face); diff --git a/pumi-meshes b/pumi-meshes index 27dca7523..93d3dadf4 160000 --- a/pumi-meshes +++ b/pumi-meshes @@ -1 +1 @@ -Subproject commit 27dca75239ed386594b7eb35d69c2d58fb28cfe3 +Subproject commit 93d3dadf4ecf73f6ec50bdf6f91a4d3ffd0f6f10 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f6d2ec6aa..9f13a3f2a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -174,6 +174,7 @@ test_exe_func(pumiLoadMesh pumiLoadMesh.cc) test_exe_func(xgc_split xgc_split.cc) test_exe_func(ma_insphere ma_insphere.cc) test_exe_func(ma_test ma_test.cc) +test_exe_func(aniso_adapt aniso_adapt.cc) test_exe_func(aniso_ma_test aniso_ma_test.cc) test_exe_func(torus_ma_test torus_ma_test.cc) test_exe_func(dg_ma_test dg_ma_test.cc) @@ -205,6 +206,7 @@ if(ENABLE_DSP) endif() if(ENABLE_SIMMETRIX) test_exe_func(curvetest curvetest.cc) + test_exe_func(aniso_adapt_sim aniso_adapt_sim.cc) test_exe_func(curve_to_bezier curve_to_bezier.cc) test_exe_func(make_all_cavities makeAllCavities.cc) test_exe_func(inClosureOf_test inClosureOf_test.cc) @@ -227,6 +229,7 @@ if(PUMI_ENABLE_CAPSTONE) util_exe_func(capGeomTest capGeomTest.cc) util_exe_func(capCheckParam capCheckParam.cc) util_exe_func(capAdapt capAdapt.cc) + test_exe_func(aniso_adapt_cap aniso_adapt_cap.cc) util_exe_func(capProbe capProbe.cc) util_exe_func(capLoadSome capLoadSome.cc) if(PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL) diff --git a/test/aniso_adapt.cc b/test/aniso_adapt.cc new file mode 100644 index 000000000..e604564f1 --- /dev/null +++ b/test/aniso_adapt.cc @@ -0,0 +1,24 @@ +#include +#include +#include +#include "aniso_adapt.h" + +int main(int argc, char** argv) +{ + PCU_ALWAYS_ASSERT(argc==3); + const char* modelFile = argv[1]; + const char* meshFile = argv[2]; + pcu::Init(&argc,&argv); + { + pcu::PCU PCUObj; + lion_set_verbosity(1); + gmi_register_mesh(); + + auto createMeshValues = [modelFile,meshFile,&PCUObj]() + { return apf::loadMdsMesh(modelFile,meshFile,&PCUObj); }; + + adaptTests(createMeshValues); + } + pcu::Finalize(); +} + diff --git a/test/aniso_adapt.h b/test/aniso_adapt.h new file mode 100644 index 000000000..9beaab9d5 --- /dev/null +++ b/test/aniso_adapt.h @@ -0,0 +1,231 @@ +#ifndef TEST_ANISO_ADAPT_H +#define TEST_ANISO_ADAPT_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +/* + Test some of the individual components in mesh adaptation to make sure that they are + functioning as intended. Right now it only tests coarsen refinement and snapping but + can be expanded to test more in the future. This has been tested with mds, simmetrix, + and capstone meshes. +*/ +class AnIso : public ma::AnisotropicFunction +{ + public: + AnIso(ma::Mesh* m, double sf1, double sf2) : + sizeFactor1(sf1), sizeFactor2(sf2) + { + average = ma::getAverageEdgeLength(m); + ma::getBoundingBox(m, lower, upper); + } + virtual void getValue(ma::Entity*, ma::Matrix& R, ma::Vector& H) + { + double h = average/sizeFactor1; + H = ma::Vector(h, h, h/sizeFactor2); + R = ma::Matrix( + 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0 + ); + } + private: + double sizeFactor1, sizeFactor2, average; + ma::Vector lower, upper; +}; + +void measureQuality(ma::Mesh* m, double& avgQuality, double& minQuality) +{ + int d = m->getDimension(); + apf::MeshEntity* elem; + apf::MeshIterator* it = m->begin(d); + ma::IdentitySizeField I(m); + while ((elem = m->iterate(it))) { + double q = ma::measureElementQuality(m, &I, elem); + avgQuality += q; + minQuality = std::min(minQuality, q); + } + m->end(it); + avgQuality = avgQuality / m->count(d); +} + +//make sure refinement is done +bool noLongEdges(ma::Adapt* a) +{ + apf::MeshEntity* edge; + apf::MeshIterator* it = a->mesh->begin(1); + while ((edge = a->mesh->iterate(it))) + if (a->sizeField->measure(edge) > ma::MAXLENGTH) + return false; + a->mesh->end(it); + return true; +} + +int countEdges(ma::Mesh* m) +{ + return m->count(1); +} + +ma::Mesh* coarsenForced(ma::Mesh* m) +{ + m->verify(); + AnIso sf(m, .5, 1); + ma::Input* in = ma::makeAdvanced(ma::configure(m, &sf)); + in->shouldForceAdaptation = true; + ma::Adapt* a = new ma::Adapt(in); + double avgQualBefore, avgQualAfter, minQualBefore, minQualAfter; + + measureQuality(m, avgQualBefore, minQualBefore); + double averageBefore = ma::getAverageEdgeLength(m); + int edgesBefore = countEdges(m); + + ma::coarsenMultiple(a); + + measureQuality(m, avgQualAfter, minQualAfter); + + //make sure that coarsening is happening and it doesn't make the mesh invalid + PCU_ALWAYS_ASSERT(edgesBefore > countEdges(m)); + PCU_ALWAYS_ASSERT(averageBefore < ma::getAverageEdgeLength(m)); + PCU_ALWAYS_ASSERT(minQualAfter >= in->validQuality); + + m->verify(); + delete a; + // cleanup input object and associated sizefield and solutiontransfer objects + if (in->ownsSizeField) + delete in->sizeField; + if (in->ownsSolutionTransfer) + delete in->solutionTransfer; + delete in; + return m; +} + +ma::Mesh* coarsenRegular(ma::Mesh* m) +{ + m->verify(); + AnIso sf(m, .5, 1); + ma::Input* in = ma::makeAdvanced(ma::configure(m, &sf)); + ma::Adapt* a = new ma::Adapt(in); + double avgQualBefore, avgQualAfter, minQualBefore, minQualAfter; + + measureQuality(m, avgQualBefore, minQualBefore); + double averageBefore = ma::getAverageEdgeLength(m); + int edgesBefore = countEdges(m); + + ma::coarsenMultiple(a); + + measureQuality(m, avgQualAfter, minQualAfter); + + //Makes sure that coarsening is happening and it isn't decreasing the quality of the mesh + PCU_ALWAYS_ASSERT(edgesBefore > countEdges(m)); + PCU_ALWAYS_ASSERT(averageBefore < ma::getAverageEdgeLength(m)); + PCU_ALWAYS_ASSERT(fabs(minQualBefore - minQualAfter) < 0.001f || minQualBefore < minQualAfter); + + m->verify(); + delete a; + // cleanup input object and associated sizefield and solutiontransfer objects + if (in->ownsSizeField) + delete in->sizeField; + if (in->ownsSolutionTransfer) + delete in->solutionTransfer; + delete in; + return m; +} + +//Makes sure that all of the vertices snapped to the model +bool allVertsOnModel(ma::Adapt* a) +{ + if (!a->input->shouldSnap) + return true; + ma::Mesh* m = a->mesh; + int dim = m->getDimension(); + ma::Iterator* it = m->begin(0); + ma::Entity* vert; + while ((vert = m->iterate(it))) { + int md = m->getModelType(m->toModel(vert)); + if (dim == 3 && md == 3) + continue; + ma::Vector modelPoint, param; + m->getPoint(vert,0,modelPoint); + m->getParam(vert,param); + ma::Model* model = m->toModel(vert); + m->snapToModel(model,param,modelPoint); + ma::Vector position = ma::getPosition(m, vert); + if (apf::areClose(modelPoint, position, 1e-12)) + continue; + } + m->end(it); + return true; +} + +ma::Mesh* refineSnapTest(ma::Mesh* m) +{ + m->verify(); + AnIso sf(m, 3, 1); + ma::Input* in = ma::makeAdvanced(ma::configure(m, &sf)); + ma::Adapt* a = new ma::Adapt(in); + int edgesBefore = countEdges(m); + double averageBefore = ma::getAverageEdgeLength(m); + + for (int i = 0; i < in->maximumIterations; ++i) + { + ma::refine(a); + ma::snap(a); + } + + double avgQualAfter, minQualAfter; + measureQuality(m, avgQualAfter, minQualAfter); + + //Makes sure that refinement is happening + PCU_ALWAYS_ASSERT(edgesBefore < countEdges(m)); + PCU_ALWAYS_ASSERT(averageBefore > ma::getAverageEdgeLength(m)); + + PCU_ALWAYS_ASSERT(noLongEdges(a)); + PCU_ALWAYS_ASSERT(allVertsOnModel(a)); + PCU_ALWAYS_ASSERT(minQualAfter >= 0); + + m->verify(); + delete a; + // cleanup input object and associated sizefield and solutiontransfer objects + if (in->ownsSizeField) + delete in->sizeField; + if (in->ownsSolutionTransfer) + delete in->solutionTransfer; + delete in; + return m; +} + +void adaptTests(const std::function& createMesh) +{ + //Mesh created multiple times to compare adaptation + ma::Mesh* meshReg = createMesh(); + apf::writeVtkFiles("startMesh", meshReg); + + refineSnapTest(meshReg); + apf::writeVtkFiles("afterRefine", meshReg); + + coarsenRegular(meshReg); + apf::writeVtkFiles("afterCoarsen", meshReg); + + ma::Mesh* meshForce = coarsenForced(refineSnapTest(createMesh())); + apf::writeVtkFiles("afterForcedCoarsen", meshForce); + + //Make sure setting to force coarsen is functioning + PCU_ALWAYS_ASSERT(countEdges(meshReg) > countEdges(meshForce)); + PCU_ALWAYS_ASSERT(ma::getAverageEdgeLength(meshReg) < ma::getAverageEdgeLength(meshForce)); + + meshReg->destroyNative(); + apf::destroyMesh(meshReg); +} + +#endif \ No newline at end of file diff --git a/test/aniso_adapt_cap.cc b/test/aniso_adapt_cap.cc new file mode 100644 index 000000000..9cc0bbec5 --- /dev/null +++ b/test/aniso_adapt_cap.cc @@ -0,0 +1,41 @@ +#include +#include +#include +#include "aniso_adapt.h" + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + { + pcu::PCU PCUObj = pcu::PCU(MPI_COMM_WORLD); + + if (argc < 2) { + if(PCUObj.Self()==0) + std::cerr << "usage: " << argv[0] + << " <.cre file> or\n" + << " <.smb file> \n"; + return EXIT_FAILURE; + } + + gmi_register_mesh(); + gmi_cap_start(); + gmi_register_cap(); + lion_set_verbosity(1); + const char* meshFile = argv[1]; + + gmi_model* model = gmi_load(meshFile); //Freed in adaptTests + ma::Mesh* apfCapMesh = apf::createCapMesh(model, &PCUObj); + apf::disownCapModel(apfCapMesh); + ma::Mesh* mesh = apf::createMdsMesh(model, apfCapMesh, true); + apf::disownMdsModel(mesh); + apf::destroyMesh(apfCapMesh); + + auto createMeshValues = [model, mesh]() + { return apf::createMdsMesh(model, mesh, true); }; + + adaptTests(createMeshValues); + apf::destroyMesh(mesh); + gmi_cap_stop(); + } + MPI_Finalize(); +} \ No newline at end of file diff --git a/test/aniso_adapt_sim.cc b/test/aniso_adapt_sim.cc new file mode 100644 index 000000000..39da0b102 --- /dev/null +++ b/test/aniso_adapt_sim.cc @@ -0,0 +1,57 @@ +#include "aniso_adapt.h" +#include +#include +#include +#include +#include "SimAdvMeshing.h" + +int main(int argc, char* argv[]) +{ + if (argc != 4) { + std::cerr << "Usage: " << argv[0] << + " " << std::endl; + return 1; + } + //Load Mesh + const char* nativefile = argv[1]; + const char* smdfile = argv[2]; + const char* smsfile = argv[3]; + MPI_Init(&argc, &argv); + pcu::PCU *PCUObj = new pcu::PCU(MPI_COMM_WORLD); + lion_set_verbosity(1); + gmi_register_mesh(); + + Sim_logOn("anisoadapt.log"); + MS_init(); + SimAdvMeshing_start(); + SimModel_start(); + Sim_readLicenseFile(NULL); + SimPartitionedMesh_start(&argc,&argv); + + gmi_sim_start(); + gmi_register_sim(); + + pProgress progress = Progress_new(); + Progress_setDefaultCallback(progress); + gmi_model* mdl_ref = gmi_sim_load(nativefile, smdfile); + pGModel model = gmi_export_sim(mdl_ref); + pParMesh mesh = PM_load(smsfile, model, progress); + ma::Mesh* mesh_ref = apf::createMesh(mesh, PCUObj); + + auto createMeshValues = [mdl_ref, mesh_ref]() + { return apf::createMdsMesh(mdl_ref, mesh_ref); }; + + adaptTests(createMeshValues); + + M_release(mesh); + Progress_delete(progress); + gmi_sim_stop(); + SimPartitionedMesh_stop(); + Sim_unregisterAllKeys(); + SimModel_stop(); + MS_exit(); + + delete PCUObj; + MPI_Finalize(); +} + diff --git a/test/aniso_ma_test.cc b/test/aniso_ma_test.cc index f1d99a360..e115d1c58 100644 --- a/test/aniso_ma_test.cc +++ b/test/aniso_ma_test.cc @@ -70,4 +70,3 @@ int main(int argc, char** argv) } pcu::Finalize(); } - diff --git a/test/capVol.cc b/test/capVol.cc index a6d0fc897..a7ccc0beb 100644 --- a/test/capVol.cc +++ b/test/capVol.cc @@ -241,7 +241,7 @@ apf::Splitter* makeSplitter(Args::Partitioner ptnr, apf::Mesh2* mesh) { default: #if defined(PUMI_HAS_ZOLTAN) return apf::makeZoltanSplitter( - adaptMesh, apf::GRAPH, apf::PARTITION + mesh, apf::GRAPH, apf::PARTITION ); #elif defined(PUMI_HAS_METIS) return apf::makeMETISsplitter(mesh); diff --git a/test/testing.cmake b/test/testing.cmake index 507420cbe..3ea560159 100644 --- a/test/testing.cmake +++ b/test/testing.cmake @@ -363,12 +363,21 @@ if(ENABLE_SIMMETRIX) "pipe_unif.smb" "pipe.smb") set_test_depends(TESTS snap_serial DEPENDS uniform_serial) + mpi_test(adapt_deformed_slab 1 + ./aniso_adapt_sim + "${MESHES}/deformedSlab/sim_model.x_t" + "${MESHES}/deformedSlab/sim_model.smd" + "${MESHES}/deformedSlab/sim_mesh.sms") endif() if(ENABLE_ZOLTAN) mpi_test(ma_serial 1 ./ma_test "${MDIR}/pipe.${GXT}" "pipe.smb") + mpi_test(adapt_cube 1 + ./aniso_adapt + "${MESHES}/cube/cube.dmg" + "${MESHES}/cube/pumi670/cube.smb") mpi_test(aniso_ma_serial 1 ./aniso_ma_test "${MESHES}/cube/cube.dmg" @@ -1016,6 +1025,7 @@ if(PUMI_ENABLE_CAPSTONE) endif() mpi_test(cap_inClosureOf 1 ./cap_inClosureOf ${MESHES}/cap/cube_surf_only.cre) mpi_test(cap_closestPoint 1 ./cap_closestPoint ${MESHES}/cap/cube_surf_only.cre) + mpi_test(adapt_curved_cube 1 ./aniso_adapt_cap "${MESHES}/curvedCube/cap_mesh.cre") if(PUMI_HAS_CAPSTONE_SIZINGMETRICTOOL) mpi_test(cap_smooth 1 ./cap_smooth) endif()