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..3236a3909 100644 --- a/apf/apfCavityOp.cc +++ b/apf/apfCavityOp.cc @@ -119,6 +119,40 @@ void CavityOp::applyToDimension(int d) sharing = 0; } +void CavityOp::applyToList(std::list& elements) +{ + do { + delete sharing; + sharing = apf::getSharing(mesh); + + isRequesting = false; + for (auto iter = elements.begin(); iter != elements.end();) { + if (sharing->isOwned(*iter)) { + if (setEntity(*iter) == OK) { + movedByDeletion=false; + apply(); + if (movedByDeletion) {iter = elements.erase(iter); 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 (auto 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..f9934a80a 100644 --- a/apf/apfCavityOp.h +++ b/apf/apfCavityOp.h @@ -14,6 +14,8 @@ #include "apfMesh.h" #include #include +#include +#include namespace apf { @@ -87,12 +89,15 @@ 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 */ + void applyToList(std::list& 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 */ void preDeletion(MeshEntity* e); /** \brief mesh pointer for convenience */ Mesh* mesh; + bool movedByDeletion; private: typedef std::vector Requests; Requests requests; @@ -102,8 +107,8 @@ class CavityOp bool tryToPull(); void applyLocallyWithModification(int d); void applyLocallyWithoutModification(int d); + bool canModify; - bool movedByDeletion; MeshIterator* iterator; protected: Sharing* sharing; 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..0be5f16db 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,8 +37,8 @@ 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; 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/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..63cdc0d6b 100644 --- a/apf/apfMesh.cc +++ b/apf/apfMesh.cc @@ -15,6 +15,7 @@ #include #include #include +#include namespace apf { @@ -221,11 +222,23 @@ void Mesh::getParamOn(ModelEntity* g, MeshEntity* e, Vector3& p) gmi_reparam(getModel(), from, &from_p[0], to, &p[0]); } +std::map, + std::pair, bool>> periodicRanges; bool Mesh::getPeriodicRange(ModelEntity* g, int axis, double range[2]) { gmi_ent* e = (gmi_ent*)g; - gmi_range(getModel(), e, axis, range); - return gmi_periodic(getModel(), e, axis); + auto it = periodicRanges.find({e, axis}); + if (it != periodicRanges.end()) { + range[0] = it->second.first[0]; + range[1] = it->second.first[1]; + return it->second.second; + } + else{ + gmi_range(getModel(), e, axis, range); + bool output = gmi_periodic(getModel(), e, axis); + periodicRanges[{e, axis}] = {{range[0], range[1]}, output}; + return output; + } } void Mesh::getClosestPoint(ModelEntity* g, Vector3 const& from, diff --git a/apf/apfVectorElement.cc b/apf/apfVectorElement.cc index fa49e6a16..c5ffd52b1 100644 --- a/apf/apfVectorElement.cc +++ b/apf/apfVectorElement.cc @@ -10,6 +10,11 @@ namespace apf { +VectorElement::VectorElement(): + ElementOf() +{ +} + VectorElement::VectorElement(VectorField* f, MeshEntity* e): ElementOf(f,e) { diff --git a/apf/apfVectorElement.h b/apf/apfVectorElement.h index de4c3c447..bc503de44 100644 --- a/apf/apfVectorElement.h +++ b/apf/apfVectorElement.h @@ -18,6 +18,7 @@ class VectorField; class VectorElement : public ElementOf { public: + VectorElement(); VectorElement(VectorField* f, MeshEntity* e); VectorElement(VectorField* f, VectorElement* p); virtual ~VectorElement() {} 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..61534ea15 --- /dev/null +++ b/ma/adapt_algorithms.tex @@ -0,0 +1,128 @@ +\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} +\usepackage{caption} +\captionsetup[figure]{justification=centering, font=small} + +\title{Mesh Adapt Pseudo Code} +\author{} +\date{} + +\begin{document} + +\maketitle + +\section{Algorithms} + + +\begin{algorithm} +\caption{Adapt Top Level} +\begin{algorithmic}[1] + \State $numIterations$ = number of splits necessary to reduce the longest edge to the desired edge length + \For{\textbf{each} $numIterations$} + \State collapse all edges below desired length if possible + \State split all edges above desired edge length + \State snap all new vertices to model boundary + \EndFor + + \State for each shape below minimum quality, attempt operations to collapse shape or improve quality + +\end{algorithmic} +\end{algorithm} + + +\begin{algorithm} +\caption{Coarsen} +\begin{algorithmic}[1] + \For {every edge} + \State if the edge length is less than the desired edge length, then add both vertices to a list + \EndFor + + \While {there are unchecked vertices in list} + \State flag vertex as checked + \State attempt to collapse adjacent edges with current vertex removed, starting with the shortest edge until one succeeds + \State if an edge collapse succeeds, then remove flag from adjacent vertices so they are processed again + \EndWhile + +\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{fig:tetTypes} +\end{figure} + + +\begin{algorithm} +\caption{Snapping} +\begin{algorithmic}[1] + \While {there are any vertices are not on model boundary} + \State reposition the vertex to the model boundary, if no adjacent element is inverted then we are done + \State otherwise use inverted adjacent elements to calculate the first problem plane (FPP) + \State now we will attempt different operations to move the vertex towards the FPP until one succeeds starting with operations that are more likely to succeed + \State attempt collapsing edges shared by inverted elements that will move the vertex to the FPP + \State attempt to collapse edges on the FPP so FPP disappears + \State attempt to collapse edges shared by inverted elements to simplify the FPP + \State evaluate collapsing any edge that will move vertex towards the FPP and select the collapse that results in the highest quality elements + \State we will determine which of the problem entities in Figure 1 is preventing movement to FPP, using the method described section 6.4.2.2 of Li's Thesis + \If {\textbf{one large dihedral angle}} + \State attempt swapping the longest edge + \State attempt split-collapse on the longest edge + \EndIf + \If {\textbf{two large dihedral angles}} + \State attempt swapping either problem edge + \State attempt split-collapse on either problem edge + \State attempt double-split-collapse on both problem edges + \ElsIf {\textbf{three large dihedral angles}} + \State attempt swapping any edge in problem face + \State attempt swapping the problem face + \State attempt split-collapsing towards the FPP + \EndIf + \EndWhile + +\end{algorithmic} +\end{algorithm} + + +\begin{algorithm} +\caption{Fix Element Shapes} +\begin{algorithmic}[1] +\While {the number of low quality entities decreases} + \State we will attempt different operations and check if they increase quality + \State attempt to collapse any short edges in entity + \State we will determine which of the low quality entities in Figure 1 matches the current entity, using the method described section 6.4.2.2 of Li's Thesis + \If {\textbf{one large dihedral angle}} + \State attempt region-collapse if the tetrahedron is on the surface + \State attempt edge-swap on the longest edge in the low quality triangle + \State attempt split-collapse on the longest edge in the low quality triangle + \State attempt edge-collapse on any edge in the low quality triangle besides the longest edge + \State reposition vertices to improve quality + \EndIf + \If {\textbf{two large dihedral angles}} + \State attempt region-collapse if the tetrahedron is on the surface + \State attempt edge-swap on $M^{1}_{j_0}$ or $M^{1}_{j_1}$ + \State attempt split-collapse on $M^{1}_{j_0}$ or $M^{1}_{j_1}$ + \State attempt double-split-collapse on $M^{1}_{j_0}$ and $M^{1}_{j_1}$ + \State reposition vertices to improve quality + \ElsIf {\textbf{three large dihedral angles}} + \State attempt edge-collapse on $M^{1}_{j_0}$, $M^{1}_{j_1}$, or $M^{1}_{j_2}$ with $M^{0}_{i_3}$ deleted + \State attempt edge-swap on any edge bounding triangle made by $M^{0}_{i_0}$, $M^{0}_{i_1}$, and $M^{0}_{i_2}$ + \State attempt face-swap on the triangle made by $M^{0}_{i_0}$, $M^{0}_{i_1}$, $M^{0}_{i_2}$ + \State attempt 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 + \State attempt 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 + \State reposition vertex $M^{0}_{i_3}$ to improve quality + \EndIf +\EndWhile + +\end{algorithmic} +\end{algorithm} + +\end{document} 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..c5ce4ca00 100644 --- a/ma/ma.cc +++ b/ma/ma.cc @@ -16,33 +16,88 @@ #include "maBalance.h" #include "maLayer.h" #include "maDBG.h" +#include "maFixShape.h" #include +#include +#include +#include "apfGeometry.h" namespace ma { +void printHistogramData(std::string name, std::vector input, double min, double max, Mesh* m) +{ + const int nbins = 10; + int count[nbins] = {0}; + const double bin_size = (max-min)/(nbins*1.0); + double inputMax = 0; + 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]; + int bin = (int)std::floor((input[i] - min)/bin_size); + if (bin >= nbins) bin = nbins - 1; + if (bin < 0) bin = 0; + count[bin] += 1; + } + + inputMin = m->getPCU()->Min(inputMin); + inputMax = m->getPCU()->Max(inputMax); + for (int i = 0; i < nbins; ++i) count[i] = m->getPCU()->Add(count[i]); + + if (m->getPCU()->Self()) return; + printf("%s Min: %f, Max: %f\n", name.c_str(), inputMin, inputMax); + for (int i = 0; i < nbins; ++i) printf("%d\n", count[i]); +} + +void printHistogramStats(Adapt* a) +{ + std::vector lengths; + std::vector qualities; + ma::stats(a->mesh, a->input->sizeField, lengths, qualities, true); + printHistogramData("\nQualities:", qualities, 0, 1, a->mesh); + printHistogramData("\nLengths:", lengths, 0, MAXLENGTH+1, a->mesh); +} + void adapt(Input* in) { + apf::writeVtkFiles("mesh_tets_start", in->mesh, 3); double t0 = pcu::Time(); print(in->mesh->getPCU(), "version 2.0 !"); validateInput(in); Adapt* a = new Adapt(in); preBalance(a); + CALLGRIND_START_INSTRUMENTATION; + CALLGRIND_STOP_INSTRUMENTATION; for (int i = 0; i < in->maximumIterations; ++i) { print(a->mesh->getPCU(), "iteration %d", i); - coarsen(a); + coarsenMultiple(a); coarsenLayer(a); midBalance(a); refine(a); snap(a); } allowSplitCollapseOutsideLayer(a); - fixElementShapes(a); + // fixElementShapes(a); + fixElementShapesNew(a); cleanupLayer(a); tetrahedronize(a); printQuality(a); postBalance(a); Mesh* m = a->mesh; + + double t1 = pcu::Time(); + print(m->getPCU(), "mesh adapted in %f seconds", t1-t0); + printHistogramStats(a); + m->verify(); + ma_dbg::useFieldInfo(a, [m] { + apf::writeVtkFiles("mesh_tets_end", m, 3); + }); + delete a; // cleanup input object and associated sizefield and solutiontransfer objects if (in->ownsSizeField) @@ -50,8 +105,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); } diff --git a/ma/maAdapt.cc b/ma/maAdapt.cc index bff23cf90..adcbc407b 100644 --- a/ma/maAdapt.cc +++ b/ma/maAdapt.cc @@ -29,7 +29,7 @@ Adapt::Adapt(Input* in) input = in; mesh = in->mesh; setupFlags(this); - setupQualityCache(this); + setupCache(this); deleteCallback = 0; buildCallback = 0; sizeField = in->sizeField; @@ -52,7 +52,7 @@ Adapt::Adapt(Input* in) Adapt::~Adapt() { clearFlags(this); - clearQualityCache(this); + clearCache(this); delete refine; delete shape; } @@ -147,25 +147,49 @@ 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); +} + +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..a1a8d77d3 100644 --- a/ma/maAdapt.h +++ b/ma/maAdapt.h @@ -50,6 +50,7 @@ class Adapt Mesh* mesh; Tag* flagsTag; Tag* qualityCache; // to avoid repeated quality computations + Tag* sizeCache; // to avoid repeated size computations DeleteCallback* deleteCallback; apf::BuildCallback* buildCallback; SizeField* sizeField; @@ -76,8 +77,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..e7b99f3a7 100644 --- a/ma/maCoarsen.cc +++ b/ma/maCoarsen.cc @@ -16,6 +16,15 @@ This file contains two coarsening alogrithms. 2. coarsen: The first will create one independent set and then collapse all of the edges in that independent set. */ +/* +This file contains two coarsening alogrithms. + 1. coarsenMultiple: will repeatedly create independent sets and collapse the mesh + until there are no more short edges left. This code is currently only used for testing + since more work needs to be done to achieve better performance and have the code work in + parallel, at which point the other coarsen algorithm can be deleted. + 2. coarsen: The first will create one independent set and then collapse all of the edges + in that independent set. +*/ #include "maCoarsen.h" #include "maAdapt.h" #include "maCollapse.h" @@ -27,60 +36,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 { @@ -319,107 +279,15 @@ bool coarsen(Adapt* a) return true; } -/* - To be used after collapsing a vertex to flag adjacent vertices as NEED_NOT_COLLAPSE. This will create a - set of collapses where no two adjacent are vertices collapsed. Will also clear adjacent vertices for - collapse in next independent set since they might succeed after a collapse. -*/ -void flagIndependentSet(Adapt* a, apf::Up& adjacent, size_t& checked) -{ - for (int adj=0; adj < adjacent.n; adj++) { - Entity* vertices[2]; - a->mesh->getDownward(adjacent.e[adj],0, vertices); - for (int v = 0; v < 2; v++) { - setFlag(a, vertices[v], NEED_NOT_COLLAPSE); - if (getFlag(a, vertices[v], CHECKED)){ - clearFlag(a, vertices[v], CHECKED); //needs to be checked again in next independent set - checked--; - } - } - } -} - -struct EdgeLength -{ - Entity* edge; - double length; - bool operator<(const EdgeLength& other) const { - return length < other.length; - } -}; - -/* - Given an iterator pointing to a vertex we will collapse the shortest adjacent edge and try the next - shorted until one succeeds and then it will expand independent set. In Li's thesis it only attempts - to collapse the shortest edge, but this gave us better results. -*/ -bool collapseShortest(Adapt* a, Collapse& collapse, std::list& shortEdgeVerts, std::list::iterator& itr, size_t& checked, apf::Up& adjacent, Tag* lengthTag) -{ - Entity* vertex = *itr; - std::vector sorted; - for (int i=0; i < adjacent.n; i++) { - double length = getLength(a, lengthTag, adjacent.e[i]); - EdgeLength measured{adjacent.e[i], length}; - if (measured.length > MINLENGTH) continue; - sorted.push_back(measured); - } - if (sorted.size() == 0) { //performance optimization, will rarely result in a missed edge - itr = shortEdgeVerts.erase(itr); - return false; - } - std::sort(sorted.begin(), sorted.end()); - for (size_t i=0; i < sorted.size(); i++) { - Entity* keepVertex = getEdgeVertOppositeVert(a->mesh, sorted[i].edge, vertex); - if (!tryCollapseEdge(a, sorted[i].edge, keepVertex, collapse, adjacent)) continue; - flagIndependentSet(a, adjacent, checked); - itr = shortEdgeVerts.erase(itr); - collapse.destroyOldElements(); - return true; - } - setFlag(a, vertex, CHECKED); - checked++; - return false; -} - -/* - Iterates through shortEdgeVerts until it finds a vertex that is adjacent to an - independent set. We want our collapses to touch the independent set in order to - reduce adjacent collapses, since collapsing adjacent vertices will result in a - lower quality mesh. -*/ -bool getAdjIndependentSet(Adapt* a, std::list& shortEdgeVerts, std::list::iterator& itr, bool& independentSetStarted, apf::Up& adjacent) -{ - size_t numItr=0; - do { - numItr++; - if (itr == shortEdgeVerts.end()) itr = shortEdgeVerts.begin(); - Entity* vertex = *itr; - if (getFlag(a, vertex, CHECKED)) {itr++; continue;} //Already tried to collapse - a->mesh->getUp(vertex, adjacent); - if (!independentSetStarted) return true; - if (getFlag(a, vertex, NEED_NOT_COLLAPSE)) {itr++; continue;} //Too close to last collapse - for (int i=0; i < adjacent.n; i++) - { - Entity* opposite = getEdgeVertOppositeVert(a->mesh, adjacent.e[i], vertex); - if (getFlag(a, opposite, NEED_NOT_COLLAPSE)) return true; //Touching independent set - } - itr++; - } while (numItr < shortEdgeVerts.size()); - for (auto v : shortEdgeVerts) clearFlag(a, v, NEED_NOT_COLLAPSE); - independentSetStarted = false; - return false; -} - -//returns a list of vertices that bound a short edge and flags them -std::list getShortEdgeVerts(Adapt* a, Tag* lengthTag) +//returns a list of vertices that bound a short edge +std::list getShortEdgeVerts(Adapt* a) { std::list shortEdgeVerts; Iterator* it = a->mesh->begin(1); Entity* edge; while ((edge = a->mesh->iterate(it))) { - double length = getLength(a, lengthTag, edge); - if (length > MINLENGTH) continue; + if (getAndCacheSize(a, edge) > MINLENGTH) continue; Entity* vertices[2]; a->mesh->getDownward(edge,0,vertices); for (int i = 0; i < 2; i++) { @@ -433,37 +301,122 @@ std::list getShortEdgeVerts(Adapt* a, Tag* lengthTag) 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) +class CollapseAll : public apf::CavityOp, public DeleteCallback { - if (!a->input->shouldCoarsen) + public: + Adapt* a; + Entity* vertex; + Collapse collapse; + std::list shortEdgeVerts; + + int success=0; + + CollapseAll(Adapt* adapt) : apf::CavityOp(adapt->mesh,true), DeleteCallback(adapt) + { + a = adapt; + collapse.Init(a); + shortEdgeVerts = getShortEdgeVerts(a); + } + + ~CollapseAll() + { + ma::clearFlagFromDimension(a, CHECKED, 0); + } + + int getTargetDimension() { return 0; } + void call(Entity* e) { movedByDeletion = true; } + 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(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; + } + }; + + /* + 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(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; - double t0 = pcu::Time(); - Tag* lengthTag = a->mesh->createDoubleTag("edge_length", 1); - std::list shortEdgeVerts = getShortEdgeVerts(a, lengthTag); + } - Collapse collapse; - collapse.Init(a); - int success = 0; - size_t checked = 0; - bool independentSetStarted = false; - std::list::iterator itr = shortEdgeVerts.begin(); - while (checked < shortEdgeVerts.size()) + 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; + } + + void apply() { apf::Up adjacent; - if (!getAdjIndependentSet(a, shortEdgeVerts, itr, independentSetStarted, adjacent)) continue; - if (collapseShortest(a, collapse, shortEdgeVerts, itr, checked, adjacent, lengthTag)) { - independentSetStarted=true; - success++; - } + a->mesh->getUp(vertex, adjacent); + setFlag(a, vertex, CHECKED); + if (collapseShortest(vertex, adjacent)) success++; } - ma::clearFlagFromDimension(a, NEED_NOT_COLLAPSE | CHECKED, 0); - a->mesh->destroyTag(lengthTag); +}; + +bool coarsenMultiple(Adapt* a) +{ + if (!a->input->shouldCoarsen) return false; + double t0 = pcu::Time(); + CollapseAll collapseAll(a); + int totalSuccess = 0; + int success = 0; + 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/maCollapse.cc b/ma/maCollapse.cc index af13ba324..4b0363018 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,65 @@ 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 +96,32 @@ 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); + adapt->mesh->setPoint(vertToCollapse, 0, target); + bool worse = false; + for (Entity* e : elementsToKeep) + if (adapt->shape->getQuality(e) < qualityToBeat) {worse = true; break;} + adapt->mesh->setPoint(vertToCollapse, 0, prev); + return worse; +} + bool Collapse::tryThisDirectionNoCancel(double qualityToBeat) { PCU_ALWAYS_ASSERT( ! adapt->mesh->isShared(vertToCollapse)); @@ -158,9 +246,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 +451,52 @@ void Collapse::computeElementSets() APF_ITERATE(Upward,adjacent,it) if ( ! elementsToCollapse.count(*it)) elementsToKeep.insert(*it); - PCU_ALWAYS_ASSERT(elementsToKeep.size()); + // PCU_ALWAYS_ASSERT(elementsToKeep.size()); TODO: change computeElementSets to bool function that fails when this is false } -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; + Entity* faceToReplace; + + 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 +509,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..67b842f61 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; + 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/maFaceSwap.cc b/ma/maFaceSwap.cc new file mode 100644 index 000000000..daa326527 --- /dev/null +++ b/ma/maFaceSwap.cc @@ -0,0 +1,236 @@ +#include "maAdapt.h" +#include "maSnapper.h" +#include "apfGeometry.h" +#include "maShapeHandler.h" +#include "maDBG.h" + +//TODO: add comments + +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; //TODO: TEST Two2Two CASE + 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; + } + } + +} \ No newline at end of file diff --git a/ma/maFaceSwap.h b/ma/maFaceSwap.h new file mode 100644 index 000000000..50b41f002 --- /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 \ No newline at end of file diff --git a/ma/maFixShape.cc b/ma/maFixShape.cc new file mode 100644 index 000000000..9db1ee436 --- /dev/null +++ b/ma/maFixShape.cc @@ -0,0 +1,440 @@ + +/****************************************************************************** + + 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 + +namespace ma { + +int markBadQualityNew(Adapt* a) +{ + Mesh* m = a->mesh; + Iterator* it; + Entity* e; + it = m->begin(m->getDimension()); + int total = 0; + while ((e = m->iterate(it))) { + 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 total; +} + +FixShape::FixShape(Adapt* adapt) : splitCollapse(adapt), doubleSplitCollapse(adapt), faceSplitCollapse(adapt), split(adapt), reposition(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]; + a->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]; + a->mesh->getDownward(edge, 0, verts); + for (int v=0; v<2; v++) { + apf::Up adjEdges; + a->mesh->getUp(verts[v], adjEdges); + for (int e=0; emesh, adjEdges.e[e], verts[v]); + if (!mesh->isOwned(keep)) continue; // TODO: improvement to requestLocality should remove use for this function + if (mesh->isShared(keep)) continue; // TODO: improvement to requestLocality should remove use for this function + 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]; + a->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]; + a->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 < .09; +} + +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]; + a->mesh->getDownward(worstTriangle, 1, edges); + Entity* longestEdge = getLongestEdge(edges); + Entity* oppositeVert = getTriVertOppositeEdge(a->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]; + a->mesh->getDownward(tet, 0, verts); + Entity* vert = verts[0]; + Entity* face = getTetFaceOppositeVert(a->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]; + a->mesh->getDownward(tet, 2, faces); + for (int i=0; i<4; i++) { + if (isLowInHigh(a->mesh, faces[i], problemEnts[0])) { + Entity* v0 = getTriVertOppositeEdge(a->mesh, faces[i], problemEnts[0]); + if (splitCollapse.run(problemEnts[0], v0)) {numEdgeSplitCollapse++; return true;} + } + else { + Entity* v0 = getTriVertOppositeEdge(a->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]; + a->mesh->getDownward(tet, 1, tetEdges); + for (int i=0; i<6; i++) + if (!isLowInHigh(a->mesh, problemEnts[0], tetEdges[i])) { + Entity* verts[2]; + a->mesh->getDownward(tetEdges[i], 0, verts); + Entity* keep = isLowInHigh(a->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]; + a->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(a->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 a->mesh->getPCU()->Add(val); +} +void FixShape::printNumOperations() +{ + print(a->mesh->getPCU(), "shape operations: \n collapses %17d\n edge swaps %16d\n reposition %16d\n split reposition %9d\n double split collapse %d\n " + "edge split collapses %5d\n face split collapses %3d\n region collapses %7d\n face swaps %12d\n ", + collect(numCollapse), collect(numEdgeSwap), collect(numReposition), collect(numSplitReposition), collect(numDoubleSplitCollapse), + collect(numEdgeSplitCollapse), collect(numFaceSplitCollapse), collect(numRegionCollapse), collect(numFaceSwap)); +} +void FixShape::printBadTypes() +{ + print(a->mesh->getPCU(), "bad shape types: \n oneShortEdge \t%d\n twoShortEdges \t%d\n threeShortEdges \t%d\n " + "moreShortEdges \t%d\n oneLargeAngle \t%d\n twoLargeAngle \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]; + a->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(a->mesh->getPCU(), "Worst is short\n"); + else if (isOneLargeAngle(problemTet, worstTri)) print(a->mesh->getPCU(), "Worst is one large angle\n"); + else if (isTwoLargeAngles(problemTet, problemEnts)) print(a->mesh->getPCU(), "Worst is two large angles\n"); + else print(a->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 = a->mesh->begin(3); + Entity* tet; + while ((tet = a->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", a->mesh, 3); + apf::writeVtkFiles("mesh_faces", a->mesh, 2); + apf::writeVtkFiles("mesh_edges", a->mesh, 1); + }); +} + +void FixShape::printNumTypes() +{ + resetCounters(); + Entity* tet; + Iterator* it = a->mesh->begin(3); + double worstQual = 1; + Entity* worstShape; + while ((tet = a->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 fixElementShapesNew(Adapt* a) +{ + if ( ! a->input->shouldFixShape) + return; + double t0 = pcu::Time(); + bool oldForce = a->input->shouldForceAdaptation; + a->input->shouldForceAdaptation = false; + int count = markBadQualityNew(a); + print(a->mesh->getPCU(), "loop %d: of shape correction loop: #bad elements %d", 0, count); + FixShape fixShape(a); + // fixShape.printNumTypes(); + int originalCount = count; + int prev_count; + int iter = 0; + do { + if (!count) break; + prev_count = count; + applyOperator(a,&fixShape); + if (a->mesh->getDimension() == 3) + snap(a); + count = markBadQualityNew(a); + midBalance(a); // balance the mesh to avoid empty parts + iter++; + print(a->mesh->getPCU(), "loop %d: bad shapes went from %d to %d", iter, prev_count, count); + } while(count < prev_count); + a->input->shouldForceAdaptation = oldForce; + double t1 = pcu::Time(); + print(a->mesh->getPCU(), "bad shapes down from %d to %d in %f seconds", + originalCount,count,t1-t0); + fixShape.printNumOperations(); + fixShape.printNumTypes(); + 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..1fbb8c1e1 --- /dev/null +++ b/ma/maFixShape.h @@ -0,0 +1,85 @@ +#ifndef MA_SHAPE_NEW +#define MA_SHAPE_NEW + +#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; + Collapse collapse; + SingleSplitCollapse splitCollapse; + DoubleSplitCollapse doubleSplitCollapse; + FaceSplitCollapse faceSplitCollapse; + RepositionVertex reposition; + EdgeSwap* edgeSwap; + Splits split; + Entity* badTet; + + 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 fixElementShapesNew(Adapt* a); +} + +#endif \ No newline at end of file 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..9ca8a8c87 100644 --- a/ma/maQuality.cc +++ b/ma/maQuality.cc @@ -20,6 +20,15 @@ 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]; @@ -156,8 +165,8 @@ double measureTetQuality(Mesh* m, SizeField* f, Entity* tet, bool useMax) m->getDownward(tet,1,e); double l[6]; for (int i=0; i < 6; ++i) - l[i] = qMeasure(m, e[i], Q); - double V = qMeasure(m, tet, Q); + l[i] = f->measure(e[i], Q); + double V = f->measure(tet, Q); double s=0; for (int i=0; i < 6; ++i) s += l[i]*l[i]; @@ -181,6 +190,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/maReposition.cc b/ma/maReposition.cc index 42be04314..db8f7120b 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,179 @@ #include #include +#include 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 (size_t 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; + mesh->setPoint(vertex, 0, target); + findInvalid(); + if (invalid.n == 0) 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; + + 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; + + } while ((left-right).getLength() > tol); + return leftValue; +} + +bool RepositionVertex::moveToImproveQuality(Entity* vertex) +{ + if (mesh->getModelType(mesh->toModel(vertex)) < 2) return false; //TODO: remove limitation + 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); }; + worstQuality = goldenSearch(getQuality, prevPosition, target, 0.000001); + 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; +} + +double RepositionVertex::findShortestEdge(Vector position) +{ + double quality = findWorstShape(position); + if (quality < startingQuality && quality < adapt->input->goodQuality) return -1; + + double shortest = 1000; + for (int i=0; isizeField->measure(adjEdges.e[i]); + if (length < shortest) shortest = length; + } + return shortest; +} + +void RepositionVertex::moveToImproveShortEdges(Entity* vertex) +{ + int dim = mesh->getModelType(mesh->toModel(vertex)); + if (dim == 0) return; + if (!init(vertex)) return; + + this->startingQuality = findWorstShape(prevPosition); + Vector center = modelCenter(); + const auto getLength = [this](const Vector& pos) { return this->findShortestEdge(pos); }; + goldenSearch(getLength, prevPosition, center, 0.000001); +} + +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..aaf0165e1 100644 --- a/ma/maReposition.h +++ b/ma/maReposition.h @@ -2,9 +2,40 @@ #define MA_REPOSITION_H #include "maMesh.h" +#include "ma.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; + + 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..1a75c248a 100644 --- a/ma/maShape.cc +++ b/ma/maShape.cc @@ -877,7 +877,7 @@ 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..db3df76af 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 */ diff --git a/ma/maSize.cc b/ma/maSize.cc index e596d0788..106fa858e 100644 --- a/ma/maSize.cc +++ b/ma/maSize.cc @@ -13,6 +13,8 @@ #include #include #include +#include "apfVectorElement.h" +#include "maAdapt.h" namespace ma { @@ -59,6 +61,14 @@ double IdentitySizeField::measure(Entity* e) return x; } +double IdentitySizeField::measure(Entity* e, Matrix const& Q) +{ + apf::MeshElement* me = apf::createMeshElement(mesh,e); + double x = apf::measure(me); + apf::destroyMeshElement(me); + return x; +} + bool IdentitySizeField::shouldSplit(Entity*) { return false; @@ -184,6 +194,7 @@ class SizeFieldIntegrator : public apf::Integrator { meshElement = me; dimension = apf::getDimension(me); + measurement = 0; } void atPoint(Vector const& p , double w, double ) { @@ -205,14 +216,48 @@ 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 = numIntegrationPoints.try_emplace(dim, countIntPoints(&me,integrationOrder)).first->second; + double w = integrationWeight.try_emplace(dim, getIntWeight(&me,integrationOrder,0)).first->second; + for (int p=0; p < np; ++p) { + Vector point = integrationPoint.try_emplace({dim, p}, getIntPoint(&me,integrationOrder,p)).first->second; + 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; + } + double measure(Entity* e, Matrix const& Q) + { + 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 = numIntegrationPoints.try_emplace(dim, countIntPoints(&me,integrationOrder)).first->second; + double w = integrationWeight.try_emplace(dim, getIntWeight(&me,integrationOrder,0)).first->second; + for (int p=0; p < np; ++p) { + Vector point = integrationPoint.try_emplace({dim, p}, getIntPoint(&me,integrationOrder,p)).first->second; + 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 +274,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 +412,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,6 +442,7 @@ struct AnisoSizeField : public MetricSizeField Vector const& xi, Matrix& Q) { + //TODO: cache these elements apf::Element* hElement = apf::createElement(hField,me); apf::Element* rElement = apf::createElement(rField,me); Vector h; @@ -454,12 +500,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 +529,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 +552,10 @@ 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); @@ -525,11 +569,11 @@ struct LogAnisoSizeField : public MetricSizeField 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 +618,7 @@ struct LogAnisoSizeField : public MetricSizeField return apf::getShape(logMField)->hasNodesIn(dimension); } apf::NewArray fieldVal; + apf::Element logMElement; apf::Field* logMField; LogMEval logMEval; }; @@ -615,10 +660,34 @@ struct IsoUserField : public IsoSizeField FieldReader reader; }; +static void clampSizeField(Mesh*m, apf::Field* sizes) +{ + Vector lower; + Vector upper; + getBoundingBox(m, lower, upper); + 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] > upper[0]) h[0] = upper[0]; + if (h[1] > upper[1]) h[1] = upper[1]; + if (h[2] > upper[2]) h[2] = upper[2]; + apf::setVector(sizes, ent, i, h); + } + } + } +} + 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 +720,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..f3e4ebfc4 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, @@ -32,6 +33,7 @@ class SizeField public: virtual ~SizeField(); virtual double measure(Entity* e) = 0; + virtual double measure(Entity* e, Matrix const& Q) = 0; virtual bool shouldSplit(Entity* edge) = 0; virtual bool shouldCollapse(Entity* edge) = 0; virtual void interpolate( @@ -57,6 +59,7 @@ struct IdentitySizeField : public SizeField { IdentitySizeField(Mesh* m); double measure(Entity* e); + double measure(Entity* e, Matrix const& Q); bool shouldSplit(Entity*); bool shouldCollapse(Entity*); void interpolate( @@ -112,6 +115,7 @@ 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); diff --git a/ma/maSnap.cc b/ma/maSnap.cc index 3c4d92b76..5aef830d2 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 @@ -39,6 +46,8 @@ static bool isCapstone(apf::Mesh* m) { #endif } +std::map> boundingBoxes; + static size_t isSurfUnderlyingFaceDegenerate( apf::Mesh* m, Model* g, // this the model entity in question @@ -53,11 +62,20 @@ static size_t isSurfUnderlyingFaceDegenerate( Vector bmin; Vector bmax; - m->boundingBox(g, bmin, bmax); + auto it = boundingBoxes.find(g); + if (it != boundingBoxes.end()) { + bmin = it->second.first; + bmax = it->second.second; + } + else{ + m->boundingBox(g, bmin, bmax); + boundingBoxes.insert(it, {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]); diff --git a/ma/maSnapper.cc b/ma/maSnapper.cc index 2f3a426d2..5f032d4cf 100644 --- a/ma/maSnapper.cc +++ b/ma/maSnapper.cc @@ -15,24 +15,37 @@ * 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); edgeSwap = makeEdgeSwap(a); + edgeSwap = makeEdgeSwap(a); vert = 0; } @@ -83,23 +96,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 @@ -290,6 +305,8 @@ int getTetStats(Adapt* a, Entity* vert, Entity* face, Entity* region, Entity* en return bit; } +static int debugprint = 0; + /* We perform this last to make sure that we have a simple region where we can determine the best operation to perform and because we want to avoid creating more vertices to snap since @@ -300,6 +317,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 +331,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 +340,20 @@ 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 (runFaceSwap(adapt, ents[0], false)) { numSwap++; return true; } + if (splitCollapse.run(ents[1], FPP->vert, qual)) { numSplitCollapse++; return true; } } return false; } @@ -560,43 +554,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 +573,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); @@ -624,7 +591,7 @@ bool Snapper::run() if (!success) FPP = getFPP(adapt, vert, snapTag, invalid); if (!success) success = tryCollapseToVertex(FPP); if (!success) success = tryReduceCommonEdges(FPP); - if (!success) success = tryCollapseTetEdges(FPP); + if (!success) success = tryCollapseTetEdges(FPP); //TODO causing holes in mesh if (!success) success = trySwapOrSplit(FPP); } @@ -634,7 +601,7 @@ bool Snapper::run() clearFlag(adapt, vert, SNAP); } #if defined(DEBUG_FPP) - if (!success && ++DEBUGFAILED == 1) printFPP(adapt, FPP); + // if (!success && ++DEBUGFAILED == 2) printFPP(adapt, FPP); #endif if (FPP) delete FPP; return success; @@ -951,4 +918,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/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/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_ma_test.cc b/test/aniso_ma_test.cc index e115d1c58..dcc4723a2 100644 --- a/test/aniso_ma_test.cc +++ b/test/aniso_ma_test.cc @@ -58,7 +58,6 @@ int main(int argc, char** argv) in->shouldRunMidParma = true; in->shouldRunPostParma = true; in->shouldRefineLayer = true; - in->goodQuality = 0.2; ma::adapt(in); m->verify(); if (logInterpolation) 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