diff --git a/CMakeLists.txt b/CMakeLists.txt index 7a2ad4c00..88623f2d4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -172,6 +172,22 @@ if(ENABLE_CGNS) add_definitions(-DHAVE_CGNS) endif() +if(DEBUG_FPP) + add_definitions(-DDEBUG_FPP) +endif() + +if(ENABLE_EIGEN) + include(FetchContent) + FetchContent_Declare( + eigen + GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git + GIT_TAG 5.0 + ) + FetchContent_MakeAvailable(eigen) + include_directories(${eigen_SOURCE_DIR}) + add_definitions(-DEIGEN_ENABLED) +endif() + configure_file(SCOREC_config.h.in SCOREC_config.h) install(FILES "${CMAKE_BINARY_DIR}/SCOREC_config.h" DESTINATION include) include_directories(PUBLIC "$") diff --git a/apf/apf.cc b/apf/apf.cc index 662a69709..f589340c0 100644 --- a/apf/apf.cc +++ b/apf/apf.cc @@ -357,6 +357,13 @@ void getIntPoint(MeshElement* e, int order, int point, Vector3& param) param = p->param; } +Vector3 getIntPoint(MeshElement* e, int order, int point) +{ + IntegrationPoint const* p = + getIntegration(e->getType())->getAccurate(order)->getPoint(point); + return p->param; +} + double getIntWeight(MeshElement* e, int order, int point) { IntegrationPoint const* p = diff --git a/apf/apf.h b/apf/apf.h index 9a9ba849a..3d03379da 100644 --- a/apf/apf.h +++ b/apf/apf.h @@ -383,6 +383,15 @@ int countIntPoints(MeshElement* e, int order); */ void getIntPoint(MeshElement* e, int order, int point, Vector3& param); + +/** \brief Get an integration point in an element. + * + * \param order The polynomial order of accuracy. + * \param point The integration point number. + * \returns The resulting local coordinate of the integration point. + */ +Vector3 getIntPoint(MeshElement* e, int order, int point); + /** \brief Get the weight of an integration point in an element. * * \details All integration point tables in APF are scaled diff --git a/apf/apfCavityOp.cc b/apf/apfCavityOp.cc index 568fd2287..957460bff 100644 --- a/apf/apfCavityOp.cc +++ b/apf/apfCavityOp.cc @@ -119,6 +119,68 @@ void CavityOp::applyToDimension(int d) sharing = 0; } +void CavityOp::prepareListMigration(EntityList& elements, EntityList::iterator& iter, std::function migrationPossible) +{ + MeshTag* migTag = mesh->createIntTag("ma_migrate_list", 1); + std::vector iterators(elements.size()); + for (iter = elements.begin(); iter != elements.end();) { + int size = (int)iterators.size(); + mesh->setIntTag(*iter, migTag, &size); + iterators.push_back(iter++); + } + mesh->onDestroy = [&](MeshEntity* e) { + if (!mesh->hasTag(e, migTag)) return; + int i = 0; + mesh->getIntTag(e, migTag, &i); + iter = elements.erase(iterators[i]); + movedByDeletion = true; + }; + + migrationPossible(); + + mesh->onDestroy = 0; + for (iter = elements.begin(); iter != elements.end();) + mesh->removeTag(*iter++, migTag); + mesh->destroyTag(migTag); +} + +void CavityOp::applyToList(EntityList& elements) +{ + auto iter = elements.begin(); + prepareListMigration(elements, iter, [&]() { + do { + delete sharing; + sharing = apf::getSharing(mesh); + + isRequesting = false; + for (iter = elements.begin(); iter != elements.end();) { + if (sharing->isOwned(*iter)) { + if (setEntity(*iter) == OK) { + movedByDeletion=false; + apply(); + if (movedByDeletion) continue; + } + } + ++iter; + } + /* request any non-local cavities. + note: it is critical that this loop + be separated from the one above for + mesh-modifying operators, since an apply() + call could change other overlapping cavities */ + isRequesting = true; + for (iter = elements.begin(); iter != elements.end();) { + if (sharing->isOwned(*iter)) + setEntity(*iter); + ++iter; + } + } while (tryToPull()); + + delete sharing; + sharing = 0; + }); +} + bool CavityOp::requestLocality(MeshEntity** entities, int count) { bool areLocal = true; diff --git a/apf/apfCavityOp.h b/apf/apfCavityOp.h index 026b2875b..86c36ced2 100644 --- a/apf/apfCavityOp.h +++ b/apf/apfCavityOp.h @@ -14,6 +14,7 @@ #include "apfMesh.h" #include #include +#include namespace apf { @@ -87,6 +88,9 @@ class CavityOp virtual void apply() = 0; /** \brief parallel collective operation over entities of one dimension */ void applyToDimension(int d); + /** \brief parallel collective operation over entities in list */ + typedef std::list EntityList; + void applyToList(EntityList& elements); /** \brief within setEntity, require that entities be made local */ bool requestLocality(MeshEntity** entities, int count); /** \brief call before deleting a mesh entity during the operation */ @@ -102,6 +106,7 @@ class CavityOp bool tryToPull(); void applyLocallyWithModification(int d); void applyLocallyWithoutModification(int d); + void prepareListMigration(EntityList& elements, EntityList::iterator& iter, std::function migrationPossible); bool canModify; bool movedByDeletion; MeshIterator* iterator; diff --git a/apf/apfElement.cc b/apf/apfElement.cc index 534e6a35c..ac107cb56 100644 --- a/apf/apfElement.cc +++ b/apf/apfElement.cc @@ -26,6 +26,10 @@ void Element::init(Field* f, MeshEntity* e, VectorElement* p) getNodeData(); } +Element::Element() +{ +} + Element::Element(Field* f, MeshEntity* e) { init(f,e,0); diff --git a/apf/apfElement.h b/apf/apfElement.h index b54cd5901..e27e33f40 100644 --- a/apf/apfElement.h +++ b/apf/apfElement.h @@ -21,6 +21,7 @@ class VectorElement; class Element { public: + Element(); Element(Field* f, MeshEntity* e); Element(Field* f, VectorElement* p); virtual ~Element(); @@ -36,12 +37,12 @@ class Element FieldShape* getFieldShape() {return field->getShape();} void getComponents(Vector3 const& xi, double* c); void getElementNodeData(NewArray& d); - protected: void init(Field* f, MeshEntity* e, VectorElement* p); + protected: void getNodeData(); Field* field; Mesh* mesh; - MeshEntity* entity; + MeshEntity* entity=0; EntityShape* shape; VectorElement* parent; int nen; diff --git a/apf/apfElementOf.h b/apf/apfElementOf.h index afedf6111..5296d14a7 100644 --- a/apf/apfElementOf.h +++ b/apf/apfElementOf.h @@ -18,6 +18,10 @@ template class ElementOf : public Element { public: + ElementOf(): + Element() + { + } ElementOf(FieldOf* f, MeshEntity* e): Element(f,e) { diff --git a/apf/apfIntegrate.cc b/apf/apfIntegrate.cc index 2925524c1..d6b8858db 100644 --- a/apf/apfIntegrate.cc +++ b/apf/apfIntegrate.cc @@ -5,6 +5,7 @@ * BSD license as described in the LICENSE file in the top-level directory. */ +#include "apfVectorElement.h" #include "apfIntegrate.h" #include "apfMesh.h" #include "apf.h" @@ -651,9 +652,8 @@ void Integrator::process(Mesh* m, int d) while ((entity = m->iterate(elements))) { if ( ! m->isOwned(entity)) continue; - MeshElement* e = createMeshElement(m,entity); - this->process(e); - destroyMeshElement(e); + MeshElement e(m->getCoordinateField(),entity); + this->process(&e); } m->end(elements); this->parallelReduce(m->getPCU()); @@ -695,9 +695,8 @@ double measure(MeshElement* e) double measure(Mesh* m, MeshEntity* e) { - MeshElement* me = createMeshElement(m,e); - double v = measure(me); - destroyMeshElement(me); + MeshElement me(m->getCoordinateField(),e); + double v = measure(&me); return v; } diff --git a/apf/apfMatrix.cc b/apf/apfMatrix.cc index 92d15f402..6f6745fb9 100644 --- a/apf/apfMatrix.cc +++ b/apf/apfMatrix.cc @@ -9,6 +9,9 @@ #include "apf2mth.h" #include #include +#ifdef EIGEN_ENABLED +#include +#endif namespace apf { @@ -69,6 +72,22 @@ int eigen(Matrix3x3 const& A, Vector<3>* eigenVectors, double* eigenValues) { +#ifdef EIGEN_ENABLED + Eigen::Matrix3d matrix; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + matrix(i, j) = A[i][j]; + + Eigen::SelfAdjointEigenSolver solver(matrix); + Eigen::Vector3d values = solver.eigenvalues(); + Eigen::Matrix3d vectors = solver.eigenvectors(); + + for (unsigned i = 0; i < 3; ++i) + eigenValues[i] = values(i); + for (unsigned i = 0; i < 3; ++i) + for (unsigned j = 0; j < 3; ++j) + eigenVectors[j][i] = vectors(i,j); +#else mth::Matrix A2 = to_mth(A); mth::Matrix L; mth::Matrix Q; @@ -79,6 +98,7 @@ int eigen(Matrix3x3 const& A, for (unsigned i = 0; i < 3; ++i) for (unsigned j = 0; j < 3; ++j) eigenVectors[j][i] = Q(i,j); +#endif return 3; } diff --git a/apf/apfMesh.cc b/apf/apfMesh.cc index 80131a6e4..1bcd86d3f 100644 --- a/apf/apfMesh.cc +++ b/apf/apfMesh.cc @@ -174,6 +174,12 @@ Mesh::~Mesh() delete coordinateField; } +void Mesh::clearModelCache() +{ + periodicRanges.clear(); + boundingBoxes.clear(); +} + int Mesh::getModelType(ModelEntity* e) { return gmi_dim(getModel(), (gmi_ent*)e); @@ -228,6 +234,21 @@ bool Mesh::getPeriodicRange(ModelEntity* g, int axis, double range[2]) return gmi_periodic(getModel(), e, axis); } +bool Mesh::getPeriodicRangeCached(ModelEntity* g, int axis, double range[2]) +{ + auto it = periodicRanges.find({g, axis}); + if (it != periodicRanges.end()) { + range[0] = it->second.first[0]; + range[1] = it->second.first[1]; + return it->second.second; + } + else{ + bool output = getPeriodicRange(g, axis, range); + periodicRanges[{g, axis}] = {{range[0], range[1]}, output}; + return output; + } +} + void Mesh::getClosestPoint(ModelEntity* g, Vector3 const& from, Vector3& to, Vector3& p) { @@ -284,6 +305,20 @@ bool Mesh::isInClosureOf(ModelEntity* g, ModelEntity* target){ return (res == 1) ? true : false; } +void Mesh::boundingBoxCached(ModelEntity* g, + Vector3& bmin, Vector3& bmax) +{ + auto it = boundingBoxes.find(g); + if (it != boundingBoxes.end()) { + bmin = it->second.first; + bmax = it->second.second; + } + else{ + boundingBox(g, bmin, bmax); + boundingBoxes.insert(it, {g, {bmin, bmax}}); + } +} + void Mesh::boundingBox(ModelEntity* g, Vector3& bmin, Vector3& bmax) { diff --git a/apf/apfMesh.h b/apf/apfMesh.h index 02d25f7b5..1f21eded2 100644 --- a/apf/apfMesh.h +++ b/apf/apfMesh.h @@ -12,11 +12,14 @@ \brief The APF Mesh interface*/ #include +#include #include +#include #include #include #include "apfVector.h" #include "apfDynamicArray.h" +#include "gmi.h" struct gmi_model; @@ -288,6 +291,8 @@ class Mesh virtual gmi_model* getModel() = 0; /** \brief set the geometric model */ virtual void setModel(gmi_model* newModel) = 0; + /** \brief clear cached model info */ + void clearModelCache(); /** \brief return the model entity dimension */ int getModelType(ModelEntity* e); /** \brief get the dimension-unique model entity identifier */ @@ -309,6 +314,11 @@ class Mesh \returns true if (g) is periodic along this axis */ bool getPeriodicRange(ModelEntity* g, int axis, double range[2]); + /** \brief get the periodic properties of a model entity and chache the value + \param range if periodic, the parametric range + \returns true if (g) is periodic along this axis */ + bool getPeriodicRangeCached(ModelEntity* g, int axis, + double range[2]); /** \brief get closest point on geometry */ void getClosestPoint(ModelEntity* g, Vector3 const& from, Vector3& to, Vector3& p); @@ -325,6 +335,8 @@ class Mesh bool isInClosureOf(ModelEntity* g, ModelEntity* target); /** \brief get the bounding box of the model entity g */ void boundingBox(ModelEntity* g, Vector3& bmin, Vector3& bmax); + /** \brief get the bounding box of the model entity g and caches the value*/ + void boundingBoxCached(ModelEntity* g, Vector3& bmin, Vector3& bmax); /** \brief checks if p is on model g */ bool isOnModel(ModelEntity* g, Vector3 p, double scale); /** \brief get the distribution of the mesh's coordinate field */ @@ -401,12 +413,16 @@ class Mesh void switchPCU(pcu::PCU *newPCU); /** \brief true if any associated fields use array storage */ bool hasFrozenFields; + std::function onDestroy=0; protected: Field* coordinateField; std::vector fields; std::vector numberings; std::vector globalNumberings; pcu::PCU* pcu_; + std::map, + std::pair, bool>> periodicRanges; + std::map> boundingBoxes; }; /** \brief run consistency checks on an apf::Mesh structure diff --git a/apf/apfMesh2.h b/apf/apfMesh2.h index 2dda18e9b..01e2b27ef 100644 --- a/apf/apfMesh2.h +++ b/apf/apfMesh2.h @@ -102,6 +102,7 @@ class Mesh2 : public Mesh downward adjacencies */ void destroy(MeshEntity* e) { + if (onDestroy) onDestroy(e); requireUnfrozen(); destroy_(e); } diff --git a/apf/apfVectorElement.cc b/apf/apfVectorElement.cc index fa49e6a16..160af5d24 100644 --- a/apf/apfVectorElement.cc +++ b/apf/apfVectorElement.cc @@ -10,6 +10,16 @@ namespace apf { +VectorElement::VectorElement(): + ElementOf() +{ +} + +VectorElement::VectorElement(Field* f, MeshEntity* e): + ElementOf(static_cast(f),e) +{ +} + VectorElement::VectorElement(VectorField* f, MeshEntity* e): ElementOf(f,e) { diff --git a/apf/apfVectorElement.h b/apf/apfVectorElement.h index de4c3c447..4847e225c 100644 --- a/apf/apfVectorElement.h +++ b/apf/apfVectorElement.h @@ -18,6 +18,8 @@ class VectorField; class VectorElement : public ElementOf { public: + VectorElement(); + VectorElement(Field* f, MeshEntity* e); VectorElement(VectorField* f, MeshEntity* e); VectorElement(VectorField* f, VectorElement* p); virtual ~VectorElement() {} diff --git a/apf_cap/apfCAP.cc b/apf_cap/apfCAP.cc index 88ce841e2..0c27453c9 100644 --- a/apf_cap/apfCAP.cc +++ b/apf_cap/apfCAP.cc @@ -616,13 +616,6 @@ void MeshCAP::setModel(gmi_model* newModel) model = newModel; } -void MeshCAP::setModelEntity(MeshEntity* e, ModelEntity* me) -{ - (void)e; - (void)me; - apf::fail("MeshCAP::setModelEntity called!\n"); -} - static MeshMG::GeometryTopoType getCapGeomType(int d) { MeshMG::GeometryTopoType gtype = MeshMG::GVERTEX; @@ -645,6 +638,14 @@ static MeshMG::GeometryTopoType getCapGeomType(int d) return gtype; } +void MeshCAP::setModelEntity(MeshEntity* e, ModelEntity* me) +{ + M_MTopo topo = fromEntity(e); + M_GTopo gtopo = fromGmiEntity(reinterpret_cast(me)); + MeshMG::GeometryTopoType gtype = getCapGeomType(getModelType(me)); + meshInterface->set_geom_entity(topo, gtype, gtopo); +} + MeshEntity* MeshCAP::createVert_(ModelEntity* me) { double xyz[3] = {0., 0., 0.}; // this will be set later by setPoint_ diff --git a/ma/CMakeLists.txt b/ma/CMakeLists.txt index e54b59b38..b645c069f 100644 --- a/ma/CMakeLists.txt +++ b/ma/CMakeLists.txt @@ -28,6 +28,8 @@ set(SOURCES maSolutionTransferHelper.cc maSnap.cc maEdgeSwap.cc + maFaceSwap.cc + maFixShape.cc maShape.cc maShapeHandler.cc maQuality.cc @@ -59,6 +61,7 @@ set(HEADERS maMesh.h maSize.h maShape.h + maFixShape.h maTables.h maSolutionTransfer.h maExtrude.h diff --git a/ma/adapt_algorithms.tex b/ma/adapt_algorithms.tex new file mode 100644 index 000000000..2552d78a5 --- /dev/null +++ b/ma/adapt_algorithms.tex @@ -0,0 +1,374 @@ +\documentclass{article} +\usepackage{graphicx} % Required for inserting images +\usepackage{algorithm} +\usepackage{algpseudocode} +\usepackage{float} % For [H] placement +\usepackage[a4paper, margin=0.5in]{geometry} +\usepackage[font=scriptsize,labelfont=bf]{caption} +\captionsetup[figure]{justification=centering, font=small} +\usepackage{hyperref} +\hypersetup{ + colorlinks = true, + linkcolor = blue, + citecolor = red +} + +\title{Mesh Adapt Pseudo Code} +\author{} +\date{} + +\begin{document} + +\maketitle + +\section{Algorithms} + + +\begin{algorithm} +\caption{Adapt Top Level} +\begin{algorithmic}[1] + \State $Input:$ desired mesh size field + \State $MetricLength=$ loop over all edges to find the longest edge in the metric space + \State $NumIterations=$ $log_2(MetricLength)+1$ rounded up + \For{$i = 0$ to $NumIterations$} + \State coarsen() \Comment {remove short edges from input mesh and created during refinement and snapping} + \State refine() \Comment {split long edges to reach desired size} + \State snap() \Comment {move vertices created in refine step to model boundary} + \EndFor + \State fixShape() \Comment {remove or improve low quality tetrahedrons} +\end{algorithmic} +\end{algorithm} + + +\textbf{Coarsen:} Remove edges below minimum edge length in metric space using Algorithm 2. This is first to reduce the peak memory use by reducing the size of the mesh. It must be re-run because refinement and snapping can create additional short edges. + +\textbf{Refine:} Loop through all edges and flag edges above the maximum edge length in the metric space. Split all these edges and their adjacent elements once. This is done in steps because it results in higher quality elements. + +\textbf{Snap:} Move new vertices created from refinement to the model boundary using Algorithm 3. + +\textbf{FixShape:} loop through all the tetrahedrons to determine which are below the desired quality in the metric space and eliminate them using Algorithm 4. This is done at the end to clean up any poor quality tetrahedrons that could have been created by the user or previous steps. + +\begin{figure}[H] + \centering + \includegraphics[width=0.5\linewidth]{adapt_coarsen_order.png} + \caption{Importance of collapsing every other vertex. Reproduced from Li's Thesis} + \label{coarsenOrder} +\end{figure} + +\textbf{Coarsen:} The following algorithm collects a list of vertices adjacent to a short edge and then collapse those short edges while avoiding a collapse of two adjacent vertices in a row. \autoref{coarsenOrder} shows the issue with collapsing two adjacent vertices. Collapsing the doted vertices in example (a) is bad because it causes the edges to be pulled in one direction as shown in example (d). Instead collapsing one of the circled vertices in example (b) results in a higher quality mesh line in example (c). + +Lines 1-11, obtain a list of vertices adjacent to short edges and the flag is being used to avoid duplicates in the list. + +Lines 12-13, flag indicate whether more collapses need to be tried. This flag is changed in line 26 because a successful collapse creates more space for collapses to succeed and it also creates an dependent set that needs to be tried in future loops. + +Lines 14-18, obtain a vertex from the list and checking if its available for collapse using a flags. There are multiple reason a vertex may or may not be available for collapse. + +Line 18, the vertex is being flagged to prevent it from being checked again. + +Line 24, the flag is removed because a collapse can increase the changes of an adjacent collapse succeeding, so those vertices need to be processed again. + +Line 25, adjacent vertices are marked as a dependent set to prevent them from collapsing until the independent set is finished. + +Lines 19-23, attempt several collapses starting with the shortest edge first. The shortest edge first to prioritize getting rid of the shortest edges and they also result in the smallest change to the resulting mesh meaning that they are more likely to succeed and creates higher quality elements. A collapse might fail if it causes adjacent element to become invalid because of improper classification or inverted or negative volume. + +\begin{algorithm}[H] +\caption{Coarsen} +\begin{algorithmic}[1] + \ForAll {edges} + \If {edge length in metric space $<$ minimum edge length} + \ForAll {vertices $M^0_i$ in edge} \Comment {edges have 2 vertices} + \If {$M^0_i$ not flagged} \Comment {flag indicates if $M^0_i$ has been added to list} + \State add $M^0_i$ to the process list + \EndIf + \State flag $M^0_i$ as checked \Comment {flag to ensure $M^0_i$ added to list once} + \EndFor + \EndIf + \EndFor + \State remove the flag from all vertices + \While {$done == false$} \Comment{flag indicates whether more independent sets are necessary} + \State $done = true$ + \While {there there is unchecked vertex $M^0_i$} + \If {$M^0_i$ is $checked$ or is $dependent$} \Comment {flags indicates if $M^0_i$ can be collapsed} + \State \textbf{continue} + \EndIf + \State flag vertex $M^0_i$ as $checked$ \Comment {ensure $M^0_i$ not processed again unless cavity changes} + \State collect a list of adjacent edges $L$ + \State sort $L$ by edge length in metric space + \ForAll {edges $L_e$ in $L$ starting with the shortest} + \State evaluate a collapse on $L_e$ with $M^0_i$ removed + \If {collapse passes checks for classification, geometry and quality} + \State remove $checked$ flag from adjacent vertices + \State flag adjacent vertices as $dependent$ + \State $done = false$ + \State rebuild the cavity + \State \textbf{break} \Comment {stop evaluating and skip $M^0_{i+1}$} + \EndIf + \EndFor + \EndWhile + \State clear dependent flags from mesh + \EndWhile +\end{algorithmic} +\end{algorithm} + +\textbf{Refine:} The following algorithm collects a list of edges that need to be split into smaller edges and then splits those edges and create new adjacent elements and transfer all of the data to the new elements. + +Lines 1-5, find edges that are longer than the max edge length and flag them. The flag is used to avoid recomputing the expensive metric length computation. + +Lines 6-20, count the number of edges, faces, and regions that need to be split while using a flag to avoid double counting. + +Line 21, create a list with space for each entity counted in lines 6-20. This list is used to keep a pointer for cleanup of data after the old entities are split. + +Lines 22-38, populate the $Old$ list with the counted entities using the flags again to avoid adding the same entity twice. + +Line 39, create a list to store the new entities created. This list is the same size as the Old list except it contains a list of entities for each old entity. This list is used for transferring data to all of the new entities as a group after refinement is done. + +Lines 40-48, iterate through the old elements, split them, transfer some of the information and then add the new elements to the list. + +Lines 49-51, transfer the rest of the necessary data to the new elements. + +Lines 52-54, destroy the allocation for the old elements. + +Line 55, destroy the list pointing to the Old and New elements. + + +\begin{algorithm}[H] +\caption{Refine} +\begin{algorithmic}[1] +\ForAll {edges $E_i$ in mesh} + \If {$E_i$ length in metric space $>$ maximum edge length} + \State flag $E_i$ with split flag + \EndIf +\EndFor +\ForAll {edges $E_i$ in mesh} \Comment {flag and count upward elements} + \If {$E_i$ has split flag} + \State add $E_i$ to count + \ForAll {faces $F_i$ adjacent to $E_i$} + \If {$F_i$ does not have split flag} \Comment {used to prevent recounting} + \State flag $F_i$ with split flag and add it to count + \ForAll {regions $R_i$ adjacent to $F_i$} + \If {$R_i$ does not have split flag} + \State flag $R_i$ with split flag and add it to count \Comment {used to prevent recounting} + \EndIf + \EndFor + \EndIf + \EndFor + \EndIf +\EndFor +\State allocate 2D array $Old$ with space for each flagged edge, face, and region +\ForAll {edges $E_i$ in mesh} \Comment {add flagged edges, faces and regions to newly allocated array} + \If {$E_i$ has split flag} + \State $Old_{i} = E_i$ + \ForAll {faces $F_j$ adjacent to $E_i$} + \If {$F_j$ has split flag} \Comment {used to prevent re-adding} + \State $Old_{j} = F_j$ + \State clear split flag + \ForAll {regions $R_k$ adjacent to $F_j$} \Comment {used to prevent re-adding} + \If {$R_k$ has split flag} + \State $Old_{k} = R_k$ + \State clear split flag from $R_k$ + \EndIf + \EndFor + \EndIf + \EndFor + \EndIf +\EndFor +\State allocate 2D array $New$ with same size as $Old$ which contains edge, face, and region +\ForAll {elements $Old_i$ in $Old$} + \If {$Old_i$ is an edge} + \State create vertex + \State transfer parametric coordinates and interpolate + \State tag vertex to snap + \Else + \State find newly created downward elements $D$ in mesh + \State build new element connected to $D$ + \EndIf + \State add new elements to $New_i$ +\EndFor +\ForAll {elements $New_i$ in $New$} + \State transfer coordinate, solution transfer and size field data from $Old_i$ to $New_i$ +\EndFor +\ForAll {elements $Old_i$ in $Old$} + \State destroy $Old_i$ +\EndFor +\State de-allocate $Old$ and $New$ + +\end{algorithmic} +\end{algorithm} + + +\begin{figure}[H] + \centering + \includegraphics[width=0.7\linewidth]{adapt_tet_types.png} + \caption{Types of low quality tetrahedron \\ + Reproduced from Wan, J., "An automated adaptive procedure for 3D metal forming simulations", Ph.D. Thesis, Rensselaer Polytechnic Institute, Troy, NY 12180 (2006)} + \label{tetTypes} +\end{figure} + +\autoref{tetTypes} shows multiple types of tetrahedron that are not desirable. Section (a) includes two shapes that are undesirable due to having one or three short edges. Section b includes three shapes that are undesirable due to their low quality. Low quality shapes fall along a continuum of shapes between one, two, or three large dihedral angles. These shapes and how to remove or improve them shall be referenced in the following algorithms. Each operator checks that the operator should pass geometric, classification and quality checks. Each shape has multiple local operations that can be used to improve the quality, and different operations are considered in the order that is most likely to succeed. After a successful operation will stop execution and continue to the next tetrahedron. + +\textbf{Fix Shape:} The following algorithm takes a mesh with low quality tetrahedrons and either find the best method to remove them or to improve their quality. + +Lines 1-6, measures the quality of the tetrahedron and counting and flagging them. The count tracks if at improving the quality of the shapes. This is necessary because several operators (like split-collapse) can remove a low quality shape and create more than one slightly better low quality shape. So the count ensure that the algorithm doesn't get stuck creating more and more bad shapes. + +Lines 7-11, start with collapsing short edges because this operator is the most likely to succeed. + +Lines 8-16, attempt some collapses that result in the short edge merging into a larger edge. + +Lines 17-37, determine the type of low quality tetrahedrons and then evaluate operations in an order that optimizes for quality and runtime. After the first operation succeeds the rest don't need to be evaluated. Collapse operations first because they always decrease the number of elements. Then compound operators second because while they are more likely to succeed, there is a high chance that they create low quality elements. Reposition last because while it is almost guaranteed to improve the shape, it also has the slowest runtime. + + +\begin{figure}[H] + \centering + \includegraphics[width=1.0\linewidth]{adapt_local_operations.png} + \caption{Local mesh modification operations. Reproduced from SCOREC papers.} + \label{meshOperations} +\end{figure} + + +\begin{algorithm}[H] +\caption{Fix Element Shapes} +\begin{algorithmic}[1] +\ForAll {tetrahedrons} + \If {quality of tetrahedron $<$ minimum quality} + \State count and flag tetrahedron to be processed + \EndIf +\EndFor +\While {the number of low quality tetrahedron decreases} \Comment {first successful evaluate skips to next tetrahedron} + \ForAll {edges in tetrahedron} + \If {length of edge $<$ minimum length} \Comment{this value is cached} + \If {collapse short edge} + \State go to next tetrahedron + \EndIf + \EndIf + \EndFor + \ForAll {edges in tetrahedron} + \If {length of edge $M^1_i$ $<$ minimum length} \Comment{this value is cached} + \If {collapsing adjacent edges that result in $M^1_i$ becoming longer} + \State go to next tetrahedron + \EndIf + \EndIf + \EndFor + \If {\textbf{one large dihedral angle}} + \If {region-collapse if the tetrahedron is on the surface \\ + \hspace{.93cm} \textbf{or} edge-swap on the longest edge in the low quality triangle \\ + \hspace{.93cm} \textbf{or} split-collapse on the longest edge in the low quality triangle \\ + \hspace{.93cm} \textbf{or} edge-collapse on any edge in the low quality triangle besides the longest edge \\ + \hspace{.93cm} \textbf{or} reposition tetrahedron vertices to improve quality} + \State go to next tetrahedron + \EndIf + \EndIf + \If {\textbf{two large dihedral angles}} + \If {evaluate region-collapse if the tetrahedron is on the surface \\ + \hspace{.93cm} \textbf{or} edge-swap on $M^{1}_{j_0}$ or $M^{1}_{j_1}$ \\ + \hspace{.93cm} \textbf{or} split-collapse on $M^{1}_{j_0}$ or $M^{1}_{j_1}$ \\ + \hspace{.93cm} \textbf{or} double-split-collapse on $M^{1}_{j_0}$ and $M^{1}_{j_1}$ \\ + \hspace{.93cm} \textbf{or} reposition tetrahedron vertices to improve quality} + \State go to next tetrahedron + \EndIf + \ElsIf {\textbf{three large dihedral angles}} + \If {edge-collapse on $M^{1}_{j_0}$, $M^{1}_{j_1}$, or $M^{1}_{j_2}$ with $M^{0}_{i_3}$ deleted \\ + \hspace{.93cm} \textbf{or} edge-swap on any edge bounding triangle made by $M^{0}_{i_0}$, $M^{0}_{i_1}$, and $M^{0}_{i_2}$ \\ + \hspace{.93cm} \textbf{or} face-swap on the triangle made by $M^{0}_{i_0}$, $M^{0}_{i_1}$, $M^{0}_{i_2}$ \\ + \hspace{.93cm} \textbf{or} split-collapse on the triangle made by $M^{0}_{i_0}$, $M^{0}_{i_1}$, $M^{0}_{i_2}$ with $M^{0}_{i_3}$ deleted \\ + \hspace{.93cm} \textbf{or} face-split-collapse on the triangle made by $M^{0}_{i_0}$, $M^{0}_{i_1}$, $M^{0}_{i_2}$ with $M^{0}_{i_3}$ deleted \\ + \hspace{.93cm} \textbf{or} reposition vertex $M^{0}_{i_3}$ to improve quality} + \State go to next tetrahedron + \EndIf + \EndIf +\EndWhile + +\end{algorithmic} +\end{algorithm} + + +\begin{figure}[H] + \centering + \includegraphics[width=0.8\linewidth]{adapt_common_edges.png} + \caption{Types of common collapses. Reproduced from Li's Thesis} + \label{commonCollapses} +\end{figure} + +\textbf{Snapping:} Is a procedure where newly created vertices are relocated to the model boundary. Part of the algorithm is similar to fix shape, but with the direction of the operations restricted towards the model boundary. + +Lines 1-7, flag all of the vertices that need to be snapped. + +Lines 8-14, try to move the vertex directly to the model boundary. Most of the time this works and then it exits in line 13. Otherwise it reverts to the original state in line 11. + +Lines 15-24, collect the first problem plane. This encompasses the entities that need to be changed before progress can be made towards the target. As shown in \autoref{commonCollapses} it can have one or multiple tetrahedron. + +Lines 25-32, collect the common edges of the tetrahedron in the first problem plane. \autoref{commonCollapses}, shows the ideal way to reach the first problem plane is to collapse one of these edges. In example (a) any of the edges can be collapsed. In example (b) only two of the edges can be collapse with out creating invalid elements. In example (c) only one edge can be collapse. In example (d), none of the edges can be collapse. + +Line 33, attempt to collapse one of the common edges. + +Line 34, attempt to reduce one of the more complicated cases (b, c, d) into one of the simpler cases by collapsing shared edges of the first problem plane. + +Line 35, evaluate adjacent collapses that might reduce the complexity of the first problem plane and still make some progress towards the first problem plane. + +Lines 36-49, are a subset of fix shape, except only operations that make progress towards the first problem plane are considered. + + +\begin{algorithm}[H] +\caption{Snapping} +\begin{algorithmic}[1] + \State \textbf{Input:} mesh with flagged new vertices on model boundary + \While {there are any vertices flagged to snap} + \State move the current vertex $M^0_i$ to the model boundary \Comment {attempt simple movement first} + \If {any adjacent elements inverted} + \State move $M^0_i$ to original position + \Else + \State remove the flag to snap from $M^0_i$ then \textbf{continue} + \EndIf + \ForAll {invalid adjacent tetrahedron $M^3_i$} \Comment {find the first problem plane} + \State $ray = $ draw a ray from the start position to the target position + \State $M^1_i =$ the face opposite to $M^0_i$ in $M^3_i$ + \State $distance = $ distance to intersection of the ray and $M^1_i$ on the infinite plane + \If {distance $<$ minimum distance} + \State list FPP = $M^3_i$ + \ElsIf {distance == minimum distance} + \State add $M^3_i$ to the FPP list + \EndIf + \EndFor + \ForAll {edges $M^1_i$ in FPP} \Comment {find common edges} + \ForAll {tetrahedrons $M^3_j$ in FPP} + \State find if $M^1_i$ in $M^3_j$ + \EndFor + \If {$M^1_i$ in every tetrahedron} + \State add $M^1_i$ to the common edge list + \EndIf + \EndFor + \If {collapse common edges that move towards FPP with $M^0_i$ preserved \\ + \hspace{.4cm} \textbf{or} collapse edges on FPP that increase common edges \\ + \hspace{.4cm} \textbf{or} collapse common edges with $M^0_i$ preserved} + \State go to $M^0_{i+1}$ + \EndIf + \If {\textbf{one large dihedral angle}} + \If {swap the longest edge \\ + \hspace{.93cm} \textbf{or} split-collapse on the longest edge} + \State go to $M^0_{i+1}$ + \EndIf + \EndIf + \If {\textbf{two large dihedral angles}} + \If {swap either $M^1_{j0}$ or $M^1_{j1}$ \\ + \hspace{.93cm} \textbf{or} split-collapse on either $M^1_{j0}$ or $M^1_{j1}$ \\ + \hspace{.93cm} \textbf{or} double-split-collapse on both $M^1_{j0}$ and $M^1_{j1}$} + \State go to $M^0_{i+1}$ + \EndIf + \ElsIf {\textbf{three large dihedral angles}} + \If {swap any edge in face made by $M^{0}_{i_0}$, $M^{0}_{i_1}$, and $M^{0}_{i_2}$ \\ + \hspace{.93cm} \textbf{or} swap the problem face made by $M^{0}_{i_0}$, $M^{0}_{i_1}$, and $M^{0}_{i_2}$ \\ + \hspace{.93cm} \textbf{or} split-collapse on any edge in face made by $M^{0}_{i_0}$, $M^{0}_{i_1}$, and $M^{0}_{i_2}$} + \State go to $M^0_{i+1}$ + \EndIf + \EndIf + \EndWhile +\end{algorithmic} +\end{algorithm} + +\textbf{Future Work} + +\begin{itemize} + \item \textbf{Collapse-Swap:} This algorithm will swap invalid edges created a by a collapse in order to make the collapse succeed. This has the potential of improving runtime since not as many collapses have to be tested. In our current algorithm we have to test every short edge adjacent to a vertex, but in Li's algorithm he only has to test the shortest edge. In addition, collapsing the shortest edge generally results in better quality meshes since collapsing the shortest edge will result in the smallest change to the mesh. + \item \textbf{Vertex Reposition:} This addition to the algorithm will improve edge lengths. As described in Li's algorithm, when a collapse creates an edge that is longer than the desired maximum edge length, then vertex repositioning can be run in order to improve those edge lengths. +\end{itemize} + +\end{document} diff --git a/ma/adapt_coarsen_order.png b/ma/adapt_coarsen_order.png new file mode 100644 index 000000000..74a1ac2c6 Binary files /dev/null and b/ma/adapt_coarsen_order.png differ diff --git a/ma/adapt_common_edges.png b/ma/adapt_common_edges.png new file mode 100644 index 000000000..585d93ab0 Binary files /dev/null and b/ma/adapt_common_edges.png differ diff --git a/ma/adapt_local_operations.png b/ma/adapt_local_operations.png new file mode 100644 index 000000000..fe31c2757 Binary files /dev/null and b/ma/adapt_local_operations.png differ diff --git a/ma/adapt_tet_types.png b/ma/adapt_tet_types.png new file mode 100644 index 000000000..4856009c6 Binary files /dev/null and b/ma/adapt_tet_types.png differ diff --git a/ma/ma.cc b/ma/ma.cc index b1aa8c3f9..f7c769968 100644 --- a/ma/ma.cc +++ b/ma/ma.cc @@ -16,7 +16,10 @@ #include "maBalance.h" #include "maLayer.h" #include "maDBG.h" +#include "maFixShape.h" #include +#include +#include "apfGeometry.h" namespace ma { @@ -43,6 +46,11 @@ void adapt(Input* in) printQuality(a); postBalance(a); Mesh* m = a->mesh; + + double t1 = pcu::Time(); + print(m->getPCU(), "mesh adapted in %f seconds", t1-t0); + apf::printStats(m); + delete a; // cleanup input object and associated sizefield and solutiontransfer objects if (in->ownsSizeField) @@ -50,9 +58,6 @@ void adapt(Input* in) if (in->ownsSolutionTransfer) delete in->solutionTransfer; delete in; - double t1 = pcu::Time(); - print(m->getPCU(), "mesh adapted in %f seconds", t1-t0); - apf::printStats(m); } void adapt(const Input* in) diff --git a/ma/maAdapt.cc b/ma/maAdapt.cc index bff23cf90..4da716cce 100644 --- a/ma/maAdapt.cc +++ b/ma/maAdapt.cc @@ -24,12 +24,22 @@ namespace ma { +//used to avoid expensive cube root operations +void setPerformanceQuality(Input* in) +{ + if (in->mesh->getDimension()==3) + in->goodQuality = std::pow(in->goodQuality, 3); + else + in->goodQuality = std::pow(in->goodQuality, 2); +} + Adapt::Adapt(Input* in) { + setPerformanceQuality(in); input = in; mesh = in->mesh; setupFlags(this); - setupQualityCache(this); + setupCache(this); deleteCallback = 0; buildCallback = 0; sizeField = in->sizeField; @@ -52,7 +62,7 @@ Adapt::Adapt(Input* in) Adapt::~Adapt() { clearFlags(this); - clearQualityCache(this); + clearCache(this); delete refine; delete shape; } @@ -60,6 +70,7 @@ Adapt::~Adapt() void setupFlags(Adapt* a) { a->flagsTag = a->mesh->createIntTag("ma_flags",1); + a->snapTag = a->mesh->createDoubleTag("ma_snap",3); } void clearFlags(Adapt* a) @@ -69,12 +80,16 @@ void clearFlags(Adapt* a) for (int d=0; d <= 3; ++d) { Iterator* it = m->begin(d); - while ((e = m->iterate(it))) + while ((e = m->iterate(it))) { if (m->hasTag(e,a->flagsTag)) m->removeTag(e,a->flagsTag); + if (m->hasTag(e,a->snapTag)) + m->removeTag(e,a->snapTag); + } m->end(it); } m->destroyTag(a->flagsTag); + m->destroyTag(a->snapTag); } int getFlags(Adapt* a, Entity* e) @@ -147,25 +162,50 @@ void clearFlagFromDimension(Adapt* a, int flag, int dimension) m->end(it); } -void setupQualityCache(Adapt* a) +void setupCache(Adapt* a) { a->qualityCache = a->mesh->createDoubleTag("ma_qual_cache",1); + a->sizeCache = a->mesh->createDoubleTag("ma_size_cache",1); } -void clearQualityCache(Adapt* a) +void clearCache(Adapt* a) { Mesh* m = a->mesh; Entity* e; - // only faces and regions can have the quality tag - for (int d=2; d <= 3; ++d) + for (int d=1; d <= 3; ++d) { Iterator* it = m->begin(d); - while ((e = m->iterate(it))) + while ((e = m->iterate(it))) { if (m->hasTag(e,a->qualityCache)) m->removeTag(e,a->qualityCache); + if (m->hasTag(e,a->sizeCache)) + m->removeTag(e,a->sizeCache); + } m->end(it); } m->destroyTag(a->qualityCache); + m->destroyTag(a->sizeCache); + m->clearModelCache(); +} + +double getCachedSize(Adapt* a, Entity* e) +{ + int type = a->mesh->getType(e); + int ed = apf::Mesh::typeDimension[type]; + PCU_ALWAYS_ASSERT(ed != 0); + if (!a->mesh->hasTag(e,a->sizeCache)) + return 0.0; //we assume 0.0 is the default value for all qualities + double size; + a->mesh->getDoubleTag(e,a->sizeCache,&size); + return size; +} + +void setCachedSize(Adapt* a, Entity* e, double size) +{ + int type = a->mesh->getType(e); + int ed = apf::Mesh::typeDimension[type]; + PCU_ALWAYS_ASSERT(ed != 0); + a->mesh->setDoubleTag(e,a->sizeCache,&size); } double getCachedQuality(Adapt* a, Entity* e) diff --git a/ma/maAdapt.h b/ma/maAdapt.h index 0b7a0324d..cf058e79d 100644 --- a/ma/maAdapt.h +++ b/ma/maAdapt.h @@ -50,6 +50,8 @@ class Adapt Mesh* mesh; Tag* flagsTag; Tag* qualityCache; // to avoid repeated quality computations + Tag* sizeCache; // to avoid repeated size computations + Tag* snapTag; DeleteCallback* deleteCallback; apf::BuildCallback* buildCallback; SizeField* sizeField; @@ -76,8 +78,10 @@ void clearFlagMatched(Adapt* a, Entity* e, int flag); void clearFlagFromDimension(Adapt* a, int flag, int dimension); -void setupQualityCache(Adapt* a); -void clearQualityCache(Adapt* a); +void setupCache(Adapt* a); +void clearCache(Adapt* a); +double getCachedSize(Adapt* a, Entity* e); +void setCachedSize(Adapt* a, Entity* e, double size); double getCachedQuality(Adapt* a, Entity* e); void setCachedQuality(Adapt* a, Entity* e, double q); diff --git a/ma/maCoarsen.cc b/ma/maCoarsen.cc index 9baa67130..39140931f 100644 --- a/ma/maCoarsen.cc +++ b/ma/maCoarsen.cc @@ -8,13 +8,8 @@ *******************************************************************************/ /* -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. + coarsen will repeatedly create independent sets and collapse the mesh + until there are no more short edges left. */ #include "maCoarsen.h" #include "maAdapt.h" @@ -27,60 +22,11 @@ This file contains two coarsening alogrithms. #include #include #include +#include "maSnapper.h" +#include "maShapeHandler.h" 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; -} -} +const double MINIMUM_QUALITY = std::pow(.1, 3); class CollapseChecker : public apf::CavityOp { @@ -223,247 +169,146 @@ int collapseAllEdges(Adapt* a, int modelDimension) return collapser.successCount; } -class MatchedEdgeCollapser : public Operator -{ - public: - MatchedEdgeCollapser(Adapt* a, int md): - modelDimension(md), - collapse(a) - { - successCount = 0; - } - virtual int getTargetDimension() {return 1;} - virtual bool shouldApply(Entity* e) - { - Adapt* a = getAdapt(); - if ( ! getFlag(a,e,COLLAPSE)) - return false; - Mesh* m = a->mesh; - int md = m->getModelType(m->toModel(e)); - if (md!=modelDimension) - return false; - collapse.setEdge(e); - return true; - } - virtual bool requestLocality(apf::CavityOp* o) - { - return collapse.requestLocality(o); - } - virtual void apply() - { - double qualityToBeat = getAdapt()->input->validQuality; - collapse.setEdges(); - if ( ! collapse.checkTopo()) - return; - if ( ! collapse.tryBothDirections(qualityToBeat)) - return; - collapse.destroyOldElements(); - ++successCount; - } - Adapt* getAdapt() {return collapse.adapt;} - int successCount; - private: - int modelDimension; - MatchedCollapse collapse; -}; -static int collapseMatchedEdges(Adapt* a, int modelDimension) +//returns a list of vertices that bound a short edge +std::list getShortEdgeVerts(Adapt* a) { - MatchedEdgeCollapser collapser(a, modelDimension); - applyOperator(a, &collapser); - return collapser.successCount; + std::list shortEdgeVerts; + Iterator* it = a->mesh->begin(1); + Entity* edge; + while ((edge = a->mesh->iterate(it))) + { + if (getAndCacheSize(a, edge) > 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; + if (!a->mesh->isOwned(vertices[i])) continue; + if (a->mesh->isShared(vertices[i])) 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; } -struct ShouldCollapse : public Predicate +class CollapseAll : public apf::CavityOp { - ShouldCollapse(Adapt* a_):a(a_) {} - bool operator()(Entity* e) - { - return a->sizeField->shouldCollapse(e); - } + public: Adapt* a; -}; + Entity* vertex; + Collapse collapse; + std::list shortEdgeVerts; -long markEdgesToCollapse(Adapt* a) -{ - ShouldCollapse p(a); - return markEntities(a, 1, p, COLLAPSE, NEED_NOT_COLLAPSE, - DONT_COLLAPSE | NEED_NOT_COLLAPSE); -} + int success=0; -bool coarsen(Adapt* a) -{ - if (!a->input->shouldCoarsen) - return false; - double t0 = pcu::Time(); - --(a->coarsensLeft); - long count = markEdgesToCollapse(a); - if ( ! count) - return false; - Mesh* m = a->mesh; - int maxDimension = m->getDimension(); - PCU_ALWAYS_ASSERT(checkFlagConsistency(a,1,COLLAPSE)); - long successCount = 0; - for (int modelDimension=1; modelDimension <= maxDimension; ++modelDimension) + CollapseAll(Adapt* adapt) : apf::CavityOp(adapt->mesh,true) { - checkAllEdgeCollapses(a,modelDimension); - findIndependentSet(a); - if (m->hasMatching()) - successCount += collapseMatchedEdges(a, modelDimension); - else - successCount += collapseAllEdges(a, modelDimension); + a = adapt; + collapse.Init(a); + shortEdgeVerts = getShortEdgeVerts(a); } - successCount = m->getPCU()->Add(successCount); - double t1 = pcu::Time(); - print(m->getPCU(), "coarsened %li edges in %f seconds", successCount,t1-t0); - return true; -} -/* + ~CollapseAll() + { + ma::clearFlagFromDimension(a, CHECKED, 0); + } + + int getTargetDimension() { return 0; } + void resetIndependentSet() { ma::clearFlagFromDimension(a, NEED_NOT_COLLAPSE, 0); } + + /* 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--; + */ + void flagIndependentSet(apf::Up& adjacent) + { + 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 } } } -} -struct EdgeLength -{ - Entity* edge; - double length; - bool operator<(const EdgeLength& other) const { - return length < other.length; - } -}; + 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); + */ + bool collapseShortest(Entity* vertex, apf::Up& adjacent) + { + double qualityToBeat = a->input->shouldForceAdaptation ? MINIMUM_QUALITY + : a->input->goodQuality; + + std::vector sorted; + for (int i=0; i < adjacent.n; i++) { + double length = getAndCacheSize(a, adjacent.e[i]); + EdgeLength measured{adjacent.e[i], length}; + if (measured.length > MINLENGTH) continue; + sorted.push_back(measured); + } + std::sort(sorted.begin(), sorted.end()); + for (size_t i=0; i < sorted.size(); i++) { + if (!collapse.run(sorted[i].edge, vertex, qualityToBeat)) continue; + flagIndependentSet(adjacent); + collapse.destroyOldElements(); + return true; + } 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; -} + Outcome setEntity(Entity* vert) + { + vertex = vert; + if (getFlag(a, vertex, CHECKED)) return SKIP; + if (getFlag(a, vertex, NEED_NOT_COLLAPSE)) return SKIP; + if (!requestLocality(&vertex, 1)) + return REQUEST; + return OK; + } -//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))) + void apply() { - 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]); - } + apf::Up adjacent; + a->mesh->getUp(vertex, adjacent); + setFlag(a, vertex, CHECKED); + if (collapseShortest(vertex, adjacent)) success++; } - 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) +bool coarsen(Adapt* a) { - if (!a->input->shouldCoarsen) - return false; + 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); + CollapseAll collapseAll(a); + int totalSuccess = 0; 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); + do { + collapseAll.success = 0; + collapseAll.applyToList(collapseAll.shortEdgeVerts); + collapseAll.resetIndependentSet(); + success = a->mesh->getPCU()->Add(collapseAll.success); + totalSuccess += success; + } while (success > 0); + double t1 = pcu::Time(); - print(a->mesh->getPCU(), "coarsened %d edges in %f seconds", success, t1-t0); + print(a->mesh->getPCU(), "coarsened %d edges in %f seconds", totalSuccess, t1-t0); return true; } diff --git a/ma/maCoarsen.h b/ma/maCoarsen.h index 2e4283c77..bf83b1e54 100644 --- a/ma/maCoarsen.h +++ b/ma/maCoarsen.h @@ -14,7 +14,6 @@ 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 af13ba324..14cdae83e 100644 --- a/ma/maCollapse.cc +++ b/ma/maCollapse.cc @@ -10,8 +10,11 @@ #include "maCollapse.h" #include "maAdapt.h" #include "maShape.h" +#include "maShapeHandler.h" #include #include +#include "apfGeometry.h" +#include "maDBG.h" namespace ma { @@ -24,6 +27,64 @@ void Collapse::Init(Adapt* a) vertToKeep = 0; } +bool Collapse::run(Entity* edge, Entity* vert, double qualityToBeat) +{ + PCU_ALWAYS_ASSERT(adapt->mesh->getType(edge) == apf::Mesh::EDGE); + PCU_ALWAYS_ASSERT(adapt->mesh->getType(vert) == apf::Mesh::VERTEX); + if (!setEdge(edge)) return false; + if (!checkClass()) return false; + if (!getFlag(adapt, vert, COLLAPSE)) { unmark(); return false; } + vertToCollapse = vert; + vertToKeep = getEdgeVertOppositeVert(adapt->mesh, edge, vert); + clearFlag(adapt, vertToKeep, COLLAPSE); + computeElementSets(); + if (elementsToKeep.size() == 0) { unmark(); return false; } + if (!isValid()) { unmark(); return false; } + + if (!checkEdgeCollapseTopology(adapt, edge)) { unmark(); return false; } + if (!getFlag(adapt, vert, COLLAPSE)) { unmark(); return false; } + PCU_ALWAYS_ASSERT(vertToCollapse == vert); + + if (!adapt->input->shouldForceAdaptation) + qualityToBeat = std::min(adapt->input->goodQuality, + std::max(getOldQuality(),adapt->input->validQuality)); + + if (anyWorseQuality(qualityToBeat)) { unmark(); return false;} + + rebuildElements(); + fitElements(); + unmark(); + + if (adapt->mesh->getDimension()==2 && !isGood2DMesh()) {destroyNewElements(); return false;} + return true; +} + +bool Collapse::run(Entity* edge, double qualityToBeat) +{ + PCU_ALWAYS_ASSERT(adapt->mesh->getType(edge) == apf::Mesh::EDGE); + if (!setEdge(edge)) return false; + if (!checkClass()) return false; + if (!checkTopo()) return false; + + computeElementSets(); + if (!adapt->input->shouldForceAdaptation) + qualityToBeat = std::min(adapt->input->goodQuality, + std::max(getOldQuality(), adapt->input->validQuality)); + + if (!isValid() || anyWorseQuality(qualityToBeat)) { + if (!getFlag(adapt, vertToKeep, COLLAPSE)) { unmark(); return false; } + std::swap(vertToKeep, vertToCollapse); + computeElementSets(); + if (!isValid() || anyWorseQuality(qualityToBeat)) { unmark(); return false; } + } + + rebuildElements(); + fitElements(); + unmark(); + return true; +} + + bool Collapse::requestLocality(apf::CavityOp* o) { /* get vertices again since this is sometimes used @@ -34,6 +95,36 @@ bool Collapse::requestLocality(apf::CavityOp* o) return o->requestLocality(v,2); } +bool Collapse::isValid() +{ + PCU_ALWAYS_ASSERT(!adapt->mesh->isShared(vertToCollapse)); + if (adapt->mesh->getDimension() == 3 && !cavity.shouldFit) return true; + Vector prev = getPosition(adapt->mesh, vertToCollapse); + Vector target = getPosition(adapt->mesh, vertToKeep); + adapt->mesh->setPoint(vertToCollapse, 0, target); + bool valid = true; + for (Entity* e : elementsToKeep) + if (!isTetValid(adapt->mesh, e)) {valid = false; break;} + adapt->mesh->setPoint(vertToCollapse, 0, prev); + return valid; +} + +bool Collapse::anyWorseQuality(double qualityToBeat) +{ + Vector prev = getPosition(adapt->mesh, vertToCollapse); + Vector target = getPosition(adapt->mesh, vertToKeep); + auto prevMatrix = adapt->sizeField->getValue(vertToCollapse); + auto targetMatrix = adapt->sizeField->getValue(vertToKeep); + adapt->mesh->setPoint(vertToCollapse, 0, target); + adapt->sizeField->setValue(vertToCollapse, targetMatrix); + bool worse = false; + for (Entity* e : elementsToKeep) + if (adapt->shape->getQuality(e) < qualityToBeat) {worse = true; break;} + adapt->mesh->setPoint(vertToCollapse, 0, prev); + adapt->sizeField->setValue(vertToCollapse, prevMatrix); + return worse; +} + bool Collapse::tryThisDirectionNoCancel(double qualityToBeat) { PCU_ALWAYS_ASSERT( ! adapt->mesh->isShared(vertToCollapse)); @@ -158,9 +249,9 @@ bool checkEdgeCollapseEdgeRings(Adapt* a, Entity* edge) Mesh* m = a->mesh; Entity* v[2]; m->getDownward(edge,0,v); - if (!getFlag(a, v[0], DONT_COLLAPSE)) //Allow collapse in one direction + if (getFlag(a, v[0], 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 + if (getFlag(a, v[1], COLLAPSE)) //Allow collapse in one direction PCU_ALWAYS_ASSERT( ! m->isShared(v[1])); apf::Up ve[2]; m->getUp(v[0],ve[0]); @@ -363,10 +454,50 @@ void Collapse::computeElementSets() APF_ITERATE(Upward,adjacent,it) if ( ! elementsToCollapse.count(*it)) elementsToKeep.insert(*it); - PCU_ALWAYS_ASSERT(elementsToKeep.size()); } -void Collapse::rebuildElements() +//Find edges and faces that can be reused in new entities after collapse +std::map Collapse::getReusableEntities() +{ + Mesh* m = adapt->mesh; + std::map reusable; + for (Entity* elm : elementsToCollapse) { + Entity* faces[4]; + m->getDownward(elm, 2, faces); + Entity* faceToKeep = 0; + Entity* faceToReplace = 0; + for (int f=0; f<4; f++) { + Entity* edges[3]; + m->getDownward(faces[f], 1, edges); + Entity* edgeToDelete=0; + Entity* edgeToKeep=0; + Entity* edgeToReplace=0; + for (int e=0; e<3; e++) { + if (edges[e] == edge) edgeToDelete = edges[e]; + else if (isInClosure(m, edges[e], vertToKeep)) edgeToKeep = edges[e]; + else if (isInClosure(m, edges[e], vertToCollapse)) edgeToReplace = edges[e]; + } + if (edgeToReplace != 0 && edgeToDelete == 0 && edgeToKeep == 0) + faceToReplace = faces[f]; + if (edgeToKeep != 0 && edgeToDelete == 0 && edgeToReplace == 0) + faceToKeep = faces[f]; + if (edgeToDelete != 0) + reusable[edgeToReplace] = edgeToKeep; + } + reusable[faceToReplace] = faceToKeep; + } + return reusable; +} + +Entity* Collapse::rebuildEntity(Mesh* m, Entity* original, Entity** downward) +{ + Entity* entity = m->createEntity(m->getType(original), m->toModel(original), downward); + if (adapt->buildCallback) adapt->buildCallback->call(entity); + if (rebuildCallback) rebuildCallback->rebuilt(entity, original); + return entity; +} + +void Collapse::rebuildElements2D() { PCU_ALWAYS_ASSERT(elementsToKeep.size()); newElements.setSize(elementsToKeep.size()); @@ -379,6 +510,47 @@ void Collapse::rebuildElements() cavity.afterBuilding(); } +void Collapse::rebuildElements() +{ + if (adapt->mesh->getDimension() < 3) return rebuildElements2D(); + PCU_ALWAYS_ASSERT(elementsToKeep.size()); + newElements.setSize(elementsToKeep.size()); + cavity.beforeBuilding(); + size_t ni=0; + + Mesh* m = adapt->mesh; + std::map rebuilt = getReusableEntities(); + APF_ITERATE(EntitySet, elementsToKeep, it) { + Entity* tetFaces[4]; + m->getDownward(*it, 2, tetFaces); + for (int f=0; f<4; f++) { + auto foundFace = rebuilt.find(tetFaces[f]); + if (foundFace != rebuilt.end()) {tetFaces[f] = foundFace->second; continue;} + if (!isInClosure(m, tetFaces[f], vertToCollapse)) continue; + + Entity* faceEdges[3]; + m->getDownward(tetFaces[f], 1, faceEdges); + for (int e=0; e<3; e++) { + auto foundEdge = rebuilt.find(faceEdges[e]); + if (foundEdge != rebuilt.end()) {faceEdges[e] = foundEdge->second; continue;} + Entity* edgeVerts[2]; + m->getDownward(faceEdges[e], 0, edgeVerts); + if (edgeVerts[0] == vertToCollapse) edgeVerts[0] = vertToKeep; + else if (edgeVerts[1] == vertToCollapse) edgeVerts[1] = vertToKeep; + else continue; + Entity* newEdge = rebuildEntity(m, faceEdges[e], edgeVerts); + rebuilt[faceEdges[e]] = newEdge; + faceEdges[e] = newEdge; + } + Entity* newFace = rebuildEntity(m, tetFaces[f], faceEdges); + rebuilt[tetFaces[f]] = newFace; + tetFaces[f] = newFace; + } + newElements[ni++] = rebuildEntity(m, *it, tetFaces); + } + cavity.afterBuilding(); +} + void Collapse::fitElements() { if(cavity.shouldFit){ diff --git a/ma/maCollapse.h b/ma/maCollapse.h index d9a449adf..c717e6ef6 100644 --- a/ma/maCollapse.h +++ b/ma/maCollapse.h @@ -24,6 +24,8 @@ class Collapse { public: void Init(Adapt* a); + bool run(Entity* edge, Entity* vert, double qualityToBeat); + bool run(Entity* edge, double qualityToBeat); bool requestLocality(apf::CavityOp* o); void destroyOldElements(); void destroyNewElements(); @@ -34,6 +36,7 @@ class Collapse void setVerts(); virtual void computeElementSets(); void rebuildElements(); + void rebuildElements2D(); void fitElements(); bool isGood2DMesh(); void cancel(); @@ -42,6 +45,8 @@ class Collapse bool tryBothDirections(double qualityToBeat); void getOldElements(EntityArray& oldElements); bool edgesGoodSize(); + bool isValid(); + bool anyWorseQuality(double qualityToBeat); double getOldQuality(); Adapt* adapt; Entity* edge; @@ -52,11 +57,15 @@ class Collapse EntityArray newElements; Cavity cavity; RebuildCallback* rebuildCallback; + private: + std::map getReusableEntities(); + Entity* rebuildEntity(Mesh* m, Entity* original, Entity** downward); }; bool checkEdgeCollapseTopology(Adapt* a, Entity* edge); bool isRequiredForAnEdgeCollapse(Adapt* adapt, Entity* vertex); bool isRequiredForMatchedEdgeCollapse(Adapt* adapt, Entity* vertex); +bool collapseEdge(Collapse& collapse, Entity* edge, double qualityToBeat); bool setupCollapse(Collapse& collapse, Entity* edge, Entity* vert); } diff --git a/ma/maDBG.cc b/ma/maDBG.cc index ed78f92df..8d8e01a9d 100644 --- a/ma/maDBG.cc +++ b/ma/maDBG.cc @@ -46,28 +46,67 @@ void writeMesh(ma::Mesh* m, apf::writeVtkFiles(fileName, m, dim); } -void addClassification(ma::Adapt* a, - const char* fieldName) +apf::Field* getField(ma::Adapt* a, int dim, const char* fieldName) { - ma::Mesh* m = a->mesh; - apf::Field* field; - field = m->findField(fieldName); + apf::Field* field = a->mesh->findField(fieldName); if (field) apf::destroyField(field); + if (dim == 0) + field = apf::createFieldOn(a->mesh, fieldName, apf::SCALAR); + else + field = apf::createField(a->mesh, fieldName, apf::SCALAR, apf::getConstant(dim)); + return 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); +void useFieldInfo(ma::Adapt* a, const std::function& funcUsingField) +{ + ma::Mesh* m = a->mesh; + apf::Field* fieldVert = getField(a, 0, "vert_classification"); + apf::Field* fieldEdge = getField(a, 1, "edge_classification"); + apf::Field* fieldFace = getField(a, 2, "face_classification"); + + ma::Entity* e; + ma::Iterator* it = m->begin(0); + while ((e = m->iterate(it))){ + double modelDimension = (double)m->getModelType(m->toModel(e)); + apf::setComponents(fieldVert, e, 0, &modelDimension); + } + m->end(it); + + it = m->begin(1); + while ((e = m->iterate(it))){ + double modelDimension = (double)m->getModelType(m->toModel(e)); + apf::setComponents(fieldEdge, e, 0, &modelDimension); } m->end(it); + + it = m->begin(2); + while ((e = m->iterate(it))){ + double modelDimension = (double)m->getModelType(m->toModel(e)); + apf::setComponents(fieldFace, e, 0, &modelDimension); + } + m->end(it); + + // measure qualities + std::vector lq_metric; + std::vector lq_no_metric; + ma::getLinearQualitiesInMetricSpace(a->mesh, a->sizeField, lq_metric); + ma::getLinearQualitiesInPhysicalSpace(a->mesh, lq_no_metric); + apf::Field* metricField = colorEntitiesOfDimWithValues(a, a->mesh->getDimension(), lq_metric, "qual_metric"); + apf::Field* physicalField = colorEntitiesOfDimWithValues(a, a->mesh->getDimension(), lq_no_metric, "qual_no_metric"); + apf::Field* snapTargetField = 0; + if (a->input->shouldSnap) snapTargetField = addTargetLocation(a, "snap_target"); + + funcUsingField(); + apf::destroyField(fieldVert); + apf::destroyField(fieldEdge); + apf::destroyField(fieldFace); + apf::destroyField(metricField); + apf::destroyField(physicalField); + if (a->input->shouldSnap) apf::destroyField(snapTargetField); } -void addTargetLocation(ma::Adapt* a, +apf::Field* addTargetLocation(ma::Adapt* a, const char* fieldName) { ma::Mesh* m = a->mesh; @@ -93,6 +132,7 @@ void addTargetLocation(ma::Adapt* a, apf::setVector(paramField , ent, 0, x); } m->end(it); + return paramField; } void addParamCoords(ma::Adapt* a, @@ -116,7 +156,60 @@ void addParamCoords(ma::Adapt* a, m->end(it); } -void colorEntitiesOfDimWithValues(ma::Adapt* a, +void flagEntity(ma::Adapt* a, int dim, const char* fieldName, ma::EntitySet entities) +{ + std::vector converted(entities.begin(), entities.end()); + flagEntity(a, dim, fieldName, &converted[0], converted.size()); +} + +void flagEntityAllDim(ma::Adapt* a, int dim, const char* fieldName, ma::Entity** entities, int size) +{ + std::vector allFaces; + std::vector allEdges; + std::vector allVerts; + + for (int i = 0; i 2){ + ma::Entity* faces[4]; + a->mesh->getDownward(entities[i], 2, faces); + for (ma::Entity* face : faces) allFaces.push_back(face); + } + if (dim > 1){ + ma::Entity* edges[6]; + a->mesh->getDownward(entities[i], 1, edges); + for (ma::Entity* edge : edges) allEdges.push_back(edge); + } + if (dim > 0){ + ma::Entity* verts[4]; + a->mesh->getDownward(entities[i], 0, verts); + for (ma::Entity* vert : verts) allVerts.push_back(vert); + } + } + + flagEntity(a, dim, fieldName, entities, size); + if (dim > 2) flagEntity(a, 2, fieldName, &allFaces[0], allFaces.size()); + if (dim > 1) flagEntity(a, 1, fieldName, &allEdges[0], allEdges.size()); + if (dim > 0) flagEntity(a, 0, fieldName, &allVerts[0], allVerts.size()); +} + +void flagEntity(ma::Adapt* a, int dim, const char* fieldName, ma::Entity** entities, int size) +{ + apf::Field* colorField = getField(a, dim, fieldName); + ma::Entity* entity; + ma::Iterator* it = a->mesh->begin(dim); + while ((entity = a->mesh->iterate(it))){ + double zero = 0.0; + apf::setComponents(colorField, entity, 0, &zero); + } + a->mesh->end(it); + + for (int i=0; i & vals, const char* fieldName) @@ -142,6 +235,7 @@ void colorEntitiesOfDimWithValues(ma::Adapt* a, index++; } m->end(it); + return colorField; } void evaluateFlags(ma::Adapt* a, diff --git a/ma/maDBG.h b/ma/maDBG.h index 16a65d882..4e59877e0 100644 --- a/ma/maDBG.h +++ b/ma/maDBG.h @@ -18,6 +18,7 @@ #include #include +#include namespace ma_dbg { @@ -27,16 +28,32 @@ void writeMesh(ma::Mesh* m, /* 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 useFieldInfo(ma::Adapt* a, const std::function& funcUsingField); -void addTargetLocation(ma::Adapt* a, +apf::Field* addTargetLocation(ma::Adapt* a, const char* fieldName); void addParamCoords(ma::Adapt* a, const char* fieldName); -void colorEntitiesOfDimWithValues(ma::Adapt* a, +void flagEntity(ma::Adapt* a, + int dim, + const char* fieldName, + ma::EntitySet entities); + +void flagEntity(ma::Adapt* a, + int dim, + const char* fieldName, + ma::Entity** entities, + int size); + +void flagEntityAllDim(ma::Adapt* a, + int dim, + const char* fieldName, + ma::Entity** entities, + int size); + +apf::Field* colorEntitiesOfDimWithValues(ma::Adapt* a, int dim, const std::vector & quals, const char* fieldName); diff --git a/ma/maFaceSplit.cc b/ma/maFaceSplit.cc index 73f6da425..2ae39ac77 100644 --- a/ma/maFaceSplit.cc +++ b/ma/maFaceSplit.cc @@ -12,6 +12,7 @@ #include "maFaceSplit.h" #include "maSolutionTransfer.h" #include "maShapeHandler.h" +#include "apfVectorElement.h" namespace ma { @@ -173,18 +174,18 @@ Entity* makeSplitVertOnFace(Adapt* a, Entity* face) SolutionTransfer* st = a->solutionTransfer; // midpoint of face (parametric space) Vector xi(1.0/3.0, 1.0/3.0, 1.0/3.0); - apf::MeshElement* me = apf::createMeshElement(m, face); + apf::MeshElement me(m->getCoordinateField(), face); Vector point; - apf::mapLocalToGlobal(me,xi,point); + apf::mapLocalToGlobal(&me,xi,point); Vector param(0,0,0); if (a->input->shouldTransferParametric) transferParametricOnTriSplit(m, face, xi, param); if (a->input->shouldTransferToClosestPoint) transferToClosestPointOnTriSplit(m, face, xi, param); Entity* vert = buildVertex(a, c, point, param); - st->onVertex(me, xi, vert); - sf->interpolate(me, xi, vert); - apf::destroyMeshElement(me); + st->onVertex(&me, xi, vert); + sf->interpolate(&me, xi, vert); + tagVertexToSnap(a, vert); return vert; } diff --git a/ma/maFaceSwap.cc b/ma/maFaceSwap.cc new file mode 100644 index 000000000..794e50bc0 --- /dev/null +++ b/ma/maFaceSwap.cc @@ -0,0 +1,243 @@ +#include "maAdapt.h" +#include "maSnapper.h" +#include "apfGeometry.h" +#include "maShapeHandler.h" +#include "maDBG.h" + +/** + * \file maFaceSwap.cc + * \brief Definition of maFaceSwap.h file. + * This file contain functions to remove the face between two tetrahedron + * and replace it with a face pointing in the opposite direction or an edge + * without affecting any of the vertices in the cavity. This is done to + * improve the quality of the cavity or to make room for vertices to move. +*/ + +namespace ma { + + Entity* findCommonEdge(Mesh* mesh, Entity* face1, Entity* face2) + { + Entity* face1Edges[3]; + mesh->getDownward(face1, 1, face1Edges); + for (int i=0; i<3; i++) + if (isLowInHigh(mesh, face2, face1Edges[i])) + return face1Edges[i]; + return face1; + } + + void fillAdjVerts(Mesh* mesh, Entity* newTetVerts[4], const Upward& adjTets) + { + apf::Up adjEdges; + mesh->getUp(newTetVerts[0], adjEdges); + int j=1; + for (int i=0; imesh), face(f), improveQuality(improve) + { + cavity.init(a); + numNewTets=0; + } + + void printCavityBefore() + { + EntityArray cavity; + for (int i=0; i<2; i++){ + cavity.append(oldTets[i]); + } + ma_dbg::createCavityMesh(adapt, cavity, "BeforeFaceSwap"); + } + + void printCavityAfter() + { + EntityArray cavity; + for (int i=0; igetDownward(face, 0, faceVerts); + Entity* tetVertsA[4] = {faceVerts[0], faceVerts[1], opp1, opp2}; + if (invalid(tetVertsA)) return true; + Entity* tetVertsB[4] = {faceVerts[1], faceVerts[2], opp1, opp2}; + if (invalid(tetVertsB)) return true; + Entity* tetVertsC[4] = {faceVerts[2], faceVerts[0], opp1, opp2}; + if (invalid(tetVertsC)) return true; + return false; + } + + bool topoCheck() + { + int modelDim = mesh->getModelType(mesh->toModel(face)); + if (modelDim == 2) + return false; + if (getFlag(adapt,face,DONT_SWAP)) + return false; + mesh->getAdjacent(face, 3, oldTets); + if (invalidSwap()) + return false; + return true; + } + + void findNumNewTets() + { + numNewTets = 3; + Entity* faceVerts[3]; + mesh->getDownward(face, 0, faceVerts); + for (int i=0; i<3; i++) { + Entity* face0 = getTetFaceOppositeVert(mesh, oldTets[0], faceVerts[i]); + Entity* face1 = getTetFaceOppositeVert(mesh, oldTets[1], faceVerts[i]); + Vector normal0 = getTriNormal(mesh, face0); + Vector normal1 = getTriNormal(mesh, face1); + if (apf::areClose(normal0, normal1, 1e-10)) { + numNewTets = 0; + commonEdge = findCommonEdge(mesh, face0, face1); + } + } + } + + void buildTwo2Two() + { + printf("\t[Warning]: found two2two face swap, this is unteasted please test\n"); + Model* model = mesh->toModel(oldTets[0]); + Entity* commEdgeVerts[2]; + mesh->getDownward(commonEdge, 0, commEdgeVerts); + + Entity* newTetVerts0[4]; + Entity* newTetVerts1[4]; + newTetVerts0[0] = commEdgeVerts[0]; + newTetVerts1[0] = commEdgeVerts[1]; + + fillAdjVerts(mesh, newTetVerts0, oldTets); + fillAdjVerts(mesh, newTetVerts1, oldTets); + + newTets[0] = buildElement(adapt, model, apf::Mesh::TET, newTetVerts0); + newTets[1] = buildElement(adapt, model, apf::Mesh::TET, newTetVerts1); + } + + void buildTwo2Three() + { + Model* model = mesh->toModel(oldTets[0]); + Entity* newTetVerts[4]; + newTetVerts[0] = getTetVertOppositeTri(mesh, oldTets[0], face); + newTetVerts[1] = getTetVertOppositeTri(mesh, oldTets[1], face); + + Entity* faceEdges[3]; + mesh->getDownward(face, 1, faceEdges); + for (int i=0; i<3; i++) { + Entity* faceEdgeVerts[2]; + mesh->getDownward(faceEdges[i], 0, faceEdgeVerts); + newTetVerts[2]=faceEdgeVerts[0]; + newTetVerts[3]=faceEdgeVerts[1]; + newTets[i] = buildElement(adapt, model, apf::Mesh::TET, newTetVerts); + } + } + + void destroyOldElements() + { + cavity.transfer(oldTets); + destroyElement(adapt, oldTets[0]); + destroyElement(adapt, oldTets[1]); + } + + void cancel() + { + for (int i=0; igetDownward(newTets[i], 1, edges); + for (int e=0; e<6; e++) + if (adapt->sizeField->measure(edges[e]) > MAXLENGTH) + return false; + } + return true; + } + + bool qualityCheck() + { + double qualityToBeat = adapt->input->validQuality; + if (improveQuality) + qualityToBeat = std::max(getWorstQuality(adapt, oldTets), adapt->input->validQuality); + + for (int i=0; ishape->getQuality(newTets[i]); + if (quality < qualityToBeat) + return false; + } + return true; + } + + private: + Adapt* adapt; + Mesh* mesh; + Entity* face; + Cavity cavity; + int numNewTets; + Upward oldTets; + Entity* newTets[3]; + bool improveQuality; + Entity* commonEdge; + }; + + bool runFaceSwap(Adapt* a, Entity* face, bool improveQuality) + { + FaceSwap faceSwap(a, face, improveQuality); + + if (faceSwap.topoCheck() + && faceSwap.geomCheck() + && faceSwap.sizeCheck() + && faceSwap.qualityCheck()) { + faceSwap.destroyOldElements(); + return true; + } + else { + faceSwap.cancel(); + return false; + } + } + +} diff --git a/ma/maFaceSwap.h b/ma/maFaceSwap.h new file mode 100644 index 000000000..d6c58a1af --- /dev/null +++ b/ma/maFaceSwap.h @@ -0,0 +1,8 @@ +#ifndef MA_FACESWAP_H +#define MA_FACESWAP_H + +namespace ma { + bool runFaceSwap(Adapt* a, Entity* face, bool improveQuality=true); +} + +#endif diff --git a/ma/maFixShape.cc b/ma/maFixShape.cc new file mode 100644 index 000000000..fbbf4aadc --- /dev/null +++ b/ma/maFixShape.cc @@ -0,0 +1,456 @@ + +/****************************************************************************** + + Copyright 2013 Scientific Computation Research Center, + Rensselaer Polytechnic Institute. All rights reserved. + + The LICENSE file included with this distribution describes the terms + of the SCOREC Non-Commercial License this program is distributed under. + +*******************************************************************************/ +#include "maFixShape.h" +#include "maSize.h" +#include "maAdapt.h" +#include "maSnap.h" +#include "maSnapper.h" +#include "maOperator.h" +#include "maEdgeSwap.h" +#include "maDoubleSplitCollapse.h" +#include "maSingleSplitCollapse.h" +#include "maFaceSplitCollapse.h" +#include "maFaceSwap.h" +#include "maShortEdgeRemover.h" +#include "maShapeHandler.h" +#include "maBalance.h" +#include "maDBG.h" +#include + +/** + * \file maFixShape.cc + * \brief Definition of maFixShape.h file. + * This file contains functions to improve the quality of tetrahedron in a mesh. + * As described in Li's thesis it will first deterimine if low quality elements + * in the mesh contain short edges and remove them through collapses. Otherwise + * it will collect some properties of the tetrahedron to determine with operations + * are most likely to succeed. +*/ + +namespace ma { + +const double GOODQUALITY2D = .09; + +int markBadQuality(Adapt* a) +{ + Mesh* m = a->mesh; + Iterator* it; + Entity* e; + it = m->begin(m->getDimension()); + int total = 0; + while ((e = m->iterate(it))) { + if (!a->mesh->isOwned(e)) continue; + if (!isSimplex(m->getType(e))) continue; + if (getAndCacheQuality(a, e) >= a->input->goodQuality) continue; + setFlag(a, e, ma::BAD_QUALITY); + total++; + } + m->end(it); + return m->getPCU()->Add(total);; +} + +FixShape::FixShape(Adapt* adapt) : splitCollapse(adapt), doubleSplitCollapse(adapt), faceSplitCollapse(adapt), reposition(adapt), split(adapt) +{ + a = adapt; + mesh = a->mesh; + collapse.Init(a); + edgeSwap = makeEdgeSwap(a); +} + +int FixShape::getTargetDimension() { return 3; } +bool FixShape::shouldApply(Entity* e) { + badTet = e; + return getFlag(a, e, BAD_QUALITY); +} +bool FixShape::requestLocality(apf::CavityOp* o) +{ + Entity* verts[4]; + mesh->getDownward(badTet, 0, verts); + return o->requestLocality(verts, 4); +} + +bool FixShape::collapseEdge(Entity* edge) +{ + if (collapse.setEdge(edge) && + collapse.checkClass() && + collapse.checkTopo() && + collapse.tryBothDirections(a->input->validQuality)) { + collapse.destroyOldElements(); + return true; + } + return false; +} + +bool FixShape::collapseToAdjacent(Entity* edge) +{ + Entity* verts[2]; + mesh->getDownward(edge, 0, verts); + for (int v=0; v<2; v++) { + apf::Up adjEdges; + mesh->getUp(verts[v], adjEdges); + for (int e=0; eisOwned(keep)) continue; + if (mesh->isShared(keep)) continue; + bool alreadyFlagged = getFlag(a, keep, DONT_COLLAPSE); + setFlag(a, keep, DONT_COLLAPSE); + bool success = collapseEdge(adjEdges.e[e]); + if (!alreadyFlagged) clearFlag(a, keep, DONT_COLLAPSE); + if (success) return true; + } + } + return false; +} + +bool FixShape::fixShortEdge(Entity* tet) +{ + Entity* edges[6]; + mesh->getDownward(tet, 1, edges); + for (int i=0; i<6; i++) + if (getAndCacheSize(a, edges[i]) < MINLENGTH) + if (collapseEdge(edges[i])) + { numCollapse++; return true; } + for (int i=0; i<6; i++) + if (getAndCacheSize(a, edges[i]) < MINLENGTH) + if (collapseToAdjacent(edges[i])) + { numCollapse++; return true; } + return false; +} + +bool FixShape::splitReposition(Entity* edge) +{ + if (!split.setEdges(&edge, 1)) + return false; + double worstQuality = getWorstQuality(a,split.getTets()); + split.makeNewElements(); + split.transfer(); + Entity* newVert = split.getSplitVert(0); + reposition.moveToImproveQuality(newVert); + + EntityArray adjacent; + mesh->getAdjacent(newVert, mesh->getDimension(), adjacent); + if (hasWorseQuality(a,adjacent,worstQuality)) { + split.cancel(); + return false; + } + split.destroyOldElements(); + return true; +} + +bool FixShape::isOneLargeAngle(Entity* tet, Entity*& worstTriangle) +{ + Entity* triangles[4]; + mesh->getDownward(tet, 2, triangles); + double worstQuality = 1; + worstTriangle = triangles[0]; + for (int i=0; i<4; i++) { + double quality = getAndCacheQuality(a, triangles[i]); + if (quality < worstQuality) { + worstQuality = quality; + worstTriangle = triangles[i]; + } + } + return worstQuality < GOODQUALITY2D; +} + +Entity* FixShape::getLongestEdge(Entity* edges[3]) +{ + double longestLength = 0; + Entity* longestEdge = edges[0]; + for (int i=0; i<3; i++) { + double length = getAndCacheSize(a, edges[i]); + if (length > longestLength) { + longestLength = length; + longestEdge = edges[i]; + } + } + return longestEdge; +} + +bool FixShape::fixOneLargeAngle(Entity* tet) +{ + if (collapseRegion(tet)) {numRegionCollapse++; return true;} + Entity* worstTriangle; + if (!isOneLargeAngle(tet, worstTriangle)) return false; + Entity* edges[3]; + mesh->getDownward(worstTriangle, 1, edges); + Entity* longestEdge = getLongestEdge(edges); + Entity* oppositeVert = getTriVertOppositeEdge(mesh, worstTriangle, longestEdge); + if (edgeSwap->run(longestEdge)) {numEdgeSwap++; return true;} + if (splitCollapse.run(longestEdge, oppositeVert)) {numEdgeSplitCollapse++; return true;} + for (int i=0; i<3; i++) + if (edges[i] != longestEdge && collapseEdge(edges[i])) + {numCollapse++; return true;} + if (reposition.improveQuality(tet)) {numReposition++; return true;} + return false; +} + +bool FixShape::collapseRegion(Entity* tet) +{ + Entity* interior[4]; + Entity* surface[4]; + Entity* faces[4]; + mesh->getDownward(tet, 2, faces); + int numSurface = 0, numInterior = 0; + for (int f=0; f<4; f++) { + if (isOnModelFace(mesh, faces[f])) + surface[numSurface++] = faces[f]; + else interior[numInterior++] = faces[f]; + } + if (numSurface == 2 && numInterior == 2) { + Entity* interiorEdge = 0; + Entity* surfaceEdge = 0; + Entity* edges[6]; + mesh->getDownward(tet, 1, edges); + for (Entity* edge : edges) { + if (isLowInHigh(mesh, surface[0], edge) && isLowInHigh(mesh, surface[1], edge)) + surfaceEdge = edge; + else if (isLowInHigh(mesh, interior[0], edge) && isLowInHigh(mesh, interior[1], edge)) + interiorEdge = edge; + } + if (interiorEdge==0 || surfaceEdge==0) return false; + if (isOnModelEdge(mesh, surfaceEdge)) return false; + + Model* modelFace = mesh->toModel(surface[0]); + mesh->destroy(tet); + mesh->destroy(surface[0]); + mesh->destroy(surface[1]); + mesh->destroy(surfaceEdge); + + mesh->setModelEntity(interior[0], modelFace); + mesh->setModelEntity(interior[1], modelFace); + mesh->setModelEntity(interiorEdge, modelFace); + return true; + } + return false; +} + +bool FixShape::isTwoLargeAngles(Entity* tet, Entity* problemEnts[4]) +{ + Entity* verts[4]; + mesh->getDownward(tet, 0, verts); + Entity* vert = verts[0]; + Entity* face = getTetFaceOppositeVert(mesh, tet, vert); + double area[4]; + int bit = getTetStats(a, vert, face, tet, problemEnts, area); + return bit==3 || bit==5 || bit==6; +} + +bool FixShape::fixTwoLargeAngles(Entity* tet, Entity* problemEnts[4]) +{ + if (collapseRegion(tet)) {numRegionCollapse++; return true;} + if (edgeSwap->run(problemEnts[0])) {numEdgeSwap++; return true;} + if (edgeSwap->run(problemEnts[1])) {numEdgeSwap++; return true;} + if (splitReposition(problemEnts[0])) {numSplitReposition++; return true;} + if (splitReposition(problemEnts[1])) {numSplitReposition++; return true;} + Entity* faces[4]; + mesh->getDownward(tet, 2, faces); + for (int i=0; i<4; i++) { + if (isLowInHigh(mesh, faces[i], problemEnts[0])) { + Entity* v0 = getTriVertOppositeEdge(mesh, faces[i], problemEnts[0]); + if (splitCollapse.run(problemEnts[0], v0)) {numEdgeSplitCollapse++; return true;} + } + else { + Entity* v0 = getTriVertOppositeEdge(mesh, faces[i], problemEnts[1]); + if (splitCollapse.run(problemEnts[1], v0)) {numEdgeSplitCollapse++; return true;} + } + } + if (doubleSplitCollapse.run(problemEnts)) {numDoubleSplitCollapse++; return true;} + if (reposition.improveQuality(tet)) {numReposition++; return true;} + return false; +} + +bool FixShape::fixThreeLargeAngles(Entity* tet, Entity* problemEnts[4]) +{ + Entity* tetEdges[6]; + mesh->getDownward(tet, 1, tetEdges); + for (int i=0; i<6; i++) + if (!isLowInHigh(mesh, problemEnts[0], tetEdges[i])) { + Entity* verts[2]; + mesh->getDownward(tetEdges[i], 0, verts); + Entity* keep = isLowInHigh(mesh, problemEnts[0], verts[0]) ? verts[0] : verts[1]; + bool alreadyFlagged = getFlag(a, keep, DONT_COLLAPSE); + setFlag(a, keep, DONT_COLLAPSE); + bool success = collapseEdge(tetEdges[i]); + if (!alreadyFlagged) clearFlag(a, keep, DONT_COLLAPSE); + if (success) {numCollapse++; return true;} + } + Entity* faceEdges[3]; + mesh->getDownward(problemEnts[0], 1, faceEdges); + if (edgeSwap->run(faceEdges[0])) {numEdgeSwap++; return true;} + if (edgeSwap->run(faceEdges[1])) {numEdgeSwap++; return true;} + if (edgeSwap->run(faceEdges[2])) {numEdgeSwap++; return true;} + // if (runFaceSwap(a, problemEnts[0], true)) {numFaceSwap++; return true;} + Entity* v = getTriVertOppositeEdge(mesh, problemEnts[2], problemEnts[1]); + if (splitCollapse.run(problemEnts[1], v)) {numEdgeSplitCollapse++; return true;} + if (faceSplitCollapse.run(problemEnts[0], tet)) {numFaceSplitCollapse++; return true;} + if (reposition.moveToImproveQuality(v)) {numReposition++; return true;} + return false; +} + +void FixShape::apply() +{ + clearFlag(a, badTet, BAD_QUALITY); + if (fixShortEdge(badTet)) return; + if (fixOneLargeAngle(badTet)) return; + Entity* problemEnts[4]; + if (isTwoLargeAngles(badTet, problemEnts)) + fixTwoLargeAngles(badTet, problemEnts); + else + fixThreeLargeAngles(badTet, problemEnts); +} + +void FixShape::resetCounters() +{ + numOneShortEdge=numTwoShortEdge=numThreeShortEdge=numMoreShortEdge=0; + numOneLargeAngle=numTwoLargeAngles=numThreeLargeAngles=0; +} +int FixShape::collect(int val) { + return mesh->getPCU()->Add(val); +} +void FixShape::printNumOperations() +{ + print(mesh->getPCU(), "shape operations: \n collapses %13s%d\n edge swaps %12s%d\n reposition %12s%d\n split reposition \t\t\t%d\n double split collapse \t%d\n " + "edge split collapses \t%d\n face split collapses \t%d\n region collapses \t\t\t%d\n face swaps %12s%d\n ", + "", collect(numCollapse), "", collect(numEdgeSwap), "", collect(numReposition), collect(numSplitReposition), collect(numDoubleSplitCollapse), + collect(numEdgeSplitCollapse), collect(numFaceSplitCollapse), collect(numRegionCollapse), "", collect(numFaceSwap)); +} +void FixShape::printBadTypes() +{ + print(mesh->getPCU(), "bad shape types: \n oneShortEdge \t\t%d\n twoShortEdges \t\t%d\n threeShortEdges \t%d\n " + "moreShortEdges \t%d\n oneLargeAngle \t\t%d\n twoLargeAngle \t\t%d\n threeLargeAngle \t%d\n", + collect(numOneShortEdge), collect(numTwoShortEdge), collect(numThreeShortEdge), + collect(numMoreShortEdge), collect(numOneLargeAngle), collect(numTwoLargeAngles), collect(numThreeLargeAngles)); +} + +bool FixShape::isShortEdge(Entity* tet) +{ + Entity* edges[6]; + mesh->getDownward(tet, 1, edges); + int numShortEdges=0; + for (int i=0; i<6; i++) + if (getAndCacheSize(a, edges[i]) < MINLENGTH) + numShortEdges++; + + if (numShortEdges==0) return false; + else if (numShortEdges==1) numOneShortEdge++; + else if (numShortEdges==2) numTwoShortEdge++; + else if (numShortEdges==3) numThreeShortEdge++; + else numMoreShortEdge++; + return true; +} + +void FixShape::printBadShape(Entity* problemTet) +{ + ma_dbg::useFieldInfo(a, [this, &problemTet] { + Entity* worstTri; + Entity* problemEnts[4]; + if (isShortEdge(problemTet)) print(mesh->getPCU(), "Worst is short\n"); + else if (isOneLargeAngle(problemTet, worstTri)) print(mesh->getPCU(), "Worst is one large angle\n"); + else if (isTwoLargeAngles(problemTet, problemEnts)) print(mesh->getPCU(), "Worst is two large angles\n"); + else print(mesh->getPCU(), "Worst is three large angles\n"); + + EntitySet bad; + bad.insert(problemTet); + EntitySet adjacent1 = getNextLayer(a, bad); + EntitySet adjacent2 = getNextLayer(a, adjacent1); + + ma_dbg::flagEntityAllDim(a, 3, "worst_tet", &problemTet, 1); + ma_dbg::flagEntity(a, 3, "worst_tet", bad); + ma_dbg::flagEntity(a, 3, "worst_tet_adj", adjacent1); + ma_dbg::flagEntity(a, 3, "worst_tet_adj2", adjacent2); + + std::vector badFaces; + Iterator* it = mesh->begin(3); + Entity* tet; + while ((tet = mesh->iterate(it))) { + if (!getFlag(a, tet, BAD_QUALITY)) continue; + Entity* faces[4]; + mesh->getDownward(tet, 2, faces); + for (Entity* face : faces) + if (isOnModelFace(mesh, face)) + badFaces.push_back(face); + } + ma_dbg::flagEntity(a, 2, "bad_surface_tets", &badFaces[0], badFaces.size()); + apf::writeVtkFiles("mesh_tets", mesh, 3); + apf::writeVtkFiles("mesh_faces", mesh, 2); + apf::writeVtkFiles("mesh_edges", mesh, 1); + }); +} + +void FixShape::printNumTypes() +{ + resetCounters(); + Entity* tet; + Iterator* it = mesh->begin(3); + // double worstQual = 1; + // Entity* worstShape; + while ((tet = mesh->iterate(it))) { + if (!getFlag(a, tet, BAD_QUALITY)) continue; + // double qual = getAndCacheQuality(a, tet); + // if (qual < worstQual) { + // worstQual = qual; + // worstShape = tet; + // } + + Entity* problemEnts[4]; + Entity* worst; + if (isShortEdge(tet)) continue; + else if (isOneLargeAngle(tet, worst)) + numOneLargeAngle++; + else if (isTwoLargeAngles(tet, problemEnts)) + numTwoLargeAngles++; + else + numThreeLargeAngles++; + } + printBadTypes(); + // printBadShape(worstShape); +} + +void fixElementShapes(Adapt* a) +{ + if ( ! a->input->shouldFixShape) + return; + double t0 = pcu::Time(); + bool oldForce = a->input->shouldForceAdaptation; + a->input->shouldForceAdaptation = false; + int count = markBadQuality(a); + print(a->mesh->getPCU(), "loop %d: of shape correction loop: #bad elements %d", 0, count); + FixShape fixShape(a); + int originalCount = count; + int prev_count; + int iter = 0; + do { + if (!count) break; + double tLoopStart = pcu::Time(); + prev_count = count; + applyOperator(a,&fixShape); + if (a->mesh->getDimension() == 3) + snap(a); + count = markBadQuality(a); + midBalance(a); // balance the mesh to avoid empty parts + iter++; + double tLoopEnd = pcu::Time(); + print(a->mesh->getPCU(), "loop %d: bad shapes went from %d to %d in %f seconds",iter,prev_count,count,tLoopEnd-tLoopStart); + } while(count < prev_count); + a->input->shouldForceAdaptation = oldForce; + double tEnd = pcu::Time(); + print(a->mesh->getPCU(), "bad shapes down from %d to %d in %f seconds",originalCount,count,tEnd-t0); + if (a->input->shouldPrintQuality) { + fixShape.printNumOperations(); + fixShape.printNumTypes(); + printHistogramStats(a); + } + clearFlagFromDimension(a, BAD_QUALITY, a->mesh->getDimension()); +} + +} \ No newline at end of file diff --git a/ma/maFixShape.h b/ma/maFixShape.h new file mode 100644 index 000000000..e806a651a --- /dev/null +++ b/ma/maFixShape.h @@ -0,0 +1,85 @@ +#ifndef MA_FIX_SHAPE +#define MA_FIX_SHAPE + +#include "maMesh.h" +#include "maTables.h" +#include "maOperator.h" +#include "maCollapse.h" +#include "maEdgeSwap.h" +#include "maDoubleSplitCollapse.h" +#include "maSingleSplitCollapse.h" +#include "maFaceSplitCollapse.h" +#include "maFaceSwap.h" +#include "maReposition.h" + +namespace ma { +class Adapt; + +class FixShape : public Operator +{ + public: + Adapt* a; + Mesh* mesh; + Entity* badTet; + Collapse collapse; + SingleSplitCollapse splitCollapse; + DoubleSplitCollapse doubleSplitCollapse; + FaceSplitCollapse faceSplitCollapse; + RepositionVertex reposition; + EdgeSwap* edgeSwap; + Splits split; + + int numCollapse=0; + int numEdgeSwap=0; + int numFaceSwap=0; + int numReposition=0; + int numSplitReposition=0; + int numRegionCollapse=0; + int numEdgeSplitCollapse=0; + int numFaceSplitCollapse=0; + int numDoubleSplitCollapse=0; + + int numOneShortEdge=0; + int numTwoShortEdge=0; + int numThreeShortEdge=0; + int numMoreShortEdge=0; + int numOneLargeAngle=0; + int numTwoLargeAngles=0; + int numThreeLargeAngles=0; + + FixShape(Adapt* adapt); + void apply(); + + int getTargetDimension(); + bool shouldApply(Entity* e); + bool requestLocality(apf::CavityOp* o); + + bool collapseEdge(Entity* edge); + bool collapseToAdjacent(Entity* edge); + double getWorstShape(EntityArray& tets, Entity*& worst); + Entity* getLongestEdge(Entity* edges[3]); + + bool splitReposition(Entity* edge); + bool collapseRegion(Entity* tet); + + bool isShortEdge(Entity* tet); + bool isOneLargeAngle(Entity* tet, Entity*& worstTriangle); + bool isTwoLargeAngles(Entity* tet, Entity* problemEnts[4]); + + bool fixShortEdge(Entity* tet); + bool fixOneLargeAngle(Entity* tet); + bool fixTwoLargeAngles(Entity* tet, Entity* problemEnts[4]); + bool fixThreeLargeAngles(Entity* tet, Entity* problemEnts[4]); + + void resetCounters(); + int collect(int val); + void printNumOperations(); + void printBadTypes(); + void printBadShape(Entity* badTet); + void printNumTypes(); +}; + +void fixElementShapes(Adapt* a); +} + +#endif diff --git a/ma/maInput.cc b/ma/maInput.cc index e153f31b9..0772184cd 100644 --- a/ma/maInput.cc +++ b/ma/maInput.cc @@ -29,16 +29,11 @@ void setDefaultValues(Input* in) in->shouldFixShape = true; in->shouldForceAdaptation = false; in->shouldPrintQuality = true; + in->goodQuality = GOODQUALITY; if (in->mesh->getDimension()==3) - { - in->goodQuality = 0.027; in->maximumEdgeRatio = 2.0; - } - else - { PCU_ALWAYS_ASSERT(in->mesh->getDimension()==2); - //old MA says .04, but that rarely kicks in. - //.2 is strict, but at least quality goes up - in->goodQuality = 0.2; + else { + PCU_ALWAYS_ASSERT(in->mesh->getDimension()==2); //this basically turns off short edge removal... in->maximumEdgeRatio = 100.0; //2D mesh adapt performs better if forceAdapt is on @@ -119,8 +114,8 @@ void validateInput(Input* in) "but shape correction does not support matching yet", in->mesh->getPCU()); if (in->goodQuality < 0.0) rejectInput("negative desired element quality", in->mesh->getPCU()); - if (in->goodQuality > 1.0) - rejectInput("desired element quality greater than one", in->mesh->getPCU()); + if (in->goodQuality > .5) + rejectInput("desired element quality too large", in->mesh->getPCU()); if (in->validQuality < 0.0) rejectInput("negative minimum element quality", in->mesh->getPCU()); if (in->maximumImbalance < 1.0) diff --git a/ma/maMesh.cc b/ma/maMesh.cc index 3e0f370e8..705026b82 100644 --- a/ma/maMesh.cc +++ b/ma/maMesh.cc @@ -328,6 +328,11 @@ bool isOnModelFace(Mesh* m, Entity* e) return m->getModelType(m->toModel(e))==2; } +bool isInModelRegion(Mesh* m, Entity* e) +{ + return m->getModelType(m->toModel(e))==3; +} + Vector getTriNormal(Mesh* m, Entity** v) { Vector x[3]; diff --git a/ma/maMesh.h b/ma/maMesh.h index f5409d2b8..cabff9763 100644 --- a/ma/maMesh.h +++ b/ma/maMesh.h @@ -100,6 +100,7 @@ Entity* findTriFromVerts(Mesh* m, Entity** v); bool isOnModelEdge(Mesh* m, Entity* e); bool isOnModelFace(Mesh* m, Entity* e); +bool isInModelRegion(Mesh* m, Entity* e); Vector getTriNormal(Mesh* m, Entity** v); Vector getTriNormal(Mesh* m, Entity* e); diff --git a/ma/maQuality.cc b/ma/maQuality.cc index 85f465116..23878625c 100644 --- a/ma/maQuality.cc +++ b/ma/maQuality.cc @@ -17,9 +17,19 @@ #include "maShapeHandler.h" #include "maShape.h" #include +#include "apfVectorElement.h" namespace ma { +bool isTetValid(Mesh* m, Entity* tet) +{ + Vector v[4]; + ma::getVertPoints(m,tet,v); + if ((cross((v[1] - v[0]), (v[2] - v[0])) * (v[3] - v[0])) < 0) + return false; + return true; +} + bool areTetsValid(Mesh* m, EntityArray& tets) { Vector v[4]; @@ -35,10 +45,9 @@ bool areTetsValid(Mesh* m, EntityArray& tets) class FixedMetricIntegrator : public apf::Integrator { public: - FixedMetricIntegrator(Mesh* inMesh, const Matrix& inQ): + FixedMetricIntegrator(const Matrix& inQ): Integrator(1), measurement(0), - mesh(inMesh), Q(inQ) { dimension = 0; @@ -47,16 +56,12 @@ class FixedMetricIntegrator : public apf::Integrator virtual void inElement(apf::MeshElement* me) { dimension = apf::getDimension(me); - elem = apf::createElement(mesh->getCoordinateField(), me); - } - virtual void outElement() - { - apf::destroyElement(elem); + this->elem = me; } virtual void atPoint(Vector const& p, double w, double) { Matrix J; - apf::getJacobian(apf::getMeshElement(elem),p,J); + apf::getJacobian(elem,p,J); /* transforms the rows of J, the differential tangent vectors, into the metric space, then uses the generalized determinant */ double dV2 = apf::getJacobianDeterminant(J*Q,dimension); @@ -64,19 +69,17 @@ class FixedMetricIntegrator : public apf::Integrator } double measurement; private: - Mesh* mesh; Matrix Q; int dimension; - apf::Element* elem; + apf::MeshElement* elem; }; static double qMeasure(Mesh* mesh, Entity* e, const Matrix& Q) { - FixedMetricIntegrator integrator(mesh, Q); - apf::MeshElement* me = apf::createMeshElement(mesh, e); - integrator.process(me); - apf::destroyMeshElement(me); + FixedMetricIntegrator integrator(Q); + apf::MeshElement me(mesh->getCoordinateField(), e); + integrator.process(&me); return integrator.measurement; } @@ -94,15 +97,14 @@ static Matrix getMetricWithMaxJacobean(Mesh* m, SizeField* sf, double maxJ = -1.0; for (int i = 0; i < nd; i++) { - apf::MeshElement* me = createMeshElement(m, dv[i]); + apf::MeshElement me(m->getCoordinateField(), dv[i]); Matrix currentQ; - sf->getTransform(me, Vector(0.0, 0.0, 0.0), currentQ); + sf->getTransform(&me, Vector(0.0, 0.0, 0.0), currentQ); double currentJ = apf::getJacobianDeterminant(currentQ, dim); if (currentJ > maxJ) { maxJ = currentJ; Q = currentQ; } - apf::destroyMeshElement(me); } return Q; } @@ -117,10 +119,9 @@ double measureTriQuality(Mesh* m, SizeField* f, Entity* tri, bool useMax) if (useMax) Q = getMetricWithMaxJacobean(m, f, tri); else { - apf::MeshElement* me = createMeshElement(m, tri); + apf::MeshElement me(m->getCoordinateField(), tri); Vector xi(1./3., 1./3., 1./3.); - f->getTransform(me, xi, Q); - apf::destroyMeshElement(me); + f->getTransform(&me, xi, Q); } Entity* e[3]; @@ -146,10 +147,9 @@ double measureTetQuality(Mesh* m, SizeField* f, Entity* tet, bool useMax) if (useMax) Q = getMetricWithMaxJacobean(m, f, tet); else { - apf::MeshElement* me = createMeshElement(m, tet); + apf::MeshElement me(m->getCoordinateField(), tet); Vector xi(0.25, 0.25, 0.25); - f->getTransform(me, xi, Q); - apf::destroyMeshElement(me); + f->getTransform(&me, xi, Q); } Entity* e[6]; @@ -181,6 +181,15 @@ double measureElementQuality(Mesh* m, SizeField* f, Entity* e, bool useMax) return table[m->getType(e)](m,f,e,useMax); } +double getAndCacheQuality(Adapt* a, Entity* e) +{ + if (a->mesh->hasTag(e, a->qualityCache)) + return getCachedQuality(a, e); + double quality = a->shape->getQuality(e); + setCachedQuality(a, e, quality); + return quality; +} + double getWorstQuality(Adapt* a, Entity** e, size_t n) { PCU_ALWAYS_ASSERT(n); diff --git a/ma/maRefine.cc b/ma/maRefine.cc index 915d8481a..155ed9379 100644 --- a/ma/maRefine.cc +++ b/ma/maRefine.cc @@ -19,6 +19,7 @@ #include "maLayer.h" #include #include +#include "apfVectorElement.h" namespace ma { @@ -135,18 +136,18 @@ Entity* makeSplitVert(Refine* r, Entity* edge) SolutionTransfer* st = a->solutionTransfer; /* midpoint of [-1,1] */ Vector xi(0,0,0); - apf::MeshElement* me = apf::createMeshElement(m,edge); + apf::MeshElement me(m->getCoordinateField(),edge); Vector point; - apf::mapLocalToGlobal(me,xi,point); + apf::mapLocalToGlobal(&me,xi,point); Vector param(0,0,0); //prevents uninitialized values if (a->input->shouldTransferParametric) transferParametricOnEdgeSplit(m,edge,0.5,param); if (a->input->shouldTransferToClosestPoint) transferToClosestPointOnEdgeSplit(m,edge,0.5,param); Entity* vert = buildVertex(a,c,point,param); - st->onVertex(me,xi,vert); - sf->interpolate(me,xi,vert); - apf::destroyMeshElement(me); + st->onVertex(&me,xi,vert); + sf->interpolate(&me,xi,vert); + tagVertexToSnap(a,vert); return vert; } diff --git a/ma/maReposition.cc b/ma/maReposition.cc index 42be04314..c60a26a28 100644 --- a/ma/maReposition.cc +++ b/ma/maReposition.cc @@ -1,4 +1,7 @@ #include "maReposition.h" +#include "maAdapt.h" +#include "maShape.h" +#include "maShapeHandler.h" #include #include #include @@ -8,9 +11,183 @@ #include #include +#include + +/** + * \file maReposition.cc + * \brief Definition of maReposition.h file. + * This file contains functions to move a vertex. If given a target it will reject the + * reposition if invalid elements are created. It also conatins a function that will move + * the vertex to a location that will localy maximize the quality of the adjactent elements. + * If the vertex is on a model boundary then it will keep that vertex on the model boundary + * as locations are using a golden-section search. +*/ namespace ma { +RepositionVertex::RepositionVertex(Adapt* a) : adapt(a), mesh(a->mesh) +{ +} + +bool RepositionVertex::init(Entity* vertex) +{ + this->invalid.n = 0; + this->adjEdges.n = 0; + this->adjacentElements.setSize(0); + this->vertex = vertex; + this->prevPosition = getPosition(mesh, vertex); + mesh->getAdjacent(vertex, mesh->getDimension(), adjacentElements); + for (size_t i = 0; i < adjacentElements.getSize(); ++i) + if (!isSimplex(mesh->getType(adjacentElements[i]))) return false; + mesh->getUp(vertex, adjEdges); + clearAdjCache(); + return true; +} + +void RepositionVertex::clearAdjCache() +{ + for (size_t i = 0; i < adjacentElements.getSize(); ++i) + mesh->removeTag(adjacentElements[i], adapt->qualityCache); + for (int i = 0; i < adjEdges.n; ++i) + mesh->removeTag(adjEdges.e[i], adapt->sizeCache); + apf::Adjacent adjTri; + mesh->getAdjacent(vertex, 2, adjTri); + for (size_t i = 0; i < adjTri.getSize(); ++i) { + mesh->removeTag(adjTri[i], adapt->qualityCache); + mesh->removeTag(adjTri[i], adapt->sizeCache); + } +} + +void RepositionVertex::findInvalid() +{ + 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(adapt, adjacentElements[i], LAYER)) continue; + + double quality; + if (mesh->getType(adjacentElements[i]) == apf::Mesh::TET && !isTetValid(mesh, adjacentElements[i])) + quality = -1; + else + quality = adapt->shape->getQuality(adjacentElements[i]); + + if (quality < adapt->input->validQuality) invalid.e[invalid.n++] = adjacentElements[i]; + } +} + +bool RepositionVertex::move(Entity* vertex, Vector target) +{ + if (!init(vertex)) return false; + Vector pos = getPosition(mesh, vertex); + if (pos[0] == target[0] && pos[1] == target[1] && pos[2] == target[2]) return true; + mesh->setPoint(vertex, 0, target); + findInvalid(); + if (invalid.n == 0) { + Vector xi(0,0,0); + me.init(mesh->getCoordinateField(), vertex, 0); + adapt->solutionTransfer->onVertex(&me, xi, vertex); + return true; + } + cancel(vertex); + return false; +} + +double RepositionVertex::findWorstShape(Vector position) +{ + Model* gm = mesh->toModel(vertex); + if (mesh->getModelType(gm) < 3 && adapt->input->shouldSnap) { + Vector p; + mesh->getParamOn(gm,vertex,p); + mesh->snapToModel(gm,p,position); + mesh->setPoint(vertex,0,position); + } + else mesh->setPoint(vertex, 0, position); + double worst = 1; + for (size_t i = 0; i < adjacentElements.getSize(); ++i) { + double quality = adapt->shape->getQuality(adjacentElements[i]); + if (quality < worst) worst = quality; + } + return worst; +} + +double goldenSearch(const std::function &f, Vector left, Vector right, double tol) +{ + const double M_GOLDEN_RATIO = (1.0 + std::sqrt(5.0)) / 2.0; + double leftValue; + int max = 50; + + do { + Vector delta_left = left + (right - left) * M_GOLDEN_RATIO; + Vector delta_right = right + (left - right) * M_GOLDEN_RATIO; + leftValue = f(delta_left); + + if (leftValue < f(delta_right)) + right = delta_right; + else + left = delta_left; + + if (--max < 0) break; + } while ((left-right).getLength() > tol); + return leftValue; +} + +bool RepositionVertex::moveToImproveQuality(Entity* vertex) +{ + if (mesh->getModelType(mesh->toModel(vertex)) < 2) return false; + if (!init(vertex)) return false; + + Vector center = modelCenter(); + Vector target = center + (prevPosition - center)/4; + const auto getQuality = [this](const Vector& pos) { return this->findWorstShape(pos); }; + double startQuality = getQuality(prevPosition); + worstQuality = goldenSearch(getQuality, prevPosition, target, 0.000001); + if (startQuality > worstQuality) { + mesh->setPoint(vertex, 0, prevPosition); + worstQuality = startQuality; + } + else { + Vector xi(0,0,0); + me.init(mesh->getCoordinateField(), vertex, 0); + adapt->solutionTransfer->onVertex(&me, xi, vertex); + } + return worstQuality > adapt->input->goodQuality; +} + +bool RepositionVertex::improveQuality(Entity* tet) +{ + Entity* verts[4]; + mesh->getDownward(tet, 0, verts); + for (int i=0; i<4; i++) + if (moveToImproveQuality(verts[i])) return true; + return false; +} + +Vector RepositionVertex::modelCenter() +{ + Model* vertModel = mesh->toModel(vertex); + Vector avg(0,0,0); + for (int i=0; itoModel(adjEdges.e[i])) continue; + Entity* opp = getEdgeVertOppositeVert(mesh, adjEdges.e[i], vertex); + avg += getPosition(mesh, opp); + } + avg = avg / adjEdges.n; + return avg; +} + +void RepositionVertex::cancel(Entity* vertex) +{ + PCU_ALWAYS_ASSERT(this->vertex == vertex); + mesh->setPoint(vertex, 0, prevPosition); +} + +apf::Up& RepositionVertex::getInvalid() +{ + return invalid; +} + typedef mth::AD AD; typedef mth::Vector ADVec; diff --git a/ma/maReposition.h b/ma/maReposition.h index 87f603926..1515c37ea 100644 --- a/ma/maReposition.h +++ b/ma/maReposition.h @@ -2,9 +2,42 @@ #define MA_REPOSITION_H #include "maMesh.h" +#include "ma.h" +#include "apfVectorElement.h" namespace ma { +class RepositionVertex +{ + public: + RepositionVertex(Adapt* a); + bool move(Entity* vertex, Vector target); + bool moveToImproveQuality(Entity* vertex); + void moveToImproveShortEdges(Entity* vertex); + bool improveQuality(Entity* tet); + void cancel(Entity* vertex); + apf::Up& getInvalid(); + + private: + Adapt* adapt; + Mesh* mesh; + Entity* vertex; + Vector prevPosition; + Upward adjacentElements; + apf::Up invalid; + apf::Up adjEdges; + double worstQuality; + double startingQuality; + apf::MeshElement me; + + void findInvalid(); + void clearAdjCache(); + Vector modelCenter(); + bool init(Entity* vertex); + double findWorstShape(Vector position); + double findShortestEdge(Vector position); +}; + bool repositionVertex(Mesh* m, Entity* v, int max_iters, double initial_speed); diff --git a/ma/maShape.cc b/ma/maShape.cc index dcd0eec17..1fdbf5aeb 100644 --- a/ma/maShape.cc +++ b/ma/maShape.cc @@ -25,129 +25,6 @@ namespace ma { -/* projects vertex 3 onto the plane - of the bottom triangle and returns - the zone in which it lands as a bit code. - Each bit indicates whether the area coordinate - of that vertex is positive. -*/ - -int getSliverCode( - Adapt* a, - Entity* tet) -{ - SizeField* sf = a->sizeField; - Mesh* m = a->mesh; - Matrix J,Q; - apf::MeshElement* me = apf::createMeshElement(m,tet); - Vector center(.25,.25,.25); - apf::getJacobian(me,center,J); - sf->getTransform(me,center,Q); - J = J*Q; //Jacobian in metric space - apf::destroyMeshElement(me); - int code = 0; - // check first face - Entity* fs[4]; - m->getDownward(tet, 2, fs); - double f0Qual = a->shape->getQuality(fs[0]); - if ((f0Qual*f0Qual*f0Qual > a->input->goodQuality*a->input->goodQuality)) { - // if its okay, use it for projection - Vector v03 = J[2]; - J[2] = apf::cross(J[0],J[1]); //face normal towards v[3] - Vector projected = v03 - apf::project(v03,J[2]); //v[3] projected to face - Matrix inverseMap = apf::invert(apf::transpose(J)); - Vector basisPoint = inverseMap * projected; - Vector areaPoint(1-basisPoint[0]-basisPoint[1], - basisPoint[0], - basisPoint[1]); - for (int i=0; i < 3; ++i) - if (areaPoint[i] > 0) - code |= (1< -0.10 && areaPoint[i] < 0.10) - code |= ((1< 0) - code |= ((1< -0.20 && areaPoint[i] < 0.20) - code |= ((1<> 6) & 1) - return table2d[(code >> 7) & 3][(code >> 9) & 3]; - else - return table[code & 7][(code >> 3) & 7]; -} - -struct IsBadQuality : public Predicate -{ - IsBadQuality(Adapt* a_):a(a_) {} - bool operator()(Entity* e) - { - return a->shape->getQuality(e) < a->input->goodQuality; - } - Adapt* a; -}; - -int markBadQuality(Adapt* a) -{ - IsBadQuality p(a); - return markEntities(a, a->mesh->getDimension(), p, BAD_QUALITY, OK_QUALITY); -} - -void unMarkBadQuality(Adapt* a) -{ - Mesh* m = a->mesh; - Iterator* it; - Entity* e; - it = m->begin(m->getDimension()); - while ((e = m->iterate(it))) { - if (getFlag(a, e, ma::BAD_QUALITY)) - clearFlag(a, e, ma::BAD_QUALITY); - } - m->end(it); -} - double getMinQuality(Adapt* a) { PCU_ALWAYS_ASSERT(a); @@ -168,528 +45,6 @@ double getMinQuality(Adapt* a) return m->getPCU()->Min(minqual); } -class ShortEdgeFixer : public Operator -{ - public: - ShortEdgeFixer(Adapt* a): - remover(a) - { - adapter = a; - mesh = a->mesh; - sizeField = a->sizeField; - shortEdgeRatio = a->input->maximumEdgeRatio; - nr = nf = 0; - element = 0; - } - virtual ~ShortEdgeFixer() - { - } - virtual int getTargetDimension() {return mesh->getDimension();} - virtual bool shouldApply(Entity* e) - { - if ( ! getFlag(adapter,e,BAD_QUALITY)) - return false; - element = e; - Downward edges; - int n = mesh->getDownward(element,1,edges); - double l[6] = {}; - for (int i=0; i < n; ++i) - l[i] = sizeField->measure(edges[i]); - double maxLength; - double minLength; - Entity* shortEdge; - maxLength = minLength = l[0]; - shortEdge = edges[0]; - for (int i=1; i < n; ++i) - { - if (l[i] > maxLength) maxLength = l[i]; - if (l[i] < minLength) - { - minLength = l[i]; - shortEdge = edges[i]; - } - } - if ((maxLength/minLength) < shortEdgeRatio) - { - clearFlag(adapter,element,BAD_QUALITY); - return false; - } - remover.setEdge(shortEdge); - return true; - } - virtual bool requestLocality(apf::CavityOp* o) - { - return remover.requestLocality(o); - } - virtual void apply() - { - if (remover.run()) - ++nr; - else - { - ++nf; - clearFlag(adapter,element,BAD_QUALITY); - } - } - private: - Adapt* adapter; - Mesh* mesh; - Entity* element; - SizeField* sizeField; - ShortEdgeRemover remover; - double shortEdgeRatio; - public: - int nr; - int nf; -}; - -class TetFixerBase -{ - public: - virtual void setTet(Entity** v) = 0; - virtual bool requestLocality(apf::CavityOp* o) = 0; - virtual bool run() = 0; -}; - -class FixBySwap : public TetFixerBase -{ - public: - FixBySwap(Adapt* a): - adapter(a) - { - mesh = a->mesh; - edgeSwap = makeEdgeSwap(a); - nes = nf = numToTry = 0; - edges[0] = 0; - edges[1] = 0; - edges[2] = 0; - } - ~FixBySwap() - { - delete edgeSwap; - } - virtual void setTet(Entity** v) - { - Entity* tet = apf::findElement(mesh, apf::Mesh::TET, v); - PCU_ALWAYS_ASSERT(tet); - match = matchSliver(adapter, tet); - Entity* dv[4]; - mesh->getDownward(tet, 0, dv); - Entity* rv[4]; - rotateTet(dv,match.rotation,rv); - - enum { EDGE_EDGE, FACE_VERT }; - if (match.code_index==EDGE_EDGE) { - Entity* ev[2]; - ev[0] = rv[0]; ev[1] = rv[2]; - edges[0] = findUpward(mesh, apf::Mesh::EDGE, ev); - ev[0] = rv[1]; ev[1] = rv[3]; - edges[1] = findUpward(mesh, apf::Mesh::EDGE, ev); - numToTry = 2; - } - else - { - PCU_ALWAYS_ASSERT(match.code_index==FACE_VERT); - apf::findTriDown(mesh,rv,edges); - numToTry = 3; - } - } - virtual bool requestLocality(apf::CavityOp* o) - { - return o->requestLocality(edges, numToTry); - } - virtual bool run() - { - for (int i=0; i < numToTry; ++i) - if (edgeSwap->run(edges[i])) - { - ++nes; - return true; - } - ++nf; - return false; - } - private: - Adapt* adapter; - Mesh* mesh; - Entity* edges[3]; - EdgeSwap* edgeSwap; - CodeMatch match; - int numToTry; - int nes; - int nf; -}; - -class EdgeVertFixer : public TetFixerBase -{ -public: - EdgeVertFixer(Adapt* a): - singleSplitCollapse(a) - { - mesh = a->mesh; - edgeSwap = makeEdgeSwap(a); - nes = nssc = nf = 0; - edge = 0; - oppVert = 0; - } - ~EdgeVertFixer() - { - delete edgeSwap; - } - virtual void setTet(Entity** v) - { - /* In this template, the edge v[0]--v[1] and vert v[3] - are too close*/ - edge = apf::findElement(mesh, apf::Mesh::EDGE, v); - oppVert = v[3]; - verts[0] = v[0]; - verts[1] = v[1]; - verts[2] = v[3]; - } - virtual bool requestLocality(apf::CavityOp* o) - { - /* by requesting locality for all the verts we can be sure - * that all the desired entities for this operator are local */ - return o->requestLocality(verts,3); - } - virtual bool run() { - if (edgeSwap->run(edge)) { - ++nes; - return true; - } - if (singleSplitCollapse.run(edge, oppVert)) - { - ++nssc; - return true; - } - ++nf; - return false; - } -private: - Mesh* mesh; - Entity* verts[3]; - Entity *edge, *oppVert; - SingleSplitCollapse singleSplitCollapse; - EdgeSwap* edgeSwap; -public: - int nes; /* number of edge swaps done */ - int nssc; /* number of SSCs done */ - int nf; /* number of failures */ -}; - -class FaceVertFixer : public TetFixerBase -{ - public: - FaceVertFixer(Adapt* a): - faceSplitCollapse(a) - { - mesh = a->mesh; - edgeSwap = makeEdgeSwap(a); - nes = nf = nfsc = 0; - edges[0] = 0; - edges[1] = 0; - edges[2] = 0; - verts[0] = 0; - verts[1] = 0; - verts[2] = 0; - verts[3] = 0; - face = 0; - oppVert = 0; - tet = 0; - } - ~FaceVertFixer() - { - delete edgeSwap; - } - virtual void setTet(Entity** v) - { -/* in this template, the bottom face and v[3] - are too close, the key edges are those that bound - face v(0,1,2) */ - apf::findTriDown(mesh,v,edges); - tet = apf::findElement(mesh, apf::Mesh::TET, v); - oppVert = v[3]; - verts[0] = v[0]; - verts[1] = v[1]; - verts[2] = v[2]; - verts[3] = v[3]; - } - virtual bool requestLocality(apf::CavityOp* o) - { - /* by requesting locality for all the verts we can be sure - * that all the desired entities for this operator are local */ - return o->requestLocality(verts,4); - } - virtual bool run() - { - for (int i=0; i < 3; ++i) - if (edgeSwap->run(edges[i])) - { - ++nes; - return true; - } - face = apf::findUpward(mesh, apf::Mesh::TRIANGLE, edges); - if (faceSplitCollapse.run(face, tet)) - { - ++nfsc; - return true; - } - ++nf; - return false; - } - private: - Mesh* mesh; - Entity* edges[3]; - Entity* verts[4]; - Entity *face, *oppVert; - Entity* tet; - FaceSplitCollapse faceSplitCollapse; - EdgeSwap* edgeSwap; - public: - int nes; /* number of edge swaps done */ - int nfsc; /* number of FSCs done */ - int nf; /* number of failures */ -}; - -class EdgeEdgeFixer : public TetFixerBase -{ - public: - EdgeEdgeFixer(Adapt* a): - doubleSplitCollapse(a) - { - mesh = a->mesh; - edgeSwap = makeEdgeSwap(a); - nes = ndsc = nf = 0; - sf = a->sizeField; - edges[0] = 0; - edges[1] = 0; - } - ~EdgeEdgeFixer() - { - delete edgeSwap; - } - virtual void setTet(Entity** v) - { -/* in this template, the v[0]-v[2] amd v[1]-v[3] - edges are too close. */ - Entity* ev[2]; - ev[0] = v[0]; ev[1] = v[2]; - edges[0] = findUpward(mesh, apf::Mesh::EDGE, ev); - ev[0] = v[1]; ev[1] = v[3]; - edges[1] = findUpward(mesh, apf::Mesh::EDGE, ev); - } - virtual bool requestLocality(apf::CavityOp* o) - { - return o->requestLocality(edges,2); - } - virtual bool run() - { - for (int i=0; i < 2; ++i) - if (edgeSwap->run(edges[i])) - { - ++nes; - return true; - } - if (doubleSplitCollapse.run(edges)) - { - ++ndsc; - return true; - } - ++nf; - return false; - } - private: - Mesh* mesh; - Entity* edges[2]; - EdgeSwap* edgeSwap; - DoubleSplitCollapse doubleSplitCollapse; - SizeField* sf; - public: - int nes; - int ndsc; - int nf; -}; - -class LargeAngleTetFixer : public Operator -{ - public: - LargeAngleTetFixer(Adapt* a): - edgeEdgeFixer(a), - edgeVertFixer(a), - faceVertFixer(a) - { - adapter = a; - mesh = a->mesh; - tet = 0; - fixer = 0; - } - virtual ~LargeAngleTetFixer() - { - } - virtual int getTargetDimension() {return 3;} - enum { EDGE_EDGE, FACE_VERT, EDGE_VERT, VERT_VERT }; - virtual bool shouldApply(Entity* e) - { - if ( ! getFlag(adapter,e,BAD_QUALITY)) - return false; - tet = e; - CodeMatch match = matchSliver(adapter,e); - if (match.code_index==EDGE_EDGE) { - fixer = &edgeEdgeFixer; - } else if (match.code_index==FACE_VERT) { - fixer = &faceVertFixer; - } else if (match.code_index==EDGE_VERT) { - fixer = &edgeVertFixer; - } else if (match.code_index==VERT_VERT) { - fixer = &faceVertFixer; - } - Entity* v[4]; - mesh->getDownward(e,0,v); - Entity* rv[4]; - rotateTet(v,match.rotation,rv); - fixer->setTet(rv); - return true; - } - virtual bool requestLocality(apf::CavityOp* o) - { - return fixer->requestLocality(o); - } - virtual void apply() - { - if ( ! fixer->run()) - clearFlag(adapter,tet,BAD_QUALITY); - } - private: - Adapt* adapter; - Mesh* mesh; - Entity* tet; - TetFixerBase* fixer; - public: - EdgeEdgeFixer edgeEdgeFixer; - EdgeVertFixer edgeVertFixer; - FaceVertFixer faceVertFixer; -}; - -class LargeAngleTetAligner : public Operator -{ - public: - LargeAngleTetAligner(Adapt* a): - fixer(a) - { - adapter = a; - mesh = a->mesh; - tet = 0; - } - virtual ~LargeAngleTetAligner() - { - } - virtual int getTargetDimension() {return 3;} - virtual bool shouldApply(Entity* e) - { - if ( ! getFlag(adapter,e,BAD_QUALITY)) - return false; - tet = e; - /* PCU_ALWAYS_ASSERT(mesh->getType(e) == apf::Mesh::TET); */ - enum { EDGE_EDGE, FACE_VERT }; - CodeMatch match = matchSliver(adapter,e); - if (match.code_index==EDGE_EDGE) { - clearFlag(adapter,tet,BAD_QUALITY); - return false; - } - /* else */ - /* { PCU_ALWAYS_ASSERT(match.code_index==FACE_VERT); */ - /* fixer = &faceVertFixer; */ - /* } */ - Entity* v[4]; - mesh->getDownward(e,0,v); - fixer.setTet(v); - return true; - } - virtual bool requestLocality(apf::CavityOp* o) - { - return fixer.requestLocality(o); - } - virtual void apply() - { - if ( ! fixer.run()) - clearFlag(adapter,tet,BAD_QUALITY); - } - private: - Adapt* adapter; - Mesh* mesh; - Entity* tet; - FixBySwap fixer; -}; - -class LargeAngleTriFixer : public Operator -{ - public: - LargeAngleTriFixer(Adapt* a) - { - adapter = a; - mesh = a->mesh; - edgeSwap = makeEdgeSwap(a); - ns = nf = 0; - tri = 0; - edge = 0; - } - virtual ~LargeAngleTriFixer() - { - delete edgeSwap; - } - virtual int getTargetDimension() {return 2;} - virtual bool shouldApply(Entity* e) - { - if ( ! getFlag(adapter,e,BAD_QUALITY)) - return false; - tri = e; - // get the metric Q for angle computations - SizeField* sf = adapter->sizeField; - Matrix Q; - apf::MeshElement* me = apf::createMeshElement(mesh, tri); - Vector center(1./3.,1./3.,1./3.); - sf->getTransform(me,center,Q); - apf::destroyMeshElement(me); - - // pick the edge opposite to the largest angle (in metric) for swap - Entity* edges[3]; - mesh->getDownward(e,1,edges); - double minCos = 1.0; - for (int i = 0; i < 3; i++) { - Entity* current = edges[i%3]; - Entity* next = edges[(i+1)%3]; - double cosAngle = apf::computeCosAngle(mesh, tri, current, next, Q); - if (cosAngle < minCos) { - minCos = cosAngle; - edge = edges[(i+2)%3]; - } - } - return true; - } - virtual bool requestLocality(apf::CavityOp* o) - { - return o->requestLocality(&edge,1); - } - virtual void apply() - { - if (edgeSwap->run(edge)) - { - ++ns; - return; - } - ++nf; - clearFlag(adapter,tri,BAD_QUALITY); - } - private: - Adapt* adapter; - Mesh* mesh; - Entity* tri; - Entity* edge; - EdgeSwap* edgeSwap; - int ns; - int nf; -}; - class QualityImprover2D : public Operator { public: @@ -737,63 +92,12 @@ class QualityImprover2D : public Operator int nf; }; -static double fixShortEdgeElements(Adapt* a) -{ - double t0 = pcu::Time(); - ShortEdgeFixer fixer(a); - applyOperator(a,&fixer); - double t1 = pcu::Time(); - return t1 - t0; -} - -static void fixLargeAngleTets(Adapt* a) -{ - LargeAngleTetFixer fixer(a); - applyOperator(a,&fixer); -} - -static void fixLargeAngleTris(Adapt* a) -{ - LargeAngleTriFixer fixer(a); - applyOperator(a,&fixer); -} - -static void alignLargeAngleTets(Adapt* a) -{ - LargeAngleTetAligner aligner(a); - applyOperator(a,&aligner); -} - -static void alignLargeAngleTris(Adapt* a) -{ - LargeAngleTriFixer aligner(a); - applyOperator(a,&aligner); -} - static void improveQualities2D(Adapt* a) { QualityImprover2D improver(a); applyOperator(a, &improver); } -static double fixLargeAngles(Adapt* a) -{ - double t0 = pcu::Time(); - if (a->mesh->getDimension()==3) - fixLargeAngleTets(a); - else - fixLargeAngleTris(a); - double t1 = pcu::Time(); - return t1 - t0; -} - -static void alignLargeAngles(Adapt* a) -{ - if (a->mesh->getDimension()==3) - alignLargeAngleTets(a); - else - alignLargeAngleTris(a); -} double improveQualities(Adapt* a) { @@ -806,78 +110,12 @@ double improveQualities(Adapt* a) return t1 - t0; } -void fixElementShapes(Adapt* a) -{ - if ( ! a->input->shouldFixShape) - return; - double t0 = pcu::Time(); - int count = markBadQuality(a); - int originalCount = count; - int prev_count; - double time; - int iter = 0; - do { - if ( ! count) - break; - prev_count = count; - print(a->mesh->getPCU(), "--iter %d of shape correction loop: #bad elements %d", iter, count); - time = fixLargeAngles(a); - /* We need to snap the new verts as soon as they are - * created (to avoid future problems). At the moment - * new verts are created only during 3D mesh adapt, so - * we only run a bulk snap operation if the mesh is 3D. - */ - if (a->mesh->getDimension() == 3) - snap(a); - count = markBadQuality(a); - print(a->mesh->getPCU(), "--fixLargeAngles in %f seconds: #bad elements %d", time,count); - time = fixShortEdgeElements(a); - count = markBadQuality(a); - print(a->mesh->getPCU(), "--fixShortEdgeElements in %f seconds: #bad elements %d", time,count); - if (count >= prev_count) - unMarkBadQuality(a); // to make sure markEntities does not complain! - // balance the mesh to avoid empty parts - midBalance(a); - print(a->mesh->getPCU(), "--percent change in number of bad elements %f", - ((double) prev_count - (double) count) / (double) prev_count); - iter++; - } while(count < prev_count); - double t1 = pcu::Time(); - print(a->mesh->getPCU(), "bad shapes down from %d to %d in %f seconds", - originalCount,count,t1-t0); -} - -void alignElements(Adapt* a) -{ - int max_iter = 5; - if ( ! a->input->shouldFixShape) - return; - double t0 = pcu::Time(); - int count = markBadQuality(a); - int originalCount = count; - int prev_count; - int i = 0; - do { - if ( ! count) - break; - prev_count = count; - alignLargeAngles(a); - count = markBadQuality(a); - ++i; - if (count >= prev_count || i >= max_iter) - unMarkBadQuality(a); - } while(count < prev_count && i < max_iter); - - double t1 = pcu::Time(); - print(a->mesh->getPCU(), "non-aligned elements down from %d to %d in %f seconds", - originalCount,count,t1-t0); -} void printQuality(Adapt* a) { if ( ! a->input->shouldPrintQuality) return; - double minqual = getMinQuality(a); + double minqual = cbrt(getMinQuality(a)); print(a->mesh->getPCU(), "worst element quality is %e", minqual); } diff --git a/ma/maShape.h b/ma/maShape.h index d72c7bbaa..725fbd0a1 100644 --- a/ma/maShape.h +++ b/ma/maShape.h @@ -20,6 +20,7 @@ class Adapt; /* quick check of positivity of volumes based on vertices */ bool areTetsValid(Mesh* m, EntityArray& tets); +bool isTetValid(Mesh* m, Entity* tets); double measureTriQuality(Mesh* m, SizeField* f, Entity* tri, bool useMax=true); double measureTetQuality(Mesh* m, SizeField* f, Entity* tet, bool useMax=true); @@ -33,6 +34,7 @@ double measureQuadraticTetQuality(Mesh* m, Entity* tet); double getWorstQuality(Adapt* a, EntityArray& e); double getWorstQuality(Adapt* a, Entity** e, size_t n); +double getAndCacheQuality(Adapt* a, Entity* e); /* has worse quality than qualityToBeat */ @@ -63,13 +65,7 @@ bool isPyramidOk(apf::Mesh* m, Entity* e, int* good_rotation = 0); bool isLayerElementOk(Mesh* m, Entity* e); -CodeMatch matchSliver( - Mesh* m, - Entity* tet); - double improveQualities(Adapt* a); -void fixElementShapes(Adapt* a); -void alignElements(Adapt* a); void printQuality(Adapt* a); } diff --git a/ma/maSize.cc b/ma/maSize.cc index e596d0788..e3a43f1fb 100644 --- a/ma/maSize.cc +++ b/ma/maSize.cc @@ -13,6 +13,8 @@ #include #include #include +#include "apfVectorElement.h" +#include "maAdapt.h" namespace ma { @@ -53,9 +55,8 @@ IdentitySizeField::IdentitySizeField(Mesh* m): double IdentitySizeField::measure(Entity* e) { - apf::MeshElement* me = apf::createMeshElement(mesh,e); - double x = apf::measure(me); - apf::destroyMeshElement(me); + apf::MeshElement me(mesh->getCoordinateField(),e); + double x = apf::measure(&me); return x; } @@ -69,6 +70,18 @@ bool IdentitySizeField::shouldCollapse(Entity*) return false; } +SizeField::Value IdentitySizeField::getValue(Entity*) +{ + Value value; + value.matrix = Matrix(1,0,0, + 0,1,0, + 0,0,1); + return value; +} +void IdentitySizeField::setValue(Entity*, Value const&) +{ +} + void IdentitySizeField::interpolate( apf::MeshElement*, Vector const&, @@ -184,6 +197,7 @@ class SizeFieldIntegrator : public apf::Integrator { meshElement = me; dimension = apf::getDimension(me); + measurement = 0; } void atPoint(Vector const& p , double w, double ) { @@ -205,14 +219,29 @@ class SizeFieldIntegrator : public apf::Integrator struct MetricSizeField : public SizeField { + MetricSizeField(Mesh* m, int o) : mesh(m), order(o) + { + } double measure(Entity* e) { - SizeFieldIntegrator sFI(this, - std::max(mesh->getShape()->getOrder(), order)+1); - apf::MeshElement* me = apf::createMeshElement(mesh, e); - sFI.process(me); - apf::destroyMeshElement(me); - return sFI.measurement; + me.init(mesh->getCoordinateField(), e, 0); + int integrationOrder = std::max(mesh->getShape()->getOrder(), order)+1; + double measurement = 0; + int dim = apf::getDimension(&me); + int np = get_or_add(numIntegrationPoints, dim, [&]() {return countIntPoints(&me,integrationOrder);}); + double w = get_or_add(integrationWeight, dim, [&]() {return getIntWeight(&me,integrationOrder,0);}); + for (int p=0; p < np; ++p) { + Vector point = get_or_add(integrationPoint, std::make_pair(dim, p), [&]() {return getIntPoint(&me,integrationOrder,p);}); + Matrix Q; + getTransform(&me,point,Q); + Matrix J; + apf::getJacobian(&me,point,J); + /* transforms the rows of J, the differential tangent vectors, + into the metric space, then uses the generalized determinant */ + double dV2 = apf::getJacobianDeterminant(J*Q,dim); + measurement += w*dV2; + } + return measurement; } bool shouldSplit(Entity* edge) { @@ -229,6 +258,10 @@ struct MetricSizeField : public SizeField } Mesh* mesh; int order; // this is the underlying sizefield order (default 1) + apf::MeshElement me; + std::map numIntegrationPoints; + std::map integrationWeight; + std::map, Vector> integrationPoint; }; AnisotropicFunction::~AnisotropicFunction() @@ -363,18 +396,14 @@ struct LogMEval : public apf::Function struct AnisoSizeField : public MetricSizeField { - AnisoSizeField() + AnisoSizeField() : MetricSizeField(0, 1) { - mesh = 0; - order = 1; } - AnisoSizeField(Mesh* m, AnisotropicFunction* f): + AnisoSizeField(Mesh* m, AnisotropicFunction* f): MetricSizeField(m, 1), bothEval(f), sizesEval(&bothEval), frameEval(&bothEval) { - mesh = m; - order = 1; hField = apf::createUserField(m, "ma_sizes", apf::VECTOR, apf::getLagrange(1), &sizesEval); rField = apf::createUserField(m, "ma_frame", apf::MATRIX, @@ -397,35 +426,48 @@ struct AnisoSizeField : public MetricSizeField Vector const& xi, Matrix& Q) { - apf::Element* hElement = apf::createElement(hField,me); - apf::Element* rElement = apf::createElement(rField,me); + if (me->getEntity() != hElement.getEntity()) + hElement.init(hField,me->getEntity(),me); + if (me->getEntity() != rElement.getEntity()) + rElement.init(rField,me->getEntity(),me); Vector h; Matrix R; - apf::getVector(hElement,xi,h); - apf::getMatrix(rElement,xi,R); - apf::destroyElement(hElement); - apf::destroyElement(rElement); + apf::getVector(&hElement,xi,h); + apf::getMatrix(&rElement,xi,R); orthogonalizeR(R); Matrix S(1/h[0],0,0, 0,1/h[1],0, 0,0,1/h[2]); Q = R*S; } + Value getValue(Entity* vert) + { + apf::MeshElement me(mesh->getCoordinateField(), vert); + if (vert != hElement.getEntity()) hElement.init(hField,vert,&me); + if (vert != rElement.getEntity()) rElement.init(rField,vert,&me); + Value value; + apf::getVector(&hElement,Vector(0.0, 0.0, 0.0),value.vector); + apf::getMatrix(&rElement,Vector(0.0, 0.0, 0.0),value.matrix); + return value; + } + void setValue(Entity* vert, Value const& value) + { + apf::setMatrix(rField,vert,0,value.matrix); + apf::setVector(hField,vert,0,value.vector); + } void interpolate( apf::MeshElement* parent, Vector const& xi, Entity* newVert) { - apf::Element* rElement = apf::createElement(rField,parent); - apf::Element* hElement = apf::createElement(hField,parent); + if (parent->getEntity() != hElement.getEntity()) hElement.init(hField,parent->getEntity(),parent); + if (parent->getEntity() != rElement.getEntity()) rElement.init(rField,parent->getEntity(),parent); Vector h; - apf::getVector(hElement,xi,h); + apf::getVector(&hElement,xi,h); Matrix R; - apf::getMatrix(rElement,xi,R); + apf::getMatrix(&rElement,xi,R); orthogonalizeR(R); this->setValue(newVert,R,h); - apf::destroyElement(hElement); - apf::destroyElement(rElement); } void setValue( Entity* vert, @@ -445,6 +487,8 @@ struct AnisoSizeField : public MetricSizeField 0,0,1), Vector(value,value,value)); } + apf::Element hElement; + apf::Element rElement; apf::Field* hField; apf::Field* rField; BothEval bothEval; @@ -454,12 +498,10 @@ struct AnisoSizeField : public MetricSizeField struct LogAnisoSizeField : public MetricSizeField { - LogAnisoSizeField() + LogAnisoSizeField() : MetricSizeField(0, 1) { - mesh = 0; - order = 1; } - LogAnisoSizeField(Mesh* m, AnisotropicFunction* f): + LogAnisoSizeField(Mesh* m, AnisotropicFunction* f): MetricSizeField(m, 1), logMEval(f) { mesh = m; @@ -485,19 +527,19 @@ struct LogAnisoSizeField : public MetricSizeField continue; it = m->begin(d); while( (ent = m->iterate(it)) ){ - int type = m->getType(ent); - int non = apf::getShape(logMField)->countNodesOn(type); - for (int i = 0; i < non; i++) { - Vector h; - Matrix f; - apf::getVector(sizes, ent, i, h); - apf::getMatrix(frames, ent, i, f); - Vector s(log(1/h[0]/h[0]), log(1/h[1]/h[1]), log(1/h[2]/h[2])); - Matrix S(s[0], 0 , 0, - 0 , s[1], 0, - 0 , 0 , s[2]); - apf::setMatrix(logMField, ent, i, f * S * transpose(f)); - } + int type = m->getType(ent); + int non = apf::getShape(logMField)->countNodesOn(type); + for (int i = 0; i < non; i++) { + Vector h; + Matrix f; + apf::getVector(sizes, ent, i, h); + apf::getMatrix(frames, ent, i, f); + Vector s(log(1/h[0]/h[0]), log(1/h[1]/h[1]), log(1/h[2]/h[2])); + Matrix S(s[0], 0 , 0, + 0 , s[1], 0, + 0 , 0 , s[2]); + apf::setMatrix(logMField, ent, i, f * S * transpose(f)); + } } m->end(it); } @@ -508,10 +550,9 @@ struct LogAnisoSizeField : public MetricSizeField Vector const& xi, Matrix& Q) { - apf::Element* logMElement = apf::createElement(logMField,me); + if (me->getEntity() != logMElement.getEntity()) logMElement.init(logMField,me->getEntity(),me); Matrix logM; - apf::getMatrix(logMElement,xi,logM); - apf::destroyElement(logMElement); + apf::getMatrix(&logMElement,xi,logM); Vector v; Matrix R; orthogonalEigenDecompForSymmetricMatrix(logM, v, R); @@ -520,16 +561,28 @@ struct LogAnisoSizeField : public MetricSizeField 0, 0, sqrt(exp(v[2]))); Q = R*S; } + Value getValue(Entity* vert) + { + apf::MeshElement me(mesh->getCoordinateField(), vert); + if (vert != logMElement.getEntity()) logMElement.init(logMField,vert,&me); + Value value; + apf::getMatrix(&logMElement,Vector(0.0, 0.0, 0.0),value.matrix); + return value; + } + void setValue(Entity* vert, Value const& value) + { + apf::setMatrix(logMField,vert,0,value.matrix); + } void interpolate( apf::MeshElement* parent, Vector const& xi, Entity* newVert) { - apf::Element* logMElement = apf::createElement(logMField,parent); + if (parent->getEntity() != logMElement.getEntity()) + logMElement.init(logMField,parent->getEntity(),parent); Matrix logM; - apf::getMatrix(logMElement,xi,logM); + apf::getMatrix(&logMElement,xi,logM); this->setValue(newVert,logM); - apf::destroyElement(logMElement); } void setValue( Entity* vert, @@ -574,6 +627,7 @@ struct LogAnisoSizeField : public MetricSizeField return apf::getShape(logMField)->hasNodesIn(dimension); } apf::NewArray fieldVal; + apf::Element logMElement; apf::Field* logMField; LogMEval logMEval; }; @@ -615,10 +669,42 @@ struct IsoUserField : public IsoSizeField FieldReader reader; }; +static void clampSizeField(Mesh*m, apf::Field* sizes) +{ + Vector lower; + Vector upper; + getBoundingBox(m, lower, upper); + double max = std::abs(lower[0]) + std::abs(upper[0]); + max = std::max(max, std::abs(lower[1]) + std::abs(upper[1])); + max = std::max(max, std::abs(lower[2]) + std::abs(upper[2])); + max = max/2; + bool clamped = false; + + Entity* ent; + Iterator* it; + for (int d = 0; d <= m->getDimension(); d++) { + it = m->begin(d); + while( (ent = m->iterate(it)) ){ + int type = m->getType(ent); + int non = sizes->getShape()->countNodesOn(type); + for (int i = 0; i < non; i++) { + Vector h; + apf::getVector(sizes, ent, i, h); + if (h[0] > max) {clamped = true; h[0] = max;} + if (h[1] > max) {clamped = true; h[1] = max;} + if (h[2] > max) {clamped = true; h[2] = max;} + apf::setVector(sizes, ent, i, h); + } + } + } + if (clamped) print(m->getPCU(), "[WARNING]: Found size field larger than the bounding box. " + "Size field has been automatically clamped to the bounding box."); +} + SizeField* makeSizeField(Mesh* m, apf::Field* sizes, apf::Field* frames, bool logInterpolation) { - // logInterpolation is "false" by default + clampSizeField(m, sizes); if (! logInterpolation) { AnisoSizeField* anisoF = new AnisoSizeField(); anisoF->init(m, sizes, frames); @@ -651,6 +737,15 @@ SizeField* makeSizeField(Mesh* m, apf::Field* size) return new IsoUserField(m, size); } +double getAndCacheSize(Adapt* a, Entity* e) +{ + if (a->mesh->hasTag(e, a->sizeCache)) + return getCachedSize(a, e); + double size = a->sizeField->measure(e); + setCachedSize(a, e, size); + return size; +} + double getAverageEdgeLength(Mesh* m) { IdentitySizeField sizeField(m); diff --git a/ma/maSize.h b/ma/maSize.h index e8f9b6058..c5b0d7977 100644 --- a/ma/maSize.h +++ b/ma/maSize.h @@ -18,6 +18,7 @@ namespace ma { +class Adapt; 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, @@ -26,14 +27,22 @@ 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; +const double GOODQUALITY = .3; class SizeField { public: + struct Value { + Matrix matrix; + Vector vector; + }; + virtual ~SizeField(); virtual double measure(Entity* e) = 0; virtual bool shouldSplit(Entity* edge) = 0; virtual bool shouldCollapse(Entity* edge) = 0; + virtual Value getValue(Entity* vert) = 0; + virtual void setValue(Entity* vert, Value const& value) = 0; virtual void interpolate( apf::MeshElement* parent, Vector const& xi, @@ -59,6 +68,8 @@ struct IdentitySizeField : public SizeField double measure(Entity* e); bool shouldSplit(Entity*); bool shouldCollapse(Entity*); + Value getValue(Entity* vert); + void setValue(Entity* vert, Value const& value); void interpolate( apf::MeshElement* parent, Vector const& xi, @@ -112,9 +123,20 @@ SizeField* makeSizeField(Mesh* m, AnisotropicFunction* f, SizeField* makeSizeField(Mesh* m, apf::Field* size); SizeField* makeSizeField(Mesh* m, IsotropicFunction* f); +double getAndCacheSize(Adapt* a, Entity* e); double getAverageEdgeLength(Mesh* m); double getMaximumEdgeLength(Mesh* m, SizeField* sf = 0); +template +typename Map::mapped_type& +get_or_add(Map& m, const Key& key, Func valueFn) +{ + auto it = m.find(key); + if (it == m.end()) + it = m.insert(std::make_pair(key, valueFn())).first; + return it->second; +} + } #endif diff --git a/ma/maSnap.cc b/ma/maSnap.cc index 3c4d92b76..b82753c65 100644 --- a/ma/maSnap.cc +++ b/ma/maSnap.cc @@ -12,6 +12,7 @@ #include "maOperator.h" #include "maSnapper.h" #include "maRefine.h" +#include "maRefine.h" #include "maMatchedSnapper.h" #include "maLayer.h" #include "maMatch.h" @@ -24,6 +25,12 @@ #include #include +#ifdef PUMI_HAS_CAPSTONE +#include "gmi_cap.h" +#endif +#include +#include + #ifdef PUMI_HAS_CAPSTONE #include "gmi_cap.h" #endif @@ -53,14 +60,14 @@ static size_t isSurfUnderlyingFaceDegenerate( Vector bmin; Vector bmax; - m->boundingBox(g, bmin, bmax); + m->boundingBoxCached(g, bmin, bmax); double tol = 1.0e-6 * (bmax - bmin).getLength(); - bool isPeriodic[2]; - double range[2][2]; + bool isPeriodic[2] = {0}; + double range[2][2] = {{0}}; int numPeriodicDims = 0; for (int i = 0; i < md; i++) { - isPeriodic[i] = m->getPeriodicRange(g,i,range[i]); + isPeriodic[i] = m->getPeriodicRangeCached(g,i,range[i]); if (isPeriodic[i]) numPeriodicDims++; if (range[i][0] > range[i][1]) @@ -163,7 +170,7 @@ static void interpolateParametricCoordinateOnEdge( Vector& p) { double range[2]; - bool isPeriodic = m->getPeriodicRange(g,0,range); + bool isPeriodic = m->getPeriodicRangeCached(g,0,range); p[0] = interpolateParametricCoordinate(t,a[0],b[0],range,isPeriodic, 0); p[1] = 0.0; p[2] = 0.0; @@ -444,8 +451,8 @@ static void interpolateParametricCoordinatesOnDegenerateFace( double d_range[2]; int p_axes = 1 - d_axes; - PCU_ALWAYS_ASSERT(m->getPeriodicRange(g, p_axes, p_range)); - PCU_ALWAYS_ASSERT(!m->getPeriodicRange(g, d_axes, d_range)); + PCU_ALWAYS_ASSERT(m->getPeriodicRangeCached(g, p_axes, p_range)); + PCU_ALWAYS_ASSERT(!m->getPeriodicRangeCached(g, d_axes, d_range)); if (numPoles == 2) { @@ -496,7 +503,7 @@ static void interpolateParametricCoordinatesOnRegularFace( int dim = m->getModelType(g); bool gface_isPeriodic = 0; for (int d=0; d < dim; ++d) { - bool isPeriodic = m->getPeriodicRange(g,d,range); + bool isPeriodic = m->getPeriodicRangeCached(g,d,range); if ((dim == 2) && (isPeriodic > 0)) gface_isPeriodic = 1; p[d] = interpolateParametricCoordinate(t,a[d],b[d],range,isPeriodic, 0); } @@ -523,7 +530,7 @@ static void interpolateParametricCoordinatesOnRegularFace( return; for (int d=0; d < dim; ++d) { - bool isPeriodic = m->getPeriodicRange(g,d,range); + bool isPeriodic = m->getPeriodicRangeCached(g,d,range); p[d] = interpolateParametricCoordinate(t,a[d],b[d],range,isPeriodic, 1); } } @@ -834,7 +841,6 @@ class SnapMatched : public Operator didAnything = didAnything || snapped; if (snapped) ++successCount; - clearFlagMatched(adapter, vert, SNAP); } int successCount; bool didAnything; @@ -853,26 +859,29 @@ bool snapMatchedVerts(Adapt* a, Tag* t, long& successCount) return a->mesh->getPCU()->Or(op.didAnything); } -long tagVertsToSnap(Adapt* a, Tag*& tag) +void tagVertexToSnap(Adapt* a, Entity* vertex) { + if (!a->input->shouldSnap) return; Mesh* m = a->mesh; int dim = m->getDimension(); + int md = m->getModelType(m->toModel(vertex)); + if (dim == 3 && md == 3) return; + Vector snapPoint; + getSnapPoint(m, vertex, snapPoint); + Vector position = getPosition(m, vertex); + setFlag(a, vertex, SNAP); + m->setDoubleTag(vertex, a->snapTag, &snapPoint[0]); +} + +long countTaggedVerts(Adapt* a) +{ + Mesh* m = a->mesh; Entity* v; long n = 0; Iterator* it = m->begin(0); while ((v = m->iterate(it))) { - int md = m->getModelType(m->toModel(v)); - if (dim == 3 && md == 3) - continue; - Vector s; - getSnapPoint(m, v, s); - Vector x = getPosition(m, v); - if (apf::areClose(s, x, 1e-12)) - continue; - m->setDoubleTag(v, tag, &s[0]); - setFlag(a, v, SNAP); - if (m->isOwned(v)) - ++n; + if (!m->hasTag(v, a->snapTag)) continue; + if (m->isOwned(v)) ++n; } m->end(it); return m->getPCU()->Add(n); @@ -887,7 +896,7 @@ long snapTaggedVerts(Adapt* a, Tag* tag, Snapper& snapper) bool snapped = true; int prevTagged = 0; while (snapped) { - int tagged = tagVertsToSnap(a, tag); + int tagged = countTaggedVerts(a); if (tagged == prevTagged) break; prevTagged = tagged; if (a->mesh->hasMatching()) @@ -915,18 +924,17 @@ void snap(Adapt* a) if (!a->input->shouldSnap) return; double t0 = pcu::Time(); - Tag* snapTag = a->mesh->createDoubleTag("ma_snap", 3); preventMatchedCavityMods(a); - int toSnap = tagVertsToSnap(a, snapTag); + int toSnap = countTaggedVerts(a); if (toSnap) { - Snapper snapper(a, snapTag); - snapTaggedVerts(a, snapTag, snapper); - snapLayer(a, snapTag); + Snapper snapper(a, a->snapTag); + snapTaggedVerts(a, a->snapTag, snapper); + snapLayer(a, a->snapTag); double t1 = pcu::Time(); - print(a->mesh->getPCU(), "ToSnap %d - Moved %d - Failed %d - CollapseToVtx %d" - " - Collapse %d - Swap %d - SplitCollapse %d" + print(a->mesh->getPCU(), "snapping %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), @@ -934,6 +942,5 @@ void snap(Adapt* a) if (a->hasLayer) checkLayerShape(a->mesh, "after snapping"); } - a->mesh->destroyTag(snapTag); } } diff --git a/ma/maSnap.h b/ma/maSnap.h index be110eed7..e91995828 100644 --- a/ma/maSnap.h +++ b/ma/maSnap.h @@ -19,7 +19,7 @@ class Adapt; void snap(Adapt* a); void visualizeGeometricInfo(Mesh* m, const char* name); -long snapTaggedVerts(Adapt* a, Tag* snapTag); +void tagVertexToSnap(Adapt* a, Entity* vertex); void interpolateParametricCoordinates( apf::Mesh* m, diff --git a/ma/maSnapper.cc b/ma/maSnapper.cc index 2f3a426d2..1cb7efb33 100644 --- a/ma/maSnapper.cc +++ b/ma/maSnapper.cc @@ -15,20 +15,32 @@ * it will collapse to simplify the region and attempt other operators such as * swap, split collapse, double split collapse. */ +/** + * \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 "maFaceSwap.h" #include "maSnap.h" #include "maDBG.h" +#include "maDBG.h" #include #include #include #include +#include "apfGeometry.cc" namespace ma { -Snapper::Snapper(Adapt* a, Tag* st) : mesh(a->mesh), splitCollapse(a), doubleSplitCollapse(a) +Snapper::Snapper(Adapt* a, Tag* st) : mesh(a->mesh), splitCollapse(a), doubleSplitCollapse(a), reposition(a) { + adapt = a; adapt = a; snapTag = st; collapse.Init(a); @@ -83,23 +95,25 @@ static void flagAndPrint(Adapt* a, Entity* ent, int dim, const char* name) #if defined(DEBUG_FPP) static void printFPP(Adapt* a, FirstProblemPlane* FPP) { - ma_dbg::addTargetLocation(a, "snap_target"); - ma_dbg::addClassification(a, "classification"); + ma_dbg::useFieldInfo(a, [a, FPP] { + apf::writeVtkFiles("FPP_Mesh", a->mesh); + EntitySet invalid; + for (auto e : FPP->problemRegions) invalid.insert(e); + ma_dbg::createCavityMesh(a, invalid, "FPP_Invalid"); - apf::writeVtkFiles("FPP_Mesh", a->mesh); - EntityArray invalid; - for (int i=0; iproblemRegions.n; i++){ - invalid.append(FPP->problemRegions.e[i]); - } - ma_dbg::createCavityMesh(a, invalid, "FPP_Invalid"); + for (auto e : FPP->commEdges) setFlag(a, e, CHECKED); + ma_dbg::dumpMeshWithFlag(a, 0, 1, CHECKED, "FPP_CommEdges", "FPP_CommEdges"); + for (auto e : FPP->commEdges) clearFlag(a, e, CHECKED); - 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"); - flagAndPrint(a, FPP->vert, 0, "FPP_Vertex"); - flagAndPrint(a, FPP->problemFace, 2, "FPP_Face"); - flagAndPrint(a, FPP->problemRegion, 3, "FPP_Region"); + EntitySet adjacent1 = getNextLayer(a, invalid); + ma_dbg::createCavityMesh(a, adjacent1, "FPP_ADJACENT_1"); + EntitySet adjacent2 = getNextLayer(a, adjacent1); + ma_dbg::createCavityMesh(a, adjacent2, "FPP_ADJACENT_2"); + }); } #endif @@ -300,6 +314,7 @@ 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 qual = adapt->input->validQuality; double min=area[0]; for(int i=1; i<4; i++) @@ -313,14 +328,8 @@ bool Snapper::trySwapOrSplit(FirstProblemPlane* FPP) 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; - } + if (edgeSwap->run(longest)) { numSwap++; return true; } + if (splitCollapse.run(longest, FPP->vert, qual)) { numSplitCollapse++; return true; } } if (ents[0] == 0) @@ -328,38 +337,19 @@ bool Snapper::trySwapOrSplit(FirstProblemPlane* FPP) // 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"); + if (edgeSwap->run(ents[0])) { numSwap++; return true; } + if (edgeSwap->run(ents[1])) { numSwap++; return true; } + if (splitCollapse.run(ents[0], FPP->vert, qual)) { numSplitCollapse++; return true; } + if (splitCollapse.run(ents[1], FPP->vert, qual)) { numSplitCollapse++; return true; } + if (doubleSplitCollapse.run(ents, qual)) { numSplitCollapse++; return true; } } // 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"); + for (int i=0; i<3; i++) + if (edgeSwap->run(edges[i])) { numSwap++; return true; } + if (splitCollapse.run(ents[1], FPP->vert, qual)) { numSplitCollapse++; return true; } } return false; } @@ -560,43 +550,12 @@ static FirstProblemPlane* getFPP(Adapt* a, Entity* vertex, Tag* snapTag, apf::Up 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 invalid; - return tryReposition(adapt, vert, snapTag, invalid); + if (!mesh->hasTag(vert, snapTag)) return true; + Vector target; + adapt->mesh->getDoubleTag(vert, snapTag, &target[0]); + return reposition.move(vert, target); } #if defined(DEBUG_FPP) static int DEBUGFAILED=0; @@ -610,8 +569,12 @@ to apply certain opperators so other algoritms were adapted from old scorec libr */ bool Snapper::run() { - apf::Up invalid; - bool success = tryReposition(adapt, vert, snapTag, invalid); + if (!mesh->hasTag(vert, snapTag)) return true; + Vector target; + adapt->mesh->getDoubleTag(vert, snapTag, &target[0]); + bool success = reposition.move(vert, target); + apf::Up& invalid = reposition.getInvalid(); + if (success) { numSnapped++; mesh->removeTag(vert,snapTag); @@ -628,11 +591,7 @@ bool Snapper::run() if (!success) success = trySwapOrSplit(FPP); } - if (!success) { - numFailed++; - mesh->removeTag(vert,snapTag); - clearFlag(adapt, vert, SNAP); - } + if (!success) numFailed++; #if defined(DEBUG_FPP) if (!success && ++DEBUGFAILED == 1) printFPP(adapt, FPP); #endif @@ -951,4 +910,20 @@ bool isLowInHigh(Mesh* mesh, Entity* highEnt, Entity* lowEnt) return false; } +EntitySet getNextLayer(Adapt* a, EntitySet& tets) +{ + EntitySet adjacent; + APF_ITERATE(ma::EntitySet,tets,it) { + Entity* faces[4]; + a->mesh->getDownward(*it, 2, faces); + for (int f=0; f<4; f++) { + apf::Up nextLayer; + a->mesh->getUp(faces[f], nextLayer); + for (int n=0; n& coords); Vector getCenter(Mesh* mesh, Entity* face); bool isLowInHigh(Mesh* mesh, Entity* highEnt, Entity* lowEnt); +EntitySet getNextLayer(Adapt* a, EntitySet& tets); } diff --git a/ma/maSolutionTransferHelper.cc b/ma/maSolutionTransferHelper.cc index d799d9b1b..750dca266 100644 --- a/ma/maSolutionTransferHelper.cc +++ b/ma/maSolutionTransferHelper.cc @@ -9,6 +9,7 @@ opyright 2013 Scientific Computation Research Center, *******************************************************************************/ #include "maSolutionTransferHelper.h" +#include "apfVectorElement.h" #include "maAffine.h" #include "maMap.h" #include @@ -76,9 +77,8 @@ void transferToNode( point = childMap * xi; } else { // else inquire the physical coordinated of local coordinate xi - apf::MeshElement* me = apf::createMeshElement(mesh,node.entity); - apf::mapLocalToGlobal(me, xi, point); - apf::destroyMeshElement(me); + apf::MeshElement me(mesh->getCoordinateField(),node.entity); + apf::mapLocalToGlobal(&me, xi, point); } Vector elemXi; int i = getBestElement(mesh,n,elems,elemInvMaps,point,elemXi); diff --git a/ma/maSplits.h b/ma/maSplits.h index 22b4e1aab..36e5ab70f 100644 --- a/ma/maSplits.h +++ b/ma/maSplits.h @@ -26,8 +26,8 @@ class Splits Entity* getSplitVert(int i); EntityArray& getTets() {return refiner->toSplit[3];} Adapt* getAdapt() {return refiner->adapt;} - private: Refine* refiner; + private: }; } diff --git a/ma/maStats.cc b/ma/maStats.cc index 479dcce21..2db4e18b4 100644 --- a/ma/maStats.cc +++ b/ma/maStats.cc @@ -6,6 +6,9 @@ */ #include "maStats.h" #include "maAdapt.h" +#include +#include +#include namespace ma { @@ -134,4 +137,69 @@ void stats(ma::Mesh* m, ma::SizeField* sf, getStatsInPhysicalSpace(m, edgeLengths, linearQualities); } +std::vector printHistogramData(std::string name, std::vector bins, std::vector input, double min, double max, Mesh* m) +{ + std::vector count(bins.size()-1, 0); + double inputMax = min; + double inputMin = max; + + for (size_t i = 0; i < input.size(); ++i) { + if (std::isnan(input[i])) continue; + if (input[i] > inputMax) inputMax = input[i]; + if (input[i] < inputMin) inputMin = input[i]; + + auto it = std::upper_bound(bins.begin(), bins.end(), input[i]); + size_t binIdx = std::distance(bins.begin(), it) - 1; + if (binIdx >= count.size()) binIdx = count.size()-1; + count[binIdx]++; + } + + inputMin = m->getPCU()->Min(inputMin); + inputMax = m->getPCU()->Max(inputMax); + for (size_t i = 0; i < count.size(); ++i) count[i] = m->getPCU()->Add(count[i]); + + if (m->getPCU()->Self()) return count; + printf("%s Min: %f, Max: %f\n", name.c_str(), inputMin, inputMax); + for (size_t i = 0; i < count.size(); ++i) printf("%d\n", count[i]); + return count; +} + +HistogramStats printHistogramStats(Adapt* a) +{ + std::vector lengths; + std::vector qualities; + ma::stats(a->mesh, a->input->sizeField, lengths, qualities, true); + std::vector qualityBins = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0}; + std::vector lengthBins; + for (double i = 0.0; i <= 10.0; i++) lengthBins.push_back(i*((MAXLENGTH+1)/10)); + std::vector qualityHist = printHistogramData("\nQualities:", qualityBins, qualities, 0, 1, a->mesh); + std::vector lengthHist = printHistogramData("\nLengths:", lengthBins, lengths, 0, MAXLENGTH+1, a->mesh); + return HistogramStats(qualityHist, lengthHist); +} + +//Compare two histograms using L2 (Euclidean) distance. +double histogramDistance(const std::vector& hist1, const std::vector& hist2, bool normalize) +{ + if (hist1.size() != hist2.size()) throw std::invalid_argument("Histograms must be the same size"); + + size_t n = hist1.size(); + std::vector h1(hist1.begin(), hist1.end()); + std::vector h2(hist2.begin(), hist2.end()); + + if (normalize) { + double sum1 = std::accumulate(h1.begin(), h1.end(), 0.0); + double sum2 = std::accumulate(h2.begin(), h2.end(), 0.0); + if (sum1 > 0) for (double& val : h1) val /= sum1; + if (sum2 > 0) for (double& val : h2) val /= sum2; + } + + double sum_sq_diff = 0.0; + for (size_t i = 0; i < n; ++i) { + double diff = h1[i] - h2[i]; + sum_sq_diff += diff * diff; + } + return std::sqrt(sum_sq_diff); +} + + } diff --git a/ma/maStats.h b/ma/maStats.h index e577d6a15..2f42e0f3a 100644 --- a/ma/maStats.h +++ b/ma/maStats.h @@ -44,5 +44,18 @@ void stats(ma::Mesh* m, ma::SizeField* sf, std::vector &linearQualities, bool inMetric); +std::vector printHistogramData(std::string name, std::vector input, + double min, double max, Mesh* m); + +struct HistogramStats { + std::vector quality; + std::vector length; + HistogramStats(std::vector q, std::vector l) : quality(q), length(l) {} +}; + +HistogramStats printHistogramStats(Adapt* a); + +double histogramDistance(const std::vector& hist1, const std::vector& hist2, bool normalize=true); + } #endif diff --git a/ma/pkg_tribits.cmake b/ma/pkg_tribits.cmake index d1db5c276..62d340dbd 100644 --- a/ma/pkg_tribits.cmake +++ b/ma/pkg_tribits.cmake @@ -26,6 +26,8 @@ set(SOURCES maSnap.cc maEdgeSwap.cc maShape.cc + maFixShape.cc + maFaceSwap.cc maShapeHandler.cc maQuality.cc maSplits.cc diff --git a/test/aniso_adapt.cc b/test/aniso_adapt.cc index beb582a23..2795e4f81 100644 --- a/test/aniso_adapt.cc +++ b/test/aniso_adapt.cc @@ -13,13 +13,10 @@ int main(int argc, char** argv) pcu::PCU PCUObj; lion_set_verbosity(1); gmi_register_mesh(); - ma::Mesh* mesh1 = apf::loadMdsMesh(modelFile,meshFile,&PCUObj); - ma::Mesh* mesh2 = apf::loadMdsMesh(modelFile,meshFile,&PCUObj); - adaptTests(mesh1, mesh2); - mesh1->destroyNative(); - mesh2->destroyNative(); - apf::destroyMesh(mesh1); - apf::destroyMesh(mesh2); + ma::Mesh* mesh = apf::loadMdsMesh(modelFile,meshFile,&PCUObj); + adaptTests(mesh, {0, 0, 0, 5168, 10638, 18363, 25244, 28600, 21121, 6330}); + mesh->destroyNative(); + apf::destroyMesh(mesh); } pcu::Finalize(); } diff --git a/test/aniso_adapt.h b/test/aniso_adapt.h index 160a88165..b1e3794af 100644 --- a/test/aniso_adapt.h +++ b/test/aniso_adapt.h @@ -15,6 +15,8 @@ #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 @@ -45,20 +47,20 @@ class AnIso : public ma::AnisotropicFunction ma::Vector lower, upper; }; -void measureQuality(ma::Mesh* m, double& avgQuality, double& minQuality) +void measureQuality(ma::Adapt* a, double& avgQuality, double& minQuality) { - int d = m->getDimension(); + ma::Mesh* m = a->mesh; apf::MeshEntity* elem; - apf::MeshIterator* it = m->begin(d); - ma::IdentitySizeField I(m); + apf::MeshIterator* it = m->begin(m->getDimension()); minQuality = 1; + avgQuality = 0; while ((elem = m->iterate(it))) { - double q = ma::measureElementQuality(m, &I, elem); + double q = cbrt(ma::getAndCacheQuality(a, elem)); avgQuality += q; minQuality = std::min(minQuality, q); } m->end(it); - avgQuality = avgQuality / m->count(d); + avgQuality = avgQuality / m->count(m->getDimension()); } //make sure refinement is done @@ -78,54 +80,59 @@ int countEdges(ma::Mesh* m) return m->count(1); } -ma::Mesh* coarsenForced(ma::Mesh* m) +ma::Mesh* fixShapeTest(ma::Mesh* m, std::vector savedQuality) { m->verify(); - AnIso sf(m, .5, 1); + AnIso sf(m, 2, 2); ma::Input* in = ma::makeAdvanced(ma::configure(m, &sf)); - in->shouldForceAdaptation = true; + in->shouldPrintQuality = false; ma::Adapt* a = new ma::Adapt(in); - double avgQualBefore=0, avgQualAfter=0, minQualBefore=0, minQualAfter=0; - measureQuality(m, avgQualBefore, minQualBefore); - double averageBefore = ma::getAverageEdgeLength(m); - int edgesBefore = countEdges(m); + for (int i = 0; i < in->maximumIterations; ++i) + { + ma::refine(a); + ma::snap(a); + ma::coarsen(a); + } - ma::coarsen(a); + double avgQualBefore=0, avgQualAfter=0, minQualBefore=0, minQualAfter=0; + measureQuality(a, avgQualBefore, minQualBefore); + ma::fixElementShapes(a); + measureQuality(a, avgQualAfter, minQualAfter); - measureQuality(m, avgQualAfter, minQualAfter); + PCU_ALWAYS_ASSERT(minQualAfter > minQualBefore); + PCU_ALWAYS_ASSERT(avgQualAfter > avgQualBefore); - //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); + ma::HistogramStats hist = ma::printHistogramStats(a); + std::vector badQualityHist(hist.quality.begin(), hist.quality.begin() + 3); + std::vector badQualityHistSaved(savedQuality.begin(), savedQuality.begin() + 3); + PCU_ALWAYS_ASSERT(ma::histogramDistance(hist.quality, savedQuality) <= .05); + PCU_ALWAYS_ASSERT(ma::histogramDistance(badQualityHist, badQualityHistSaved) <= .05); 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; + if (in->ownsSizeField) delete in->sizeField; + if (in->ownsSolutionTransfer) delete in->solutionTransfer; delete in; return m; } -ma::Mesh* coarsenRegular(ma::Mesh* m) +ma::Mesh* coarsenTest(ma::Mesh* m) { m->verify(); AnIso sf(m, .5, 1); ma::Input* in = ma::makeAdvanced(ma::configure(m, &sf)); + in->shouldPrintQuality = false; ma::Adapt* a = new ma::Adapt(in); double avgQualBefore=0, avgQualAfter=0, minQualBefore=0, minQualAfter=0; - measureQuality(m, avgQualBefore, minQualBefore); + measureQuality(a, avgQualBefore, minQualBefore); double averageBefore = ma::getAverageEdgeLength(m); int edgesBefore = countEdges(m); ma::coarsen(a); - measureQuality(m, avgQualAfter, minQualAfter); + measureQuality(a, avgQualAfter, minQualAfter); //Makes sure that coarsening is happening and it isn't decreasing the quality of the mesh PCU_ALWAYS_ASSERT(edgesBefore > countEdges(m)); @@ -174,6 +181,7 @@ ma::Mesh* refineSnapTest(ma::Mesh* m) m->verify(); AnIso sf(m, 3, 1); ma::Input* in = ma::makeAdvanced(ma::configure(m, &sf)); + in->shouldPrintQuality = false; ma::Adapt* a = new ma::Adapt(in); int edgesBefore = countEdges(m); double averageBefore = ma::getAverageEdgeLength(m); @@ -185,7 +193,7 @@ ma::Mesh* refineSnapTest(ma::Mesh* m) } double avgQualAfter=0, minQualAfter=0; - measureQuality(m, avgQualAfter, minQualAfter); + measureQuality(a, avgQualAfter, minQualAfter); //Makes sure that refinement is happening PCU_ALWAYS_ASSERT(edgesBefore < countEdges(m)); @@ -206,22 +214,21 @@ ma::Mesh* refineSnapTest(ma::Mesh* m) return m; } -void adaptTests(ma::Mesh* meshReg, ma::Mesh* meshForced) +void adaptTests(ma::Mesh* meshReg, std::vector savedQuality) { apf::writeVtkFiles("startMesh", meshReg); + printf("\n==REFINE TEST==\n"); refineSnapTest(meshReg); apf::writeVtkFiles("afterRefine", meshReg); - coarsenRegular(meshReg); + printf("\n==COARSEN TEST==\n"); + coarsenTest(meshReg); apf::writeVtkFiles("afterCoarsen", meshReg); - coarsenForced(refineSnapTest(meshForced)); - apf::writeVtkFiles("afterForcedCoarsen", meshForced); - - //Make sure setting to force coarsen is functioning - PCU_ALWAYS_ASSERT(countEdges(meshReg) > countEdges(meshForced)); - PCU_ALWAYS_ASSERT(ma::getAverageEdgeLength(meshReg) < ma::getAverageEdgeLength(meshForced)); + printf("\n==FIX SHAPE TEST==\n"); + fixShapeTest(meshReg, savedQuality); + apf::writeVtkFiles("afterFixShape", meshReg); } #endif \ No newline at end of file diff --git a/test/aniso_adapt_cap.cc b/test/aniso_adapt_cap.cc index 464c9ed89..3063d0231 100644 --- a/test/aniso_adapt_cap.cc +++ b/test/aniso_adapt_cap.cc @@ -25,18 +25,14 @@ int main(int argc, char** argv) gmi_model* model = gmi_load(meshFile); ma::Mesh* apfCapMesh = apf::createCapMesh(model, &PCUObj); apf::disownCapModel(apfCapMesh); - ma::Mesh* mesh1 = apf::createMdsMesh(model, apfCapMesh); - apf::disownMdsModel(mesh1); - ma::Mesh* mesh2 = apf::createMdsMesh(model, apfCapMesh); - apf::disownMdsModel(mesh2); + ma::Mesh* mesh = apf::createMdsMesh(model, apfCapMesh); + apf::disownMdsModel(mesh); - adaptTests(apfCapMesh, mesh2); + adaptTests(apfCapMesh, {0, 0, 1, 2448, 5307, 9044, 13173, 14190, 10931, 3499}); - mesh1->destroyNative(); - mesh2->destroyNative(); + mesh->destroyNative(); apfCapMesh->destroyNative(); - apf::destroyMesh(mesh1); - apf::destroyMesh(mesh2); + apf::destroyMesh(mesh); apf::destroyMesh(apfCapMesh); gmi_cap_stop(); } diff --git a/test/aniso_adapt_sim.cc b/test/aniso_adapt_sim.cc index 74e472b1e..5991c2926 100644 --- a/test/aniso_adapt_sim.cc +++ b/test/aniso_adapt_sim.cc @@ -37,18 +37,14 @@ int main(int argc, char* argv[]) pGModel model = gmi_export_sim(mdl_ref); pParMesh mesh = PM_load(smsfile, model, progress); ma::Mesh* mesh_ref = apf::createMesh(mesh, PCUObj); - ma::Mesh* mesh1 = apf::createMdsMesh(mdl_ref, mesh_ref); - ma::Mesh* mesh2 = apf::createMdsMesh(mdl_ref, mesh_ref); - apf::disownMdsModel(mesh1); - apf::disownMdsModel(mesh2); + ma::Mesh* meshMDS = apf::createMdsMesh(mdl_ref, mesh_ref); + apf::disownMdsModel(meshMDS); - adaptTests(mesh1, mesh2); + adaptTests(meshMDS, {0, 0, 1, 89, 141, 205, 157, 169, 97, 29}); - mesh1->destroyNative(); - mesh2->destroyNative(); + meshMDS->destroyNative(); mesh_ref->destroyNative(); - apf::destroyMesh(mesh1); - apf::destroyMesh(mesh2); + apf::destroyMesh(meshMDS); apf::destroyMesh(mesh_ref); gmi_destroy(mdl_ref); Progress_delete(progress); diff --git a/test/testing.cmake b/test/testing.cmake index 67f75b948..8dbe931be 100644 --- a/test/testing.cmake +++ b/test/testing.cmake @@ -398,7 +398,7 @@ if(ENABLE_ZOLTAN) "${MESHES}/2dlayersNoAdapt/mesh.smb" "doNotAdapt" "8") set_tests_properties(ma_2dLayersOff PROPERTIES - PASS_REGULAR_EXPRESSION "number of triangle 18698") + PASS_REGULAR_EXPRESSION "number of triangle 18696") endif() mpi_test(tet_serial 1 ./tetrahedronize