From 26846d6991f1d53179f16a98c92583aac4f2b985 Mon Sep 17 00:00:00 2001 From: regislebrun Date: Mon, 25 Sep 2023 11:23:27 +0200 Subject: [PATCH] Improved MeshBoundaryFactory --- lib/src/Base/Geom/MeshBoundaryFactory.cxx | 216 +++++++++++++--------- 1 file changed, 130 insertions(+), 86 deletions(-) diff --git a/lib/src/Base/Geom/MeshBoundaryFactory.cxx b/lib/src/Base/Geom/MeshBoundaryFactory.cxx index e063e87ec0f..4687b3da2ca 100644 --- a/lib/src/Base/Geom/MeshBoundaryFactory.cxx +++ b/lib/src/Base/Geom/MeshBoundaryFactory.cxx @@ -18,9 +18,11 @@ * along with this library. If not, see . * */ +#include + #include "openturns/PersistentObjectFactory.hxx" #include "openturns/MeshBoundaryFactory.hxx" -#include +#include "openturns/TBBImplementation.hxx" BEGIN_NAMESPACE_OPENTURNS @@ -87,6 +89,107 @@ struct FaceCompare } }; // FaceCompare +struct ComputeNormalPolicy +{ + const Sample & vertices_; + const Collection< Indices > & boundaryFaces_; + const UnsignedInteger dimension_; + const UnsignedInteger baseVertexIndex_; + const Scalar offset_; + Sample & boundaryVertices_; + IndicesCollection & boundarySimplices_; + + ComputeNormalPolicy(const Sample & vertices, + const Collection< Indices > & boundaryFaces, + const UnsignedInteger baseVertexIndex, + const Scalar offset, + Sample & boundaryVertices, + IndicesCollection & boundarySimplices) + : vertices_(vertices) + , boundaryFaces_(boundaryFaces) + , dimension_(boundaryVertices.getDimension()) + , baseVertexIndex_(baseVertexIndex) + , offset_(offset) + , boundaryVertices_(boundaryVertices) + , boundarySimplices_(boundarySimplices) + { + // Nothing to do + } + + inline void operator()(const TBBImplementation::BlockedRange & r) const + { + Indices face; + Matrix U; + Matrix VT; + SquareMatrix M(dimension_); + SquareMatrix simplexMatrix(dimension_ + 1); + + for (UnsignedInteger i = r.begin(); i != r.end(); ++ i) + { + UnsignedInteger newVertexIndex = baseVertexIndex_ + i; + face = boundaryFaces_[i]; + const UnsignedInteger remainingVertexIndex = face[dimension_]; + // Build the face matrix & compute the boundary vertex to thicken the face + // In order to avoid roundoff we first compute the center of the face, + // then we compute the hyperplane of the face translated to the center + Point center(dimension_); + for (UnsignedInteger j = 0; j < dimension_; ++j) + { + const UnsignedInteger vertexIndex = face[j]; + for (UnsignedInteger k = 0; k < dimension_; ++k) + center[k] += boundaryVertices_(vertexIndex, k); + } // j + center /= dimension_; + // Now the centered matrix + for (UnsignedInteger j = 0; j < dimension_; ++j) + { + const UnsignedInteger vertexIndex = face[j]; + for (UnsignedInteger k = 0; k < dimension_; ++k) + M(j, k) = boundaryVertices_(vertexIndex, k) - center[k]; + } // j + // Add a new vertex + try + { + // Compute the face normal + (void) M.computeSVD(U, VT, true, false); + Point normal(dimension_); + // The normal is the last row of VT + for (UnsignedInteger j = 0; j < dimension_; ++j) + normal[j] = VT(dimension_ - 1, j); + // Check if the vertex removed from the simplex to form a face is on the same side + // of the plane as the one pointed by the normal + if (normal.dot(vertices_[remainingVertexIndex] - center) < 0.0) + center += normal * offset_; + else + center -= normal * offset_; + } + catch (const Exception &) + { + // The face was in fact degenerated, not a (dimension - 1) entity + // Nothing to do + } + boundaryVertices_[newVertexIndex] = center; + boundarySimplices_(i, dimension_) = newVertexIndex; + // Fix the simplex orientation + for (UnsignedInteger j = 0; j <= dimension_; ++j) + { + const UnsignedInteger vertexJ = boundarySimplices_(i, j); + for (UnsignedInteger k = 0; k < dimension_; ++k) simplexMatrix(k, j) = boundaryVertices_(vertexJ, k); + simplexMatrix(dimension_, j) = 1.0; + } + Scalar sign = 0.0; + (void) simplexMatrix.computeLogAbsoluteDeterminant(sign, false); + // In odd dimension the positive orientation is for a negative determinant of + // the simplex matrix + if ((sign > 0.0) != (dimension_ % 2 == 1)) + { + IndicesCollection::iterator cit = boundarySimplices_.begin_at(i); + std::swap(*cit, *(cit + 1)); + } + } // i + } // operator () +}; // struct ComputeNormalPolicy + } // Anonymous namespace /* @@ -154,33 +257,31 @@ Mesh MeshBoundaryFactory::build(const Mesh & mesh, const UnsignedInteger boundaryVerticesReserve = (offset == 0.0 ? 0 : boundaryFaces.getSize()); Sample boundaryVertices(vertices.getSize() + boundaryVerticesReserve, dimension); UnsignedInteger newVertexIndex = 0; - // This unordered map is used to compute the new indices of the boundary vertices - std::unordered_map oldToNewIndices; - // Preallocate enough space - oldToNewIndices.reserve(vertices.getSize()); - std::unordered_map::const_iterator itIndex; - for (UnsignedInteger i = 0; i < boundaryFaces.getSize(); ++i) { - for (UnsignedInteger j = 0; j < dimension; ++j) - { - // Get the old vertex index of the jth vertex of the ith boundary face - const UnsignedInteger oldVertexIndex = boundaryFaces[i][j]; - itIndex = oldToNewIndices.find(oldVertexIndex); - // If the vertex has not been seen so far, insert it into the sample of boundary vertices - if (itIndex == oldToNewIndices.end()) - { - boundaryFaces[i][j] = newVertexIndex; - for (UnsignedInteger k = 0; k < dimension; ++k) - boundaryVertices(newVertexIndex, k) = vertices(oldVertexIndex, k); - oldToNewIndices.insert({oldVertexIndex, newVertexIndex}); - ++newVertexIndex; - } - else + const UnsignedInteger nbVertices = vertices.getSize(); + Indices oldToNewIndices(nbVertices, nbVertices); + for (UnsignedInteger i = 0; i < boundaryFaces.getSize(); ++i) { - boundaryFaces[i][j] = itIndex->second; - } // if - } //j - } // i + for (UnsignedInteger j = 0; j < dimension; ++j) + { + // Get the old vertex index of the jth vertex of the ith boundary face + const UnsignedInteger oldVertexIndex = boundaryFaces[i][j]; + // If the vertex has not been seen so far, insert it into the sample of boundary vertices + if (oldToNewIndices[oldVertexIndex] == nbVertices) + { + boundaryFaces[i][j] = newVertexIndex; + for (UnsignedInteger k = 0; k < dimension; ++k) + boundaryVertices(newVertexIndex, k) = vertices(oldVertexIndex, k); + oldToNewIndices[oldVertexIndex] = newVertexIndex; + ++newVertexIndex; + } + else + { + boundaryFaces[i][j] = oldToNewIndices[oldVertexIndex]; + } // if + } //j + } // i + } // In a dedicated scope to allow for the liberation of oldToNewIndices // Resize the boundary vertices boundaryVertices.erase(boundaryVertices.getImplementation()->begin() + newVertexIndex + boundaryVerticesReserve, boundaryVertices.getImplementation()->end()); // Now, create the face in the boundary mesh @@ -189,7 +290,6 @@ Mesh MeshBoundaryFactory::build(const Mesh & mesh, // The last index is repeated in order to indicate that it is a surface mesh boundarySimplices(i, dimension) = boundarySimplices(i, dimension - 1); - Mesh result; // III. If offset is not zero, compute the normal of each face // using a SVD decomposition of the vertices matrix. // Here we use the initial vertices to get the vertex of the simplex @@ -198,67 +298,11 @@ Mesh MeshBoundaryFactory::build(const Mesh & mesh, // domain. if (offset != 0.0) { - // Preallocate the matrix used for the SVD computation - SquareMatrix M(dimension); - newVertexIndex = boundaryVertices.getSize() - boundaryVerticesReserve; - for (UnsignedInteger i = 0; i < boundaryFaces.getSize(); ++i) - { - face = boundaryFaces[i]; - const UnsignedInteger remainingVertexIndex = face[dimension]; - // Build the face matrix & compute the boundary vertex to thicken the face - // In order to avoid roundoff we first compute the center of the face, - // then we compute the hyperplane of the face translated to the center - Point center(dimension); - for (UnsignedInteger j = 0; j < dimension; ++j) - { - const UnsignedInteger vertexIndex = face[j]; - for (UnsignedInteger k = 0; k < dimension; ++k) - center[k] += boundaryVertices(vertexIndex, k); - } // j - center /= dimension; - // Now the centered matrix - for (UnsignedInteger j = 0; j < dimension; ++j) - { - const UnsignedInteger vertexIndex = face[j]; - for (UnsignedInteger k = 0; k < dimension; ++k) - M(j, k) = boundaryVertices(vertexIndex, k) - center[k]; - } // j - // Add a new vertex - try - { - // Compute the face normal - Matrix U; - Matrix VT; - (void) M.computeSVD(U, VT, true, false); - Point normal(dimension); - // The normal is the last row of VT - for (UnsignedInteger j = 0; j < dimension; ++j) - normal[j] = VT(dimension - 1, j); - // Check if the vertex removed from the simplex to form a face is on the same side - // of the plane as the one pointed by the normal - if (normal.dot(vertices[remainingVertexIndex] - center) < 0.0) - center += normal * offset; - else - center -= normal * offset; - } - catch (const Exception &) - { - // The face was in fact degenerated, not a (dimension - 1) entity - // Nothing to do - } - boundaryVertices[newVertexIndex] = center; - boundarySimplices(i, dimension) = newVertexIndex; - ++newVertexIndex; - } // i + const ComputeNormalPolicy policy(vertices, boundaryFaces, boundaryVertices.getSize() - boundaryVerticesReserve, offset, boundaryVertices, boundarySimplices); + TBBImplementation::ParallelFor(0, boundaryFaces.getSize(), policy); } // offset != 0 - result = Mesh(boundaryVertices, boundarySimplices); - - // Fix the orientation - SquareMatrix matrix(dimension + 1); - for (UnsignedInteger i = 0; i < boundarySimplices.getSize(); ++i) - result.fixOrientation(i, matrix); // Return the boundary mesh - return result; + return Mesh(boundaryVertices, boundarySimplices, false); } END_NAMESPACE_OPENTURNS