Added MeshBoundaryFactory

This class allows to build the mesh of the boundary of a given mesh, with an optional offset.

lib/src/Base/Geom/CMakeLists.txt

@@ -11,6 +11,7 @@ ot_add_source_file (LevelSet.cxx)
 ot_add_source_file (LevelSetMesher.cxx)
 ot_add_source_file (Mesh.cxx)
 ot_add_source_file (MeshDomain.cxx)
+ot_add_source_file (MeshBoundaryFactory.cxx)
 ot_add_source_file (RegularGrid.cxx)
 ot_add_source_file (DomainComplement.cxx)
 ot_add_source_file (DomainIntersection.cxx)
@@ -27,6 +28,7 @@ ot_install_header_file (LevelSet.hxx)
 ot_install_header_file (LevelSetMesher.hxx)
 ot_install_header_file (Mesh.hxx)
 ot_install_header_file (MeshDomain.hxx)
+ot_install_header_file (MeshBoundaryFactory.hxx)
 ot_install_header_file (RegularGrid.hxx)
 ot_install_header_file (DomainComplement.hxx)
 ot_install_header_file (DomainIntersection.hxx)

lib/src/Base/Geom/MeshBoundaryFactory.cxx

@@ -0,0 +1,263 @@
+//                                               -*- C++ -*-
+/**
+ *  @brief Boundary extraction algorithm for meshes
+ *
+ *  Copyright 2005-2023 Airbus-EDF-IMACS-ONERA-Phimeca
+ *
+ *  This library is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This library is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Lesser General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public License
+ *  along with this library.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+#include "openturns/PersistentObjectFactory.hxx"
+#include "openturns/MeshBoundaryFactory.hxx"
+
+BEGIN_NAMESPACE_OPENTURNS
+
+CLASSNAMEINIT(MeshBoundaryFactory)
+static const Factory<MeshBoundaryFactory> Factory_MeshBoundaryFactory;
+
+
+/* Default constructor */
+MeshBoundaryFactory::MeshBoundaryFactory()
+  : PersistentObject()
+{
+  // Nothing to do
+}
+
+/* Virtual constructor */
+MeshBoundaryFactory * MeshBoundaryFactory::clone() const
+{
+  return new MeshBoundaryFactory(*this);
+}
+
+/* String converter */
+String MeshBoundaryFactory::__repr__() const
+{
+  OSS oss(true);
+  oss << "class=" << MeshBoundaryFactory::GetClassName();
+  return oss;
+}
+
+/* String converter */
+String MeshBoundaryFactory::__str__(const String & ) const
+{
+  return __repr__();
+}
+
+/* Here is the interface that all derived class must implement */
+
+
+namespace
+{
+struct FaceHasher
+{
+  UnsignedInteger operator()(const Indices &face) const
+  {
+    UnsignedInteger seed = face.getSize() - 1;
+    for (UnsignedInteger i = 0; i < face.getSize() - 1; ++i)
+    {
+      UnsignedInteger x = face[i];
+      x = ((x >> 16) ^ x) * 0x45d9f3b;
+      x = ((x >> 16) ^ x) * 0x45d9f3b;
+      x = (x >> 16) ^ x;
+      seed ^= x + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+    }
+    return seed;
+  }
+}; // FaceHasher
+
+struct FaceCompare
+{
+  Bool operator()(const Indices &face1, const Indices &face2) const
+  {
+    // The two faces are of size (dimension+1), and we test
+    // only the dimension first components for equality
+    return std::equal(face1.begin(), face1.end() - 1, face2.begin());
+  }
+}; // FaceCompare
+
+} // Anonymous namespace
+
+/*
+  Amélioration des performances:
+  1) Etape de préparation des faces.
+  + utiliser une IndicesCollection, pour chaque simplex trier ses numéros de noeuds
+    sans création d'Indices temporaires puis injecter ces numéros de noeuds
+ */
+Mesh MeshBoundaryFactory::build(const Mesh & mesh,
+                                const Scalar offset) const
+{
+  const UnsignedInteger dimension = mesh.getDimension();
+  // I. build the list of unique faces
+  const UnsignedInteger simplicesNumber = mesh.getSimplicesNumber();
+  IndicesCollection simplices(mesh.getSimplices());
+  // a) generate the faces given a simplex: for one simplex we get (dimension+1) faces
+  IndicesCollection simplex2Faces(dimension + 1, dimension + 1);
+  for (UnsignedInteger i = 0; i <= dimension; ++i)
+  {
+    for (UnsignedInteger j = 0; j < i; ++j)
+      simplex2Faces(i, j) = j;
+    for (UnsignedInteger j = i; j < dimension; ++j)
+      simplex2Faces(i, j) = j + 1;
+    simplex2Faces(i, dimension) = i;
+  }
+  // b) for all the faces of all the simplices, count how many time they appear
+  //    using an unordered map. Here we reserve the maximum possible size in order
+  //    to avoid costly memory reallocation and rehash
+  std::unordered_map<Indices, UnsignedInteger, FaceHasher, FaceCompare> faces;
+  faces.reserve(simplices.getSize() * (dimension + 1));
+  Indices face(dimension + 1);
+  std::unordered_map<Indices, UnsignedInteger, FaceHasher, FaceCompare>::iterator itFaces;
+  for (UnsignedInteger i = 0; i < simplices.getSize(); ++i)
+  {
+    for (UnsignedInteger j = 0; j < dimension + 1; ++j)
+    {
+      for (UnsignedInteger k = 0; k < dimension + 1; ++k)
+        face[k] = simplices(i, simplex2Faces(j, k));
+      // Sort only the indices related to the face vertices, not the last one
+      std::sort(face.begin(), face.end() - 1);
+      itFaces = faces.find(face);
+      if (itFaces == faces.end())
+        faces[face] = 1;
+      else
+        itFaces->second += 1;
+    } // j
+  } // i
+  // c) now we can detect the boundary faces: they have a count equal to 1 in the map
+  // We use a collection of indices instead of an IndicesCollection here because the size
+  // of the collection is not known in advance
+  Collection< Indices > boundaryFaces(simplicesNumber * (dimension + 1), Indices(dimension + 1));
+  UnsignedInteger boundaryFaceIndex = 0;
+  for (std::unordered_map<Indices, UnsignedInteger, FaceHasher, FaceCompare>::const_iterator it = faces.begin(); it != faces.end(); ++it)
+    if (it->second == 1)
+    {
+      boundaryFaces[boundaryFaceIndex] = it->first;
+      ++boundaryFaceIndex;
+    } // if
+  // Remove the unused space
+  boundaryFaces.erase(boundaryFaces.begin() + boundaryFaceIndex, boundaryFaces.end());
+  // II. Create the boundary simplices and vertices
+  const Sample vertices(mesh.getVertices());
+  // a) Preallocate enough space to store the boundary vertices and the offset vertex (if asked for)
+  //    Compute the renumbering of the boundary vertices on the fly
+  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
+      {
+        boundaryFaces[i][j] = itIndex->second;
+      } // if
+    } //j
+  } // i
+  // Resize the boundary vertices
+  boundaryVertices.erase(boundaryVertices.getImplementation()->begin() + newVertexIndex + boundaryVerticesReserve, boundaryVertices.getImplementation()->end());
+  // Now, create the face in the boundary mesh
+  IndicesCollection boundarySimplices(boundaryFaces);
+  for (UnsignedInteger i = 0; i < boundarySimplices.getSize(); ++i)
+    // 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
+  //      from which the face has been extracted not in the face. It allows us
+  //      to make a distinction between the interior and the exterior of the
+  //      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
+  } // 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;
+}
+
+END_NAMESPACE_OPENTURNS

lib/src/Base/Geom/openturns/MeshBoundaryFactory.hxx

@@ -0,0 +1,73 @@
+//                                               -*- C++ -*-
+/**
+ *  @brief Boundary extraction algorithm for meshes
+ *
+ *  Copyright 2005-2023 Airbus-EDF-IMACS-ONERA-Phimeca
+ *
+ *  This library is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This library is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Lesser General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public License
+ *  along with this library.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+#ifndef OPENTURNS_MESHBOUNDARYFACTORY_HXX
+#define OPENTURNS_MESHBOUNDARYFACTORY_HXX
+
+#include "openturns/LevelSet.hxx"
+#include "openturns/Interval.hxx"
+#include "openturns/Mesh.hxx"
+#include "openturns/OptimizationAlgorithm.hxx"
+#include "openturns/AbdoRackwitz.hxx"
+
+BEGIN_NAMESPACE_OPENTURNS
+
+/**
+ * @class MeshBoundaryFactory
+ */
+class OT_API MeshBoundaryFactory
+  : public PersistentObject
+{
+  CLASSNAME
+public:
+
+  /** Default constructor */
+  MeshBoundaryFactory();
+
+  /** Virtual constructor */
+  MeshBoundaryFactory * clone() const override;
+
+  /** String converter */
+  String __repr__() const override;
+
+  /** String converter */
+  String __str__(const String & offset = "") const override;
+
+  /* Here is the interface that all derived class must implement */
+  /** Build the boundary of the given mesh */
+  virtual Mesh build(const Mesh & mesh,
+                     const Scalar offset = 0.0) const;
+
+protected:
+
+private:
+
+  /* Discretization in each dimension */
+  Indices discretization_;
+
+  /* Optimization solver used to project the vertices */
+  mutable OptimizationAlgorithm solver_;
+
+}; /* class MeshBoundaryFactory */
+
+END_NAMESPACE_OPENTURNS
+
+#endif /* OPENTURNS_MESHBOUNDARYFACTORY_HXX */
+

lib/src/Base/Geom/openturns/OTGeom.hxx

@@ -29,6 +29,7 @@
 #include "openturns/LevelSetMesher.hxx"
 #include "openturns/Mesh.hxx"
 #include "openturns/MeshDomain.hxx"
+#include "openturns/MeshBoundaryFactory.hxx"
 #include "openturns/RegularGrid.hxx"
 #include "openturns/DomainComplement.hxx"
 #include "openturns/DomainIntersection.hxx"