Skip to content

Commit

Permalink
Improved MeshBoundaryFactory
Browse files Browse the repository at this point in the history
  • Loading branch information
@regislebrun
regislebrun committed Sep 25, 2023
1 parent 0c3b2ab commit 26846d6
Showing 1 changed file with 130 additions and 86 deletions.
216 changes: 130 additions & 86 deletions lib/src/Base/Geom/MeshBoundaryFactory.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,11 @@
* along with this library. If not, see <http://www.gnu.org/licenses/>.
*
*/
#include <unordered_map>

#include "openturns/PersistentObjectFactory.hxx"
#include "openturns/MeshBoundaryFactory.hxx"
#include <unordered_map>
#include "openturns/TBBImplementation.hxx"

BEGIN_NAMESPACE_OPENTURNS

Expand Down Expand Up @@ -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<UnsignedInteger> & 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

/*
Expand Down Expand Up @@ -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<UnsignedInteger, UnsignedInteger> oldToNewIndices;
// Preallocate enough space
oldToNewIndices.reserve(vertices.getSize());
std::unordered_map<UnsignedInteger, UnsignedInteger>::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
Expand All @@ -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
Expand All @@ -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

0 comments on commit 26846d6

Please sign in to comment.