From 3d59f750151734a369695081f81083f012bbd288 Mon Sep 17 00:00:00 2001 From: regislebrun Date: Wed, 13 Sep 2023 16:14:47 +0200 Subject: [PATCH] Added BoundaryMesher This class allows to build the mesh of the boundary of a given mesh, with an optional offset. --- lib/src/Base/Geom/BoundaryMesher.cxx | 325 ++++++++++++++++++ lib/src/Base/Geom/CMakeLists.txt | 2 + .../Base/Geom/openturns/BoundaryMesher.hxx | 73 ++++ lib/src/Base/Geom/openturns/OTGeom.hxx | 1 + python/doc/pyplots/BoundaryMesher.py | 49 +++ python/doc/user_manual/base_objects.rst | 1 + python/src/BoundaryMesher.i | 11 + python/src/BoundaryMesher_doc.i.in | 56 +++ python/src/CMakeLists.txt | 1 + python/src/experimental_module.i | 4 + python/test/CMakeLists.txt | 1 + python/test/t_BoundaryMesher_std.py | 31 ++ 12 files changed, 555 insertions(+) create mode 100644 lib/src/Base/Geom/BoundaryMesher.cxx create mode 100644 lib/src/Base/Geom/openturns/BoundaryMesher.hxx create mode 100644 python/doc/pyplots/BoundaryMesher.py create mode 100644 python/src/BoundaryMesher.i create mode 100644 python/src/BoundaryMesher_doc.i.in create mode 100644 python/test/t_BoundaryMesher_std.py diff --git a/lib/src/Base/Geom/BoundaryMesher.cxx b/lib/src/Base/Geom/BoundaryMesher.cxx new file mode 100644 index 00000000000..2eda82853a6 --- /dev/null +++ b/lib/src/Base/Geom/BoundaryMesher.cxx @@ -0,0 +1,325 @@ +// -*- 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 . + * + */ +#include + +#include "openturns/PersistentObjectFactory.hxx" +#include "openturns/BoundaryMesher.hxx" +#include "openturns/TBBImplementation.hxx" + +BEGIN_NAMESPACE_OPENTURNS + +CLASSNAMEINIT(BoundaryMesher) +static const Factory Factory_BoundaryMesher; + + +/* Default constructor */ +BoundaryMesher::BoundaryMesher() + : PersistentObject() +{ + // Nothing to do +} + +/* Virtual constructor */ +BoundaryMesher * BoundaryMesher::clone() const +{ + return new BoundaryMesher(*this); +} + +/* String converter */ +String BoundaryMesher::__repr__() const +{ + OSS oss(true); + oss << "class=" << BoundaryMesher::GetClassName(); + return oss; +} + +/* String converter */ +String BoundaryMesher::__str__(const String & ) const +{ + return __repr__(); +} + +/* Here is the interface that all derived class must implement */ + + +namespace +{ +struct FaceHasher +{ + // This hash function is taken from + // https://stackoverflow.com/questions/20511347/a-good-hash-function-for-a-vector + // It is adapted to the fact that only the (size-1) first components have to be + // considered. + 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 + +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 + +/* The faces defining the boundary of a mesh of simplices are the faces of + simplices not shared by two simplices. + The faces of a simplex are obtained by removing one of its vertices, so in + dimension d a simplex has (d+1) faces. + In order to compare two faces, which are made of vertices indices, we sort the + indices. + As we want to keep track of the removed vertex in order to know on which side + of the face is the interior of the domain, we store a face as a (d+1) list of + integers, the d first ones being the indices of the face vertices, the last + one is the index of the removed vertex. A consequence is that we must provide + a dedicated hash function and comparison operator to the std::unordered_map + The whole algorithm is made of 3 main steps: + 1) Detect the unique faces + 2) Select the corresponding vertices and compact the corresponding indices + 3) (optional) If one ask for a thick boundary, compute a normal vector for + each face of the boundary and create a new vertex on the correct side of + the face (inside of the domain if offset<0, outside of the domain if + offset>0) +*/ +Mesh BoundaryMesher::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 faces; + faces.reserve(simplices.getSize() * (dimension + 1)); + Indices face(dimension + 1); + std::unordered_map::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::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; + { + const UnsignedInteger nbVertices = vertices.getSize(); + Indices oldToNewIndices(nbVertices, nbVertices); + 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]; + // 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 + 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); + + // 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) + { + const ComputeNormalPolicy policy(vertices, boundaryFaces, boundaryVertices.getSize() - boundaryVerticesReserve, offset, boundaryVertices, boundarySimplices); + TBBImplementation::ParallelFor(0, boundaryFaces.getSize(), policy); + } // offset != 0 + // Return the boundary mesh + return Mesh(boundaryVertices, boundarySimplices, false); +} + +END_NAMESPACE_OPENTURNS diff --git a/lib/src/Base/Geom/CMakeLists.txt b/lib/src/Base/Geom/CMakeLists.txt index 5cd266e8586..cf8d27ae5d4 100644 --- a/lib/src/Base/Geom/CMakeLists.txt +++ b/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 (BoundaryMesher.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 (BoundaryMesher.hxx) ot_install_header_file (RegularGrid.hxx) ot_install_header_file (DomainComplement.hxx) ot_install_header_file (DomainIntersection.hxx) diff --git a/lib/src/Base/Geom/openturns/BoundaryMesher.hxx b/lib/src/Base/Geom/openturns/BoundaryMesher.hxx new file mode 100644 index 00000000000..68a171057c8 --- /dev/null +++ b/lib/src/Base/Geom/openturns/BoundaryMesher.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 . + * + */ +#ifndef OPENTURNS_BOUNDARYMESHER_HXX +#define OPENTURNS_BOUNDARYMESHER_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 BoundaryMesher + */ +class OT_API BoundaryMesher + : public PersistentObject +{ + CLASSNAME +public: + + /** Default constructor */ + BoundaryMesher(); + + /** Virtual constructor */ + BoundaryMesher * 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 BoundaryMesher */ + +END_NAMESPACE_OPENTURNS + +#endif /* OPENTURNS_BOUNDARYMESHER_HXX */ + diff --git a/lib/src/Base/Geom/openturns/OTGeom.hxx b/lib/src/Base/Geom/openturns/OTGeom.hxx index 033368ddfce..3619a3150ba 100644 --- a/lib/src/Base/Geom/openturns/OTGeom.hxx +++ b/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/BoundaryMesher.hxx" #include "openturns/RegularGrid.hxx" #include "openturns/DomainComplement.hxx" #include "openturns/DomainIntersection.hxx" diff --git a/python/doc/pyplots/BoundaryMesher.py b/python/doc/pyplots/BoundaryMesher.py new file mode 100644 index 00000000000..2c35f3be02c --- /dev/null +++ b/python/doc/pyplots/BoundaryMesher.py @@ -0,0 +1,49 @@ +import openturns as ot +from openturns.viewer import View +import openturns.experimental as otexp +import matplotlib.pyplot as plt + +# Define the vertices of the mesh +vertices = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [1.5, 1.0], [3.0, 1.5], [0.5, 1.5]] +# Define the simplices of the mesh +simplices = [[0, 1, 2], [1, 2, 3], [2, 3, 4], [2, 4, 5], [0, 2, 5]] +# Create the Mesh +mesh2D = ot.Mesh(vertices, simplices) +# Build the Mesh boundary +mesh2DBoundary = otexp.BoundaryMesher().build(mesh2D) +# Build a thick outside Mesh boundary +mesh2DBoundaryOutside = otexp.BoundaryMesher().build(mesh2D, 0.05) +# Build a thick inside Mesh boundary +mesh2DBoundaryInside = otexp.BoundaryMesher().build(mesh2D, -0.05) +# Create a Graph +graph = ot.Graph("", "", "", True, "bottomright") +graph.add(mesh2D.draw()) +# Then, draw it +fig = plt.figure(figsize=(16, 4)) +axis = fig.add_subplot(141) +View(graph, figure=fig, axes=[axis], add_legend=True) +axis.set_xlim(auto=True) + +# Create a Graph +graph = ot.Graph("", "", "", True, "bottomright") +graph.add(mesh2DBoundary.draw()) +# Then, draw it +axis = fig.add_subplot(142) +View(graph, figure=fig, axes=[axis], add_legend=True) +axis.set_xlim(auto=True) + +# Create a Graph +graph = ot.Graph("", "", "", True, "bottomright") +graph.add(mesh2DBoundaryInside.draw()) +# Then, draw it +axis = fig.add_subplot(143) +View(graph, figure=fig, axes=[axis], add_legend=True) +axis.set_xlim(auto=True) + +# Create a Graph +graph = ot.Graph("A 2D mesh with its boundary, its inside thick boundary and its outside thick boundary", "", "", True, "bottomright") +graph.add(mesh2DBoundaryOutside.draw()) +# Then, draw it +axis = fig.add_subplot(144) +View(graph, figure=fig, axes=[axis], add_legend=True) +axis.set_xlim(auto=True) diff --git a/python/doc/user_manual/base_objects.rst b/python/doc/user_manual/base_objects.rst index 4eb013950bb..e2f5ac1aae6 100644 --- a/python/doc/user_manual/base_objects.rst +++ b/python/doc/user_manual/base_objects.rst @@ -48,6 +48,7 @@ Domains :template: classWithPlot.rst_t Mesh + experimental.BoundaryMesher Matrices ======== diff --git a/python/src/BoundaryMesher.i b/python/src/BoundaryMesher.i new file mode 100644 index 00000000000..dafc45d05c9 --- /dev/null +++ b/python/src/BoundaryMesher.i @@ -0,0 +1,11 @@ +// SWIG file BoundaryMesher.i + +%{ +#include "openturns/BoundaryMesher.hxx" +%} + +%include BoundaryMesher_doc.i + +%include openturns/BoundaryMesher.hxx + +namespace OT {%extend BoundaryMesher {BoundaryMesher(const BoundaryMesher & other){return new OT::BoundaryMesher(other);}}} diff --git a/python/src/BoundaryMesher_doc.i.in b/python/src/BoundaryMesher_doc.i.in new file mode 100644 index 00000000000..abeb8ae79a8 --- /dev/null +++ b/python/src/BoundaryMesher_doc.i.in @@ -0,0 +1,56 @@ +%feature("docstring") OT::BoundaryMesher +"Creation of the boundary mesh of a given mesh.. + +.. warning:: + This class is experimental and likely to be modified in future releases. + To use it, import the ``openturns.experimental`` submodule. + +Available constructor: + BoundaryMesher() + +Notes +----- +The boundary extraction is based on the remark that a :math:`(d-1)` face of a +:math:`d` dimensional simplex of a mesh is a boundary face if and only if this +face is not shared by any other simplex of the mesh. By convention, a face is +represented as a *flat* simplex, by repeating the :math:`d` th vertex so the +two last indices of the simplex are equal. +A *thick* version of the boundary can also be built, where each face is replaced +by a simplex with a last vertex at distance *offset* from the face hyperplane. +If *offset* is positive, this new vertex is in the open half-space not containing +the simplex associated to the face, otherwise it is in the open half-space +containing the simplex. + +Examples +-------- +Create a mesh: + +>>> import openturns as ot +>>> import openturns.experimental as otexp +>>> mesher = ot.LevelSetMesher([5, 10]) +>>> level = 1.0 +>>> function = ot.SymbolicFunction(['x0', 'x1'], ['x0^2+x1^2']) +>>> levelSet = ot.LevelSet(function, ot.LessOrEqual(), level) +>>> mesh = mesher.build(levelSet, ot.Interval([-2.0]*2, [2.0]*2)) +>>> boundaryFactory = otexp.BoundaryMesher() +>>> boundaryMesh = boundaryFactory.build(mesh) +>>> thickBoundaryMeshOutside = boundaryFactory.build(mesh, 1e-5) +>>> thickBoundaryMeshInside = boundaryFactory.build(mesh, -1e-5) +" + +// --------------------------------------------------------------------- + +%feature("docstring") OT::BoundaryMesher::build +"Build the boundary of the given mesh. + +Parameters +---------- +mesh : :class:`~openturns.Mesh` + The mesh from which the boundary is extracted. +offset : float + The thickness of the boundary to be generated, see Notes. + +Returns +------- +mesh : :class:`~openturns.Mesh` + The mesh built." diff --git a/python/src/CMakeLists.txt b/python/src/CMakeLists.txt index f16e0b86e26..f5e390c8648 100644 --- a/python/src/CMakeLists.txt +++ b/python/src/CMakeLists.txt @@ -1030,6 +1030,7 @@ ot_add_python_module(experimental experimental_module.i IntegrationExpansion.i IntegrationExpansion_doc.i.in UniformOverMesh.i UniformOverMesh_doc.i.in GeneralizedExtremeValueValidation.i GeneralizedExtremeValueValidation_doc.i.in + BoundaryMesher.i BoundaryMesher_doc.i.in ) set (OPENTURNS_PYTHON_MODULES ${OPENTURNS_PYTHON_MODULES} PARENT_SCOPE) # for the docstring test diff --git a/python/src/experimental_module.i b/python/src/experimental_module.i index 9fb3ab5b4de..116f829d28e 100644 --- a/python/src/experimental_module.i +++ b/python/src/experimental_module.i @@ -27,6 +27,9 @@ %include BaseFuncCollection.i %import base_module.i +/* Base/Geom */ +%include BoundaryMesher.i + /* Uncertainty/Model */ /* Uncertainty/Distribution */ %import model_copula_module.i @@ -44,6 +47,7 @@ %import metamodel_module.i %include simulation_module.i + /* Uncertainty/Algorithm/Metamodel */ %include UserDefinedMetropolisHastings.i %include FieldFunctionalChaosResult.i diff --git a/python/test/CMakeLists.txt b/python/test/CMakeLists.txt index d92f032a9fe..54d88c7c820 100644 --- a/python/test/CMakeLists.txt +++ b/python/test/CMakeLists.txt @@ -261,6 +261,7 @@ ot_pyinstallcheck_test (LevelSet_std) ot_pyinstallcheck_test (LevelSetMesher_std) ot_pyinstallcheck_test (Mesh_std) ot_pyinstallcheck_test (Mesh_draw IGNOREOUT) +ot_pyinstallcheck_test (BoundaryMesher_std IGNOREOUT) ot_pyinstallcheck_test (DomainComplement_std) ot_pyinstallcheck_test (DomainDifference_std) ot_pyinstallcheck_test (DomainIntersection_std) diff --git a/python/test/t_BoundaryMesher_std.py b/python/test/t_BoundaryMesher_std.py new file mode 100644 index 00000000000..ed9893e3fb8 --- /dev/null +++ b/python/test/t_BoundaryMesher_std.py @@ -0,0 +1,31 @@ +import openturns as ot +import openturns.testing as ott +import openturns.experimental as otexp + +def compute_length(mesh2D): + length = 0.0 + for ind in mesh2D.getSimplices(): + length += (mesh2D.getVertex(ind[0]) - mesh2D.getVertex(ind[1])).norm() + return length + + +# Define the vertices of the mesh +vertices = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [1.5, 1.0], [3.0, 1.5], [0.5, 1.5]] +# Define the simplices of the mesh +simplices = [[0, 1, 2], [1, 2, 3], [2, 3, 4], [2, 4, 5], [0, 2, 5]] +# Create the Mesh +mesh2D = ot.Mesh(vertices, simplices) +# Build the Mesh boundary +mesh2DBoundary = otexp.BoundaryMesher().build(mesh2D) +volume = mesh2DBoundary.getVolume() +ott.assert_almost_equal(volume, 0.0) +length = compute_length(mesh2DBoundary) +ott.assert_almost_equal(length, 7.780311648918275) +# Build a thick outside Mesh boundary +mesh2DBoundaryOutside = otexp.BoundaryMesher().build(mesh2D, 0.05) +volume = mesh2DBoundaryOutside.getVolume() +ott.assert_almost_equal(volume, 0.194508) +# Build a thick inside Mesh boundary +mesh2DBoundaryInside = otexp.BoundaryMesher().build(mesh2D, -0.05) +volume = mesh2DBoundaryInside.getVolume() +ott.assert_almost_equal(volume, 0.194508)