OPTION(MEDCOUPLING_ENABLE_PYTHON "Build PYTHON bindings." ON)
OPTION(MEDCOUPLING_ENABLE_PARTITIONER "Build MEDPartitioner." ON)
OPTION(MEDCOUPLING_ENABLE_RENUMBER "Build Renumber." ON)
+OPTION(MEDCOUPLING_ENABLE_SHAPERECOGN "Build ShapeRecogn" OFF)
OPTION(MEDCOUPLING_WITH_FILE_EXAMPLES "Install examples of files containing meshes and fields of different formats." ON)
OPTION(MEDCOUPLING_USE_MPI "(Use MPI containers) - For MED this triggers the build of ParaMEDMEM." OFF)
OPTION(MEDCOUPLING_BUILD_TESTS "Build MEDCoupling C++ tests." ON)
SALOME_LOG_OPTIONAL_PACKAGE(Boost MEDCOUPLING_ENABLE_RENUMBER)
ENDIF(MEDCOUPLING_ENABLE_RENUMBER)
+IF(MEDCOUPLING_ENABLE_SHAPERECOGN)
+ FIND_PACKAGE(BLAS REQUIRED)
+ FIND_PACKAGE(LAPACK REQUIRED)
+ FIND_LIBRARY(LAPACKE_LIB NAMES lapacke REQUIRED)
+ SET(LAPACK_LIBRARIES ${LAPACKE_LIB} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
+ FIND_PATH(LAPACKE_INCLUDE_DIRS NAMES lapacke.h HINTS ${LAPACK_LIBRARIES})
+ IF(LAPACK_FOUND)
+ MESSAGE(STATUS "Lapacke libraries: ${LAPACK_LIBRARIES}")
+ MESSAGE(STATUS "Lapacke include dirs: ${LAPACKE_INCLUDE_DIRS}")
+ ELSE()
+ MESSAGE(FATAL_ERROR "Error in Lapacke detection ! lapacke not found !")
+ ENDIF(LAPACK_FOUND)
+ SALOME_LOG_OPTIONAL_PACKAGE(Lapack MEDCOUPLING_ENABLE_SHAPERECOGN)
+ENDIF(MEDCOUPLING_ENABLE_SHAPERECOGN)
+
IF(MEDCOUPLING_ENABLE_PYTHON)
FIND_PACKAGE(SalomePythonInterp)
FIND_PACKAGE(SalomePythonLibs)
ENDIF()
ENDIF()
ENDIF()
+
+IF(MEDCOUPLING_ENABLE_SHAPERECOGN)
+ LIST(APPEND _${PROJECT_NAME}_exposed_targets shaperecogn)
+ENDIF(MEDCOUPLING_ENABLE_SHAPERECOGN)
# Add all targets to the build-tree export set
test_MED_MAIL.sauv
castem17_result_xdr.sauv
castem17_result_ascii.sauv
+ ShapeRecognPlane.med
+ ShapeRecognCylinder.med
+ ShapeRecognCone.med
+ ShapeRecognSphere.med
+ ShapeRecognTorus.med
)
SET(MED_RESOURCES_FILES ${MED_test_fig_files})
ENDIF(MEDCOUPLING_BUILD_TESTS)
ENDIF(MEDCOUPLING_USE_MPI)
+IF(MEDCOUPLING_ENABLE_SHAPERECOGN)
+ ADD_SUBDIRECTORY(ShapeRecogn)
+ENDIF(MEDCOUPLING_ENABLE_SHAPERECOGN)
+
# Application tests
configure_file(CTestTestfileInstall.cmake.in "CTestTestfileST.cmake" @ONLY)
install(FILES ${CMAKE_CURRENT_BINARY_DIR}/CTestTestfileST.cmake DESTINATION ${MEDCOUPLING_INSTALL_TESTS} RENAME CTestTestfile.cmake)
--- /dev/null
+#include "Areas.hxx"
+#include "MathOps.hxx"
+
+#include <random>
+#include <cblas.h>
+
+using namespace MEDCoupling;
+
+Areas::Areas(const Nodes *nodes) : nodes(nodes), areaIdByNodes(nodes->getNbNodes(), -1)
+{
+}
+
+mcIdType Areas::addArea(PrimitiveType primitive)
+{
+ Area area;
+ area.primitive = primitive;
+ areas.push_back(area);
+ return areas.size() - 1;
+}
+
+void Areas::cleanArea(mcIdType areaId)
+{
+ cleanArea(areaId, -1);
+}
+
+void Areas::cancelArea(mcIdType areaId, PrimitiveType primitive)
+{
+ cleanArea(areaId, -2);
+ Area &area = areas[areaId];
+ area.primitive = primitive;
+}
+
+void Areas::removeArea(mcIdType areaId)
+{
+ for (size_t nodeId = 0; nodeId < areaIdByNodes.size(); ++nodeId)
+ {
+ if (areaIdByNodes[nodeId] == areaId)
+ areaIdByNodes[nodeId] = -1;
+ else if (areaIdByNodes[nodeId] > areaId)
+ areaIdByNodes[nodeId] -= 1;
+ }
+ areas.erase(areas.begin() + areaId);
+}
+
+void Areas::addNode(mcIdType areaId, mcIdType nodeId)
+{
+ removeNode(nodeId);
+ areaIdByNodes[nodeId] = areaId;
+ Area &area = areas[areaId];
+ area.nodeIds.push_back(nodeId);
+ size_t nbNodes = area.nodeIds.size();
+ area.k1 = ((double)(nbNodes - 1) * area.k1 + nodes->getK1(nodeId)) / (double)nbNodes;
+ area.k2 = ((double)(nbNodes - 1) * area.k2 + nodes->getK2(nodeId)) / (double)nbNodes;
+ area.adimK1 = ((double)(nbNodes - 1) * area.adimK1 + nodes->getAdimK1(nodeId)) / (double)nbNodes;
+ area.adimK2 = ((double)(nbNodes - 1) * area.adimK2 + nodes->getAdimK2(nodeId)) / (double)nbNodes;
+ area.adimKdiff0 = ((double)(nbNodes - 1) * area.adimKdiff0 + nodes->getAdimKdiff0(nodeId)) / (double)nbNodes;
+}
+
+void Areas::cleanInvalidNodeAreas()
+{
+ for (mcIdType nodeId = 0; nodeId < (mcIdType)areaIdByNodes.size(); ++nodeId)
+ {
+ if (areaIdByNodes[nodeId] <= -2)
+ areaIdByNodes[nodeId] = -1;
+ }
+}
+
+mcIdType Areas::getAreaId(mcIdType nodeId) const
+{
+ return areaIdByNodes[nodeId];
+}
+
+bool Areas::isEmpty(mcIdType areaId) const
+{
+ return areas[areaId].nodeIds.empty();
+}
+
+size_t Areas::getNumberOfAreas() const
+{
+ return areas.size();
+}
+
+size_t Areas::getNumberOfNodes(mcIdType areaId) const
+{
+ return areas[areaId].nodeIds.size();
+}
+
+bool Areas::isNodeCompatible(mcIdType areaId, mcIdType nodeId) const
+{
+ PrimitiveType areaType = getPrimitiveType(areaId);
+ PrimitiveType nodeType = nodes->getPrimitiveType(nodeId);
+ return (
+ (areaType != PrimitiveType::Cylinder ||
+ nodeType != PrimitiveType::Plane) &&
+ (areaType != PrimitiveType::Sphere ||
+ nodeType != PrimitiveType::Plane));
+}
+
+PrimitiveType Areas::getPrimitiveType(mcIdType areaId) const
+{
+ const Area &area = areas[areaId];
+ if (area.primitive != PrimitiveType::Unknown)
+ return area.primitive;
+ else if (area.nodeIds.empty())
+ return PrimitiveType::Unknown;
+ else
+ return nodes->getPrimitiveType(area.nodeIds[0]);
+}
+
+std::string Areas::getPrimitiveTypeName(mcIdType areaId) const
+{
+ return convertPrimitiveToString(getPrimitiveType(areaId));
+}
+
+int Areas::getPrimitiveTypeInt(mcIdType areaId) const
+{
+ return convertPrimitiveToInt(getPrimitiveType(areaId));
+}
+
+const std::vector<mcIdType> &Areas::getNodeIds(mcIdType areaId) const
+{
+ return areas[areaId].nodeIds;
+}
+
+double Areas::getAdimK1(mcIdType areaId) const
+{
+
+ return areas[areaId].adimK1;
+}
+
+double Areas::getAdimK2(mcIdType areaId) const
+{
+ return areas[areaId].adimK2;
+}
+
+double Areas::getAdimKdiff0(mcIdType areaId) const
+{
+ return areas[areaId].adimKdiff0;
+}
+
+void Areas::computeProperties(mcIdType areaId)
+{
+ switch (getPrimitiveType(areaId))
+ {
+ case PrimitiveType::Plane:
+ computePlaneProperties(areaId);
+ break;
+ case PrimitiveType::Sphere:
+ computeSphereProperties(areaId);
+ break;
+ case PrimitiveType::Cylinder:
+ computeCylinderProperties(areaId);
+ break;
+ case PrimitiveType::Cone:
+ computeConeProperties(areaId);
+ break;
+ case PrimitiveType::Torus:
+ computeTorusProperties(areaId);
+ break;
+ case PrimitiveType::Unknown:
+ default:
+ break;
+ }
+}
+
+double Areas::getMinorRadius(mcIdType areaId) const
+{
+ const Area &area = areas[areaId];
+ return area.minorRadius;
+}
+
+double Areas::getRadius(mcIdType areaId) const
+{
+ const Area &area = areas[areaId];
+ return area.radius;
+}
+
+double Areas::getAngle(mcIdType areaId) const
+{
+ return areas[areaId].angle;
+}
+
+const std::array<double, 3> &Areas::getNormal(mcIdType areaId) const
+{
+ return areas[areaId].normal;
+}
+
+const std::array<double, 3> &Areas::getCenter(mcIdType areaId) const
+{
+ const Area &area = areas[areaId];
+ return area.center;
+}
+
+const std::array<double, 3> &Areas::getAxis(mcIdType areaId) const
+{
+ return areas[areaId].axis;
+}
+
+const std::array<double, 3> &Areas::getAxisPoint(mcIdType areaId) const
+{
+ return areas[areaId].axisPoint;
+}
+
+const std::array<double, 3> &Areas::getApex(mcIdType areaId) const
+{
+ return areas[areaId].apex;
+}
+
+std::array<double, 3> Areas::getAffinePoint(mcIdType areaId) const
+{
+ const Area &area = areas[areaId];
+ if (area.nodeIds.empty())
+ return {0.0, 0.0, 0.0};
+ else
+ return nodes->getCoordinates(area.nodeIds[0]);
+}
+
+void Areas::cleanArea(mcIdType areaId, mcIdType newAreaId = -1)
+{
+ Area &area = areas[areaId];
+ for (mcIdType nodeId : area.nodeIds)
+ areaIdByNodes[nodeId] = newAreaId;
+ area.primitive = PrimitiveType::Unknown;
+ area.k1 = 0.0;
+ area.k2 = 0.0;
+ area.adimK1 = 0.0;
+ area.adimK2 = 0.0;
+ area.adimKdiff0 = 0.0;
+ area.minorRadius = 0.0;
+ area.radius = 0.0;
+ area.angle = 0.0;
+ area.normal = {0.0, 0.0, 0.0};
+ area.center = {0.0, 0.0, 0.0};
+ area.axis = {0.0, 0.0, 0.0};
+ area.axisPoint = {0.0, 0.0, 0.0};
+ area.apex = {0.0, 0.0, 0.0};
+ area.nodeIds.clear();
+}
+
+void Areas::removeNode(mcIdType nodeId)
+{
+ mcIdType areaId = areaIdByNodes[nodeId];
+ if (areaId > -1)
+ {
+ Area &area = areas[areaId];
+ area.nodeIds.erase(
+ std::remove(
+ area.nodeIds.begin(),
+ area.nodeIds.end(),
+ nodeId),
+ area.nodeIds.end());
+ areaIdByNodes[nodeId] = -1;
+ // TODO: Update the parameters of the area ?
+ }
+}
+
+void Areas::computePlaneProperties(mcIdType areaId)
+{
+ Area &area = areas[areaId];
+ const std::vector<double> &normals = nodes->getNormals();
+ mcIdType nbNodes = area.nodeIds.size();
+ area.normal = {0.0, 0.0, 0.0};
+ for (mcIdType nodeId : area.nodeIds)
+ {
+ for (size_t i = 0; i < 3; ++i)
+ area.normal[i] += normals[3 * nodeId + i];
+ }
+ for (size_t i = 0; i < 3; ++i)
+ area.normal[i] /= (double)nbNodes;
+}
+
+void Areas::computeSphereProperties(mcIdType areaId)
+{
+ Area &area = areas[areaId];
+ const std::vector<double> &normals = nodes->getNormals();
+ area.radius = (2 / (area.k1 + area.k2));
+ std::array<double, 3> center{0.0, 0.0, 0.0};
+ if (!area.nodeIds.empty())
+ {
+ size_t nbNodes = area.nodeIds.size();
+ for (mcIdType nodeId : area.nodeIds)
+ {
+ std::array<double, 3> nodeCoords = nodes->getCoordinates(nodeId);
+ for (size_t i = 0; i < 3; ++i)
+ center[i] += nodeCoords[i] - area.radius * normals[3 * nodeId + i];
+ }
+ for (size_t i = 0; i < 3; ++i)
+ center[i] /= (double)nbNodes;
+ }
+ area.center = center;
+ area.radius = fabs(area.radius);
+}
+
+void Areas::computeCylinderProperties(mcIdType areaId)
+{
+ Area &area = areas[areaId];
+ size_t nbNodes = area.nodeIds.size();
+ // Project the nodes to the central axis of the cylinder
+ std::vector<double> projectedNodes(3 * nbNodes, 0.0);
+ area.radius = 0;
+ area.axisPoint.fill(0.0);
+ for (size_t i = 0; i < nbNodes; ++i)
+ {
+ mcIdType nodeId = area.nodeIds[i];
+ std::array<double, 3> nodeCoords = nodes->getCoordinates(nodeId);
+ double approxRadius = 1.0 / nodes->getKdiff0(nodeId);
+ for (size_t j = 0; j < 3; ++j)
+ {
+ projectedNodes[3 * i + j] =
+ nodeCoords[j] - approxRadius * nodes->getNormals()[3 * nodeId + j];
+ area.axisPoint[j] += projectedNodes[3 * i + j];
+ }
+ area.radius += approxRadius;
+ }
+ // Axis point is the mean of the projected nodes
+ for (size_t i = 0; i < 3; ++i)
+ area.axisPoint[i] /= (double)nbNodes;
+ // Radius of the cylinder is the mean of the approximate radius of each node
+ area.radius = fabs(area.radius / (double)nbNodes);
+ // Compute the axis of the cylinder
+ area.axis = MathOps::computePCAFirstAxis(projectedNodes);
+}
+
+void Areas::computeConeProperties(mcIdType areaId)
+{
+ Area &area = areas[areaId];
+ size_t nbNodes = area.nodeIds.size();
+ // Project the nodes to the central axis of the cone
+ std::vector<double> projectedNodes(3 * nbNodes, 0.0);
+ std::vector<double> radiusNodes(nbNodes, 0.0);
+ area.axisPoint.fill(0.0);
+ for (size_t i = 0; i < nbNodes; ++i)
+ {
+ mcIdType nodeId = area.nodeIds[i];
+ std::array<double, 3> nodeCoords = nodes->getCoordinates(nodeId);
+ radiusNodes[i] = 1.0 / nodes->getKdiff0(nodeId);
+ for (size_t j = 0; j < 3; ++j)
+ {
+ projectedNodes[3 * i + j] =
+ nodeCoords[j] - radiusNodes[i] * nodes->getNormals()[3 * nodeId + j];
+ area.axisPoint[j] += projectedNodes[3 * i + j];
+ }
+ }
+ // Axis point is the mean of the projected nodes
+ for (size_t i = 0; i < 3; ++i)
+ area.axisPoint[i] /= (double)nbNodes;
+ // Compute the axis of the cone
+ area.axis = MathOps::computePCAFirstAxis(projectedNodes);
+ double normAxis = MathOps::computeNorm(area.axis);
+ for (size_t i = 0; i < 3; ++i)
+ area.axis[i] /= normAxis;
+ // Compute the angle of the cone
+ const std::vector<double> &weakDirections = nodes->getWeakDirections();
+ std::vector<double> weakDirectionNodes(3 * nbNodes, 0.0);
+ for (size_t i = 0; i < nbNodes; ++i)
+ {
+ mcIdType nodeId = area.nodeIds[i];
+ for (size_t j = 0; j < 3; ++j)
+ weakDirectionNodes[3 * i + j] = weakDirections[3 * nodeId + j];
+ }
+ std::vector<double> angles = MathOps::computeAngles(
+ weakDirectionNodes, area.axis);
+ // Correct the angles > pi/2 with the supplementary angle
+ for (size_t i = 0; i < angles.size(); ++i)
+ {
+ if (angles[i] > M_PI_2)
+ angles[i] = M_PI - angles[i];
+ }
+ area.angle = MathOps::mean(angles);
+ // Compute the radius of the cone which is the mean of the approximate radius of each node
+ area.radius = MathOps::mean(radiusNodes);
+ // Select extrem nodes
+ double q1 = MathOps::computeQuantile(radiusNodes, 0.1);
+ double q2 = MathOps::computeQuantile(radiusNodes, 0.9);
+ std::vector<mcIdType> q1_indices;
+ for (mcIdType idx = 0; idx < (mcIdType)radiusNodes.size(); ++idx)
+ {
+ if (radiusNodes[idx] < q1)
+ q1_indices.push_back(idx);
+ }
+ std::vector<mcIdType> q2_indices;
+ for (mcIdType idx = 0; idx < (mcIdType)radiusNodes.size(); ++idx)
+ {
+ if (radiusNodes[idx] > q2)
+ q2_indices.push_back(idx);
+ }
+ std::random_device rd;
+ std::mt19937 g(rd());
+ std::shuffle(q2_indices.begin(), q2_indices.end(), g);
+ // Compute the height of the cone
+ // std::vector<double> heights(q1_indices.size(), 0.0);
+ // std::vector<double> distancesToApex(q1_indices.size(), 0.0);
+ std::array<double, 3> p{0.0, 0.0, 0.0};
+ size_t min_q_size = std::min<size_t>(q1_indices.size(), q2_indices.size());
+ for (size_t i = 0; i < min_q_size; ++i)
+ {
+ for (size_t j = 0; j < 3; ++j)
+ p[j] = (projectedNodes[3 * q1_indices[i] + j] - projectedNodes[3 * q2_indices[i] + j]);
+ double height = MathOps::computeNorm(p);
+ double distanceToApex = height * radiusNodes[q1_indices[i]] / (radiusNodes[q2_indices[i]] - radiusNodes[q1_indices[i]]);
+ double orientationAxis = MathOps::dot(p, area.axis);
+ for (size_t j = 0; j < 3; ++j)
+ {
+ area.apex[j] += projectedNodes[3 * q1_indices[i] + j];
+ if (orientationAxis >= 0)
+ area.apex[j] += distanceToApex * area.axis[j];
+ else
+ area.apex[j] -= distanceToApex * area.axis[j];
+ }
+ }
+ for (size_t j = 0; j < 3; ++j)
+ area.apex[j] /= (double)min_q_size;
+}
+
+void Areas::computeTorusProperties(mcIdType areaId)
+{
+ Area &area = areas[areaId];
+ size_t n = area.nodeIds.size();
+ if (n == 0)
+ return;
+ std::vector<double> areaNodesK1(n, 0.0);
+ std::vector<double> areaNodesK2(n, 0.0);
+ for (size_t i = 0; i < n; ++i)
+ {
+ areaNodesK1[i] = nodes->getK1(area.nodeIds[i]);
+ areaNodesK2[i] = nodes->getK2(area.nodeIds[i]);
+ }
+ double var1 = MathOps::computeVariance(areaNodesK1);
+ double var2 = MathOps::computeVariance(areaNodesK2);
+ double minorCurvature;
+ if (var1 > var2)
+ minorCurvature = MathOps::mean(areaNodesK2);
+ else
+ minorCurvature = MathOps::mean(areaNodesK1);
+ area.minorRadius = 1.0 / minorCurvature;
+ std::vector<double> majorRadiusNodes(3 * n, 0.0);
+ for (size_t i = 0; i < n; ++i)
+ {
+ std::array<double, 3> coords = nodes->getCoordinates(area.nodeIds[i]);
+ std::array<double, 3> normal = nodes->getNormal(area.nodeIds[i]);
+ for (size_t j = 0; j < 3; ++j)
+ majorRadiusNodes[3 * i + j] = coords[j] - normal[j] * area.minorRadius;
+ }
+ std::array<double, 3> meanMajorRadiusNodes = MathOps::meanCoordinates(majorRadiusNodes);
+ for (size_t i = 0; i < n; ++i)
+ {
+ for (size_t j = 0; j < 3; ++j)
+ majorRadiusNodes[3 * i + j] -= meanMajorRadiusNodes[j];
+ }
+ std::array<double, 3> normal = MathOps::computePCAThirdAxis(majorRadiusNodes);
+ std::array<double, 6> base2d = MathOps::computeBaseFromNormal(normal);
+ std::vector<double> projectedMajorRadiusNodes(2 * n, 0.0);
+ cblas_dgemm(
+ CBLAS_LAYOUT::CblasRowMajor,
+ CBLAS_TRANSPOSE::CblasNoTrans,
+ CBLAS_TRANSPOSE::CblasTrans,
+ (int)n, 2, 3,
+ 1.0,
+ majorRadiusNodes.data(), 3,
+ base2d.data(), 3,
+ 0.0, projectedMajorRadiusNodes.data(), 2);
+ std::vector<double> A(3 * n, 1.0);
+ std::vector<double> B(n, 0.0);
+ for (size_t i = 0; i < n; ++i)
+ {
+ for (size_t j = 0; j < 2; ++j)
+ A[3 * i + j] = projectedMajorRadiusNodes[2 * i + j];
+ B[i] = pow(projectedMajorRadiusNodes[2 * i], 2) +
+ pow(projectedMajorRadiusNodes[2 * i + 1], 2);
+ }
+ std::vector<double> fit = MathOps::lstsqRow(A, B);
+ double a = fit[0];
+ double b = fit[1];
+ double c = fit[2];
+ double xc = a / 2.0;
+ double yc = b / 2.0;
+ area.radius = sqrt(4.0 * c + pow(a, 2) + pow(b, 2)) / 2.0;
+ for (size_t i = 0; i < 3; ++i)
+ area.center[i] = xc * base2d[i] + yc * base2d[3 + i] + meanMajorRadiusNodes[i];
+}
+
+const std::vector<mcIdType> &Areas::getAreaIdByNodes() const
+{
+ return areaIdByNodes;
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __AREAS_HXX__
+#define __AREAS_HXX__
+
+#include "PrimitiveType.hxx"
+#include "Nodes.hxx"
+
+#include <array>
+
+namespace MEDCoupling
+{
+
+ struct Area
+ {
+ PrimitiveType primitive = PrimitiveType::Unknown;
+ double k1 = 0.0;
+ double k2 = 0.0;
+ double adimK1 = 0.0;
+ double adimK2 = 0.0;
+ double adimKdiff0 = 0.0;
+ double minorRadius = 0.0;
+ double radius = 0.0;
+ double angle = 0.0;
+ std::array<double, 3> normal{0.0, 0.0, 0.0};
+ std::array<double, 3> center{0.0, 0.0, 0.0};
+ std::array<double, 3> axis{0.0, 0.0, 0.0};
+ std::array<double, 3> axisPoint{0.0, 0.0, 0.0};
+ std::array<double, 3> apex{0.0, 0.0, 0.0};
+ std::vector<mcIdType> nodeIds;
+ };
+
+ class Areas
+ {
+ public:
+ Areas(const Nodes *nodes);
+
+ mcIdType addArea(PrimitiveType primitive = PrimitiveType::Unknown);
+ void cleanArea(mcIdType areaId);
+ void cancelArea(mcIdType areaId, PrimitiveType primitive = PrimitiveType::Unknown);
+ void removeArea(mcIdType areaId);
+ void addNode(mcIdType areaId, mcIdType nodeId);
+ void cleanInvalidNodeAreas();
+
+ mcIdType getAreaId(mcIdType nodeId) const;
+ const std::vector<mcIdType> &getAreaIdByNodes() const;
+
+ bool isEmpty(mcIdType areaId) const;
+ size_t getNumberOfAreas() const;
+ size_t getNumberOfNodes(mcIdType areaId) const;
+
+ bool isNodeCompatible(mcIdType areaId, mcIdType nodeId) const;
+
+ PrimitiveType getPrimitiveType(mcIdType areaId) const;
+ std::string getPrimitiveTypeName(mcIdType areaId) const;
+ int getPrimitiveTypeInt(mcIdType areaId) const;
+
+ const std::vector<mcIdType> &getNodeIds(mcIdType areaId) const;
+
+ double getAdimK1(mcIdType areaId) const;
+ double getAdimK2(mcIdType areaId) const;
+ double getAdimKdiff0(mcIdType areaId) const;
+
+ void computeProperties(mcIdType areaId);
+ double getMinorRadius(mcIdType areaId) const;
+ double getRadius(mcIdType areaId) const;
+ double getAngle(mcIdType areaId) const;
+ const std::array<double, 3> &getNormal(mcIdType areaId) const;
+ const std::array<double, 3> &getCenter(mcIdType areaId) const;
+ const std::array<double, 3> &getAxis(mcIdType areaId) const;
+ const std::array<double, 3> &getAxisPoint(mcIdType areaId) const;
+ const std::array<double, 3> &getApex(mcIdType areaId) const;
+ std::array<double, 3> getAffinePoint(mcIdType areaId) const;
+
+ private:
+ void cleanArea(mcIdType areaId, mcIdType newAreaId);
+ void removeNode(mcIdType nodeId);
+
+ void computePlaneProperties(mcIdType areaId);
+ void computeSphereProperties(mcIdType areaId);
+ void computeCylinderProperties(mcIdType areaId);
+ void computeConeProperties(mcIdType areaId);
+ void computeTorusProperties(mcIdType areaId);
+
+ std::vector<Area> areas;
+ const Nodes *nodes;
+ std::vector<mcIdType> areaIdByNodes;
+ };
+};
+
+#endif // __AREAS_HXX__
--- /dev/null
+#include "AreasBuilder.hxx"
+#include "MathOps.hxx"
+#include "ShapeRecongConstants.hxx"
+
+#include <unordered_set>
+
+using namespace MEDCoupling;
+
+AreasBuilder::AreasBuilder(const Nodes *nodes) : nodes(nodes), areas(new Areas(nodes))
+{
+ size_t nbNodes = nodes->getNbNodes();
+ threshold = std::max<size_t>(THRESHOLD_MIN_NB_NODES, nbNodes / THRESHOLD_MAX_NB_AREAS);
+}
+
+void AreasBuilder::build()
+{
+ explore();
+ expand();
+ rebuild();
+}
+
+void AreasBuilder::explore()
+{
+ exploreAreas();
+ areas->cleanInvalidNodeAreas();
+ filterHighPass();
+}
+
+void AreasBuilder::expand()
+{
+ expandAreas();
+ filterHighPass();
+}
+
+void AreasBuilder::rebuild()
+{
+ rebuildInvalidAreas();
+ filterHighPass();
+ expandAreasByType(PrimitiveType::Cone);
+ filterHighPass();
+}
+
+Areas *AreasBuilder::getAreas() const
+{
+ return areas;
+}
+
+void AreasBuilder::exploreAreas()
+{
+ // int nbNodesExplored = 0;
+ std::vector<bool> exploredNodeIds(nodes->getNbNodes(), false);
+ std::unordered_set<mcIdType> nodesToExplore;
+ // Reserve a set with the size of nodes to avoid reallocation for each insert/erase
+ // TODO: Improve the size ? Nb of Nodes is too much ?
+ nodesToExplore.reserve(nodes->getNbNodes());
+ mcIdType areaId = -1;
+ for (mcIdType nodeId = 0; nodeId < nodes->getNbNodes(); ++nodeId)
+ {
+ if (!exploredNodeIds[nodeId] &&
+ nodes->getPrimitiveType(nodeId) != PrimitiveType::Unknown)
+ {
+ exploredNodeIds[nodeId] = true;
+ if (areaId != -1 && areas->getNumberOfNodes(areaId) < threshold)
+ areas->cancelArea(areaId, nodes->getPrimitiveType(nodeId));
+ else
+ areaId = areas->addArea(nodes->getPrimitiveType(nodeId));
+ areas->addNode(areaId, nodeId);
+ const std::vector<mcIdType> neighbors = nodes->getNeighbors(nodeId);
+ for (mcIdType neighborId : neighbors)
+ {
+ if (nodes->getPrimitiveType(neighborId) == areas->getPrimitiveType(areaId) &&
+ areas->getAreaId(neighborId) <= -1)
+ nodesToExplore.insert(neighborId);
+ }
+ // Explore all the neighbors matching the area
+ while (!nodesToExplore.empty())
+ {
+ mcIdType neighborId = *nodesToExplore.begin();
+ nodesToExplore.erase(neighborId);
+ if (doesItMatch(areaId, neighborId))
+ {
+ exploredNodeIds[neighborId] = true;
+ areas->addNode(areaId, neighborId);
+ const std::vector<mcIdType> neighborsOfNeighbor = nodes->getNeighbors(neighborId);
+ for (mcIdType neighborIdOfNeighbor : neighborsOfNeighbor)
+ {
+ if (!exploredNodeIds[neighborIdOfNeighbor] &&
+ areas->getAreaId(neighborIdOfNeighbor) <= -1 &&
+ // Already in doesItMatch but avoid useless insertion
+ nodes->getPrimitiveType(neighborIdOfNeighbor) == areas->getPrimitiveType(areaId))
+ nodesToExplore.insert(neighborIdOfNeighbor);
+ }
+ }
+ }
+ }
+ // if (!exploredNodeIds[nodeId])
+ // nbNodesExplored += 1;
+ exploredNodeIds[nodeId] = true;
+ }
+}
+
+void AreasBuilder::expandAreas()
+{
+ // Expand by topological order
+ expandAreasByType(PrimitiveType::Plane);
+ expandAreasByType(PrimitiveType::Sphere);
+ expandAreasByType(PrimitiveType::Cylinder);
+ expandAreasByType(PrimitiveType::Cone);
+ expandAreasByType(PrimitiveType::Torus);
+}
+
+void AreasBuilder::expandAreasByType(PrimitiveType primitive)
+{
+ std::unordered_set<mcIdType> nodesToExplore;
+ // Reserve a set with the size of nodes to avoid reallocation for each insert/erase
+ // TODO: Improve the size ? Nb of Nodes is too much ?
+ nodesToExplore.reserve(nodes->getNbNodes());
+ for (mcIdType areaId = 0; areaId < (mcIdType)areas->getNumberOfAreas(); ++areaId)
+ {
+ if (areas->getPrimitiveType(areaId) == primitive)
+ {
+ std::vector<bool> exploredNodeIds(nodes->getNbNodes(), false);
+ areas->computeProperties(areaId);
+ const std::vector<mcIdType> &nodeIds = areas->getNodeIds(areaId);
+ for (mcIdType nodeId : nodeIds)
+ {
+ exploredNodeIds[nodeId] = true;
+ nodesToExplore.insert(nodeId);
+ }
+ while (!nodesToExplore.empty())
+ {
+ mcIdType nodeId = *nodesToExplore.begin();
+ nodesToExplore.erase(nodeId);
+ if (doesItBelong(areaId, nodeId))
+ {
+ // TODO: Is the properties need to be updated after adding a node ?
+ // It gives bad results for the cone and the cylinder
+ // mcIdType oldAreaId = areas->getAreaId(nodeId);
+ areas->addNode(areaId, nodeId);
+ // areas->computeProperties(areaId);
+ // areas->computeProperties(oldAreaId);
+ const std::vector<mcIdType> neighborIds = nodes->getNeighbors(nodeId);
+ for (mcIdType neighborId : neighborIds)
+ {
+ if (!exploredNodeIds[neighborId])
+ nodesToExplore.insert(neighborId);
+ }
+ }
+ exploredNodeIds[nodeId] = true;
+ }
+ }
+ }
+}
+
+void AreasBuilder::rebuildInvalidAreas()
+{
+ std::vector<mcIdType> exploredNodeIds(nodes->getNbNodes(), false);
+ std::vector<bool> isIinvalidNodes(nodes->getNbNodes(), false);
+ std::unordered_set<mcIdType> nodesToExplore;
+ // Reserve a set with the size of nodes to avoid reallocation for each insert/erase
+ // TODO: Improve the size ? Nb of Nodes is too much ?
+ nodesToExplore.reserve(nodes->getNbNodes());
+ for (mcIdType nodeId = 0; nodeId < nodes->getNbNodes(); ++nodeId)
+ isIinvalidNodes[nodeId] = isInvalidCylinderNode(nodeId);
+ for (mcIdType nodeId = 0; nodeId < nodes->getNbNodes(); ++nodeId)
+ {
+ if (isIinvalidNodes[nodeId] && !exploredNodeIds[nodeId])
+ {
+ mcIdType areaId = areas->addArea(PrimitiveType::Cone);
+ areas->addNode(areaId, nodeId);
+ const std::vector<mcIdType> neighbors = nodes->getNeighbors(nodeId);
+ for (mcIdType neighborId : neighbors)
+ nodesToExplore.insert(neighborId);
+ while (!nodesToExplore.empty())
+ {
+ mcIdType neighborId = *nodesToExplore.begin();
+ nodesToExplore.erase(neighborId);
+ if (
+ (areas->getAreaId(neighborId) == -1 ||
+ isIinvalidNodes[neighborId]))
+ {
+ exploredNodeIds[neighborId] = true;
+ areas->addNode(areaId, neighborId);
+ const std::vector<mcIdType> neighborsOfNeighbor = nodes->getNeighbors(neighborId);
+ for (mcIdType neighborIdOfNeighbor : neighborsOfNeighbor)
+ {
+ if (!exploredNodeIds[neighborIdOfNeighbor])
+ nodesToExplore.insert(neighborIdOfNeighbor);
+ }
+ }
+ }
+ }
+ }
+}
+
+void AreasBuilder::filterHighPass()
+{
+ mcIdType nbAreas = areas->getNumberOfAreas();
+ for (mcIdType areaId = (nbAreas - 1); areaId >= 0; --areaId)
+ {
+ if (areas->getNumberOfNodes(areaId) < threshold)
+ areas->removeArea(areaId);
+ }
+}
+
+bool AreasBuilder::doesItMatch(mcIdType areaId, mcIdType nodeId) const
+{
+ PrimitiveType areaType = areas->getPrimitiveType(areaId);
+ PrimitiveType neighborType = nodes->getPrimitiveType(nodeId);
+ bool isMatching = false;
+ if (areaType == neighborType)
+ {
+ switch (areaType)
+ {
+ case PrimitiveType::Plane:
+ case PrimitiveType::Torus:
+ isMatching = true;
+ break;
+ case PrimitiveType::Sphere:
+ {
+ double kmoy = fabs(
+ (areas->getAdimK1(areaId) + areas->getAdimK2(areaId)) / 2.0);
+ double nodeKmoy = fabs(
+ (nodes->getAdimK1(nodeId) + nodes->getAdimK2(nodeId)) / 2.0);
+ isMatching = fabs((nodeKmoy - kmoy)) < TOL_MATCH_SPHERE;
+ }
+ break;
+ case PrimitiveType::Cylinder:
+ isMatching = fabs(areas->getAdimKdiff0(areaId) -
+ nodes->getAdimKdiff0(nodeId)) < TOL_MATCH_CYLINDER;
+ break;
+ case PrimitiveType::Cone:
+ case PrimitiveType::Unknown:
+ default:
+ break;
+ }
+ }
+ return isMatching;
+}
+
+bool AreasBuilder::doesItBelong(mcIdType areaId, mcIdType nodeId) const
+{
+ bool isClose = false;
+ if (areas->isNodeCompatible(areaId, nodeId))
+ {
+ PrimitiveType areaType = areas->getPrimitiveType(areaId);
+ switch (areaType)
+ {
+ case PrimitiveType::Plane:
+ {
+ isClose = distanceToPlane(
+ nodes->getCoordinates(nodeId),
+ areas->getAffinePoint(areaId),
+ areas->getNormal(areaId)) < DELTA_PLANE;
+ }
+ break;
+ case PrimitiveType::Sphere:
+ {
+ double distanceToCenter = distanceToSphere(
+ nodes->getCoordinates(nodeId),
+ areas->getCenter(areaId));
+ double radius = areas->getRadius(areaId);
+ isClose = fabs((distanceToCenter - radius) / radius) < DELTA_SPHERE;
+ }
+ break;
+ case PrimitiveType::Cylinder:
+ {
+ double distance = distanceToCylinder(
+ nodes->getCoordinates(nodeId),
+ areas->getAxis(areaId),
+ areas->getAxisPoint(areaId));
+ double radius = areas->getRadius(areaId);
+ isClose = fabs((distance - radius) / radius) < DELTA_CYLINDER;
+ }
+ break;
+ case PrimitiveType::Cone:
+ {
+ double radius = areas->getRadius(areaId);
+ isClose = distanceToCone(
+ nodes->getCoordinates(nodeId),
+ areas->getAxis(areaId),
+ areas->getApex(areaId),
+ areas->getAngle(areaId)) /
+ fabs(radius) <
+ DELTA_CONE;
+ }
+ break;
+ case PrimitiveType::Torus:
+ case PrimitiveType::Unknown:
+ default:
+ break;
+ }
+ }
+ return isClose;
+}
+
+bool AreasBuilder::isInvalidCylinderNode(mcIdType nodeId) const
+{
+ mcIdType areaId = areas->getAreaId(nodeId);
+ if (areaId != -1 &&
+ nodes->getPrimitiveType(nodeId) == PrimitiveType::Cylinder &&
+ areas->getPrimitiveType(areaId) == PrimitiveType::Cylinder)
+ {
+ double angle = MathOps::computeAngle(
+ nodes->getWeakDirection(nodeId),
+ areas->getAxis(areaId));
+ return angle >= THETA_MAX_CYLINDER && angle <= (M_PI - THETA_MAX_CYLINDER);
+ }
+ return false;
+}
+
+double AreasBuilder::distanceToPlane(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &point,
+ const std::array<double, 3> &normal)
+{
+ std::array<double, 3> vec{0.0, 0.0, 0.0};
+ for (size_t i = 0; i < nodeCoords.size(); ++i)
+ vec[i] = nodeCoords[i] - point[i];
+ return fabs(MathOps::dot(vec, normal)) / MathOps::computeNorm(normal);
+}
+
+double AreasBuilder::distanceToSphere(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> ¢er)
+{
+ return MathOps::computeNorm(
+ std::array<double, 3>{
+ nodeCoords[0] - center[0],
+ nodeCoords[1] - center[1],
+ nodeCoords[2] - center[2]});
+}
+
+double AreasBuilder::distanceToCylinder(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &axis,
+ const std::array<double, 3> &axisPoint)
+{
+
+ std::array<double, 3> pa = {
+ axisPoint[0] - nodeCoords[0],
+ axisPoint[1] - nodeCoords[1],
+ axisPoint[2] - nodeCoords[2]};
+ double innerProduct = MathOps::dot(pa, axis);
+ return MathOps::computeNorm(
+ std::array<double, 3>({pa[0] - innerProduct * axis[0],
+ pa[1] - innerProduct * axis[1],
+ pa[2] - innerProduct * axis[2]}));
+}
+
+double AreasBuilder::distanceToCone(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &axis,
+ const std::array<double, 3> &apex,
+ double angle)
+{
+ std::array<double, 3> ps{
+ apex[0] - nodeCoords[0],
+ apex[1] - nodeCoords[1],
+ apex[2] - nodeCoords[2]};
+ std::array<double, 3> v(axis);
+ if (MathOps::dot(axis, ps) <= 0)
+ {
+ for (size_t i = 0; i < 3; ++i)
+ v[i] *= -1;
+ }
+ double a = MathOps::computeNorm(
+ MathOps::cross(ps, v));
+ double b = MathOps::dot(ps, v);
+ return fabs(a * cos(angle) - b * sin(angle));
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __AREASBUILDER_HXX__
+#define __AREASBUILDER_HXX__
+
+#include "Nodes.hxx"
+#include "Areas.hxx"
+
+namespace MEDCoupling
+{
+ class AreasBuilder
+ {
+ public:
+ AreasBuilder(const Nodes *nodes);
+
+ void build();
+
+ Areas *getAreas() const;
+
+ private:
+ void explore();
+ void expand();
+ void rebuild();
+ void exploreAreas();
+ void expandAreas();
+ void expandAreasByType(PrimitiveType primitive);
+ void rebuildInvalidAreas();
+ void filterHighPass();
+ bool doesItMatch(mcIdType areaId, mcIdType nodeId) const;
+ bool doesItBelong(mcIdType areaId, mcIdType nodeId) const;
+ bool isInvalidCylinderNode(mcIdType nodeId) const;
+
+ static double distanceToPlane(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &point,
+ const std::array<double, 3> &normal);
+ static double distanceToSphere(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> ¢er);
+ static double distanceToCylinder(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &axis,
+ const std::array<double, 3> &axisPoint);
+ static double distanceToCone(
+ const std::array<double, 3> &nodeCoords,
+ const std::array<double, 3> &axis,
+ const std::array<double, 3> &apex,
+ double angle);
+
+ const Nodes *nodes;
+ Areas *areas;
+
+ size_t threshold = 5;
+ };
+};
+
+#endif // __AREASBUILDER_HXX__
\ No newline at end of file
--- /dev/null
+# Copyright (C) 2012-2024 CEA, EDF
+#
+# 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 2.1 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, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+
+ADD_DEFINITIONS(${HDF5_DEFINITIONS} ${MEDFILE_DEFINITIONS})
+
+IF (NOT DEFINED MSVC)
+ ADD_DEFINITIONS(-Wsign-compare -Wconversion)
+ENDIF()
+
+IF(MEDCOUPLING_ENABLE_PYTHON)
+ ADD_SUBDIRECTORY(Swig)
+ENDIF(MEDCOUPLING_ENABLE_PYTHON)
+
+IF(MEDCOUPLING_BUILD_TESTS)
+ ADD_SUBDIRECTORY(Test)
+ENDIF(MEDCOUPLING_BUILD_TESTS)
+
+INCLUDE_DIRECTORIES(
+ ${MEDFILE_INCLUDE_DIRS}
+ ${HDF5_INCLUDE_DIRS}
+ ${LAPACKE_INCLUDE_DIRS}
+ ${CMAKE_CURRENT_SOURCE_DIR}/../MEDLoader
+ ${CMAKE_CURRENT_SOURCE_DIR}/../MEDCoupling
+ ${CMAKE_CURRENT_SOURCE_DIR}/../INTERP_KERNEL
+ ${CMAKE_CURRENT_SOURCE_DIR}/../INTERP_KERNEL/Bases
+)
+
+SET(shaperecogn_SOURCES
+ MathOps.cxx
+ Nodes.cxx
+ NodesBuilder.cxx
+ Areas.cxx
+ AreasBuilder.cxx
+ ShapeRecognMeshBuilder.cxx
+ ShapeRecognMesh.cxx
+)
+
+ADD_LIBRARY(shaperecogn ${shaperecogn_SOURCES})
+SET_TARGET_PROPERTIES(shaperecogn PROPERTIES COMPILE_FLAGS "")
+TARGET_LINK_LIBRARIES(shaperecogn medcouplingcpp medloader ${MEDFILE_C_LIBRARIES} ${HDF5_LIBRARIES} ${LAPACK_LIBRARIES})
+INSTALL(TARGETS shaperecogn EXPORT ${PROJECT_NAME}TargetGroup DESTINATION ${MEDCOUPLING_INSTALL_LIBS})
+
+FILE(GLOB shaperecogn_HEADERS_HXX "${CMAKE_CURRENT_SOURCE_DIR}/*.hxx")
+INSTALL(FILES ${shaperecogn_HEADERS_HXX} DESTINATION ${MEDCOUPLING_INSTALL_HEADERS})
+
+# To allow usage as SWIG dependencies:
+SET(shaperecogn_HEADERS_HXX PARENT_SCOPE)
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#include "MathOps.hxx"
+
+#include <algorithm>
+#include <lapacke.h>
+#include <cblas.h>
+#include <iostream>
+#include <cfloat>
+
+using namespace MEDCoupling;
+
+std::vector<double> MathOps::lstsq(
+ std::vector<double> &a, const std::vector<double> &b)
+{
+ int m = (int)b.size();
+ int n = 3;
+ int nrhs = 1;
+ return lstsq(a, b, m, n, nrhs);
+}
+
+std::vector<double> MathOps::lstsq(
+ std::vector<double> &a,
+ const std::vector<double> &b,
+ int m, int n, int nrhs)
+{
+ int ldb = std::max<int>(m, n);
+ int lds = std::min<int>(m, n);
+ std::vector<double> x(ldb, 0.0);
+ for (size_t i = 0; i < b.size(); ++i)
+ x[i] = b[i];
+ std::vector<double> s(lds, 0.0);
+ double rcond = DBL_EPSILON * (double)ldb; // same value as numpy.linalg.lstsq
+ int rank = 0;
+ int info = LAPACKE_dgelsd(
+ LAPACK_COL_MAJOR,
+ m, n, nrhs,
+ a.data(), m,
+ x.data(), ldb,
+ s.data(),
+ rcond,
+ &rank);
+ return x;
+}
+
+std::vector<double> MathOps::lstsqRow(
+ std::vector<double> &a,
+ const std::vector<double> &b)
+{
+ int m = b.size();
+ int n = 3;
+ int nrhs = 1;
+ int ldb = std::max<int>(m, n);
+ int lds = std::min<int>(m, n);
+ std::vector<double> x(ldb, 0.0);
+ for (size_t i = 0; i < b.size(); ++i)
+ x[i] = b[i];
+ std::vector<double> s(lds, 0.0);
+ double rcond = DBL_EPSILON * (double)ldb; // same value as numpy.linalg.lstsq
+ int rank = 0;
+ int info = LAPACKE_dgelsd(
+ LAPACK_ROW_MAJOR,
+ m, n, nrhs,
+ a.data(), n,
+ x.data(), nrhs,
+ s.data(),
+ rcond,
+ &rank);
+ return x;
+}
+std::array<double, 3> MathOps::cross(const std::array<double, 3> &a, const std::array<double, 3> &b)
+{
+ std::array<double, 3> c{0.0, 0.0, 0.0};
+ c[0] = a[1] * b[2] - a[2] * b[1];
+ c[1] = a[2] * b[0] - a[0] * b[2];
+ c[2] = a[0] * b[1] - a[1] * b[0];
+ return c;
+}
+
+std::array<double, 3> MathOps::normalize(const std::array<double, 3> &a)
+{
+ std::array<double, 3> an{0.0, 0.0, 0.0};
+ double n = computeNorm(a);
+ for (size_t i = 0; i < 3; i++)
+ an[i] = a[i] / n;
+ return an;
+}
+
+double MathOps::computeNorm(const std::array<double, 3> &a)
+{
+ double n = 0;
+ for (size_t i = 0; i < 3; i++)
+ n += a[i] * a[i];
+ return sqrt(n);
+}
+
+std::vector<double> MathOps::computeNorm(
+ const std::vector<double> &a)
+{
+ size_t s = a.size() / 3;
+ std::vector<double> n(s, 0.0);
+ for (size_t i = 0; i < s; i++)
+ {
+ for (size_t j = 0; j < 3; j++)
+ n[i] += a[3 * i + j] * a[3 * i + j];
+ n[i] = sqrt(n[i]);
+ }
+ return n;
+}
+
+double MathOps::dot(const std::array<double, 3> &a, const std::array<double, 3> &b)
+{
+ double d = 0.0;
+ for (size_t i = 0; i < 3; i++)
+ d += a[i] * b[i];
+ return d;
+}
+
+std::vector<double> MathOps::dot(const std::vector<double> &a, const std::array<double, 3> &b)
+{
+ size_t nbNodes = a.size() / 3;
+ std::vector<double> d(nbNodes, 0.0);
+ cblas_dgemv(
+ CBLAS_LAYOUT::CblasRowMajor,
+ CBLAS_TRANSPOSE::CblasNoTrans,
+ (int)nbNodes, 3, 1.0,
+ a.data(), 3,
+ b.data(), 1,
+ 0.0, d.data(), 1);
+ return d;
+}
+
+double MathOps::mean(const std::vector<double> &values)
+{
+ double mean = 0.0;
+ for (double value : values)
+ mean += value;
+ return mean / values.size();
+}
+
+std::array<double, 3> MathOps::meanCoordinates(const std::vector<double> &coordinates)
+{
+ std::array<double, 3> coordsMean{0.0, 0.0, 0.0};
+ size_t nbNodes = coordinates.size() / 3;
+ for (size_t nodeId = 0; nodeId < nbNodes; ++nodeId)
+ {
+ coordsMean[0] += coordinates[3 * nodeId];
+ coordsMean[1] += coordinates[3 * nodeId + 1];
+ coordsMean[2] += coordinates[3 * nodeId + 2];
+ }
+ coordsMean[0] /= nbNodes;
+ coordsMean[1] /= nbNodes;
+ coordsMean[2] /= nbNodes;
+ return coordsMean;
+}
+
+std::array<double, 9> MathOps::computeCov(
+ const std::vector<double> &coordinates)
+{
+ std::array<double, 9> covMatrix;
+ covMatrix.fill(0);
+ size_t nbNodes = coordinates.size() / 3;
+ // Center the coordinates
+ std::vector<double> coordsCentered(coordinates);
+ std::array<double, 3> coordsMean = meanCoordinates(coordinates);
+ for (size_t nodeId = 0; nodeId < nbNodes; ++nodeId)
+ {
+ coordsCentered[3 * nodeId] -= coordsMean[0];
+ coordsCentered[3 * nodeId + 1] -= coordsMean[1];
+ coordsCentered[3 * nodeId + 2] -= coordsMean[2];
+ }
+ cblas_dgemm(
+ CBLAS_LAYOUT::CblasColMajor,
+ CBLAS_TRANSPOSE::CblasNoTrans,
+ CBLAS_TRANSPOSE::CblasTrans,
+ 3, 3, (int)nbNodes, 1.0,
+ coordsCentered.data(), 3,
+ coordsCentered.data(), 3,
+ 0.0, covMatrix.data(), 3);
+ for (size_t i = 0; i < 9; ++i)
+ {
+ covMatrix[i] /= (double)nbNodes - 1;
+ }
+ return covMatrix;
+}
+
+std::array<double, 9> MathOps::computePCA(const std::vector<double> &coordinates)
+{
+ std::array<double, 9> covMatrix = computeCov(coordinates);
+ std::array<double, 3> eigenValues{0.0, 0.0, 0.0};
+ LAPACKE_dsyevd(
+ LAPACK_COL_MAJOR,
+ 'V',
+ 'U',
+ 3,
+ covMatrix.data(),
+ 3,
+ eigenValues.data());
+ return covMatrix;
+}
+
+std::array<double, 3> MathOps::computePCAFirstAxis(const std::vector<double> &coordinates)
+{
+ std::array<double, 9> pca = computePCA(coordinates);
+ // The eignvalues are in ascending order so the first axis correspond to the last vector
+ return {pca[6], pca[7], pca[8]};
+}
+
+std::array<double, 3> MathOps::computePCAThirdAxis(const std::vector<double> &coordinates)
+{
+ std::array<double, 9> pca = computePCA(coordinates);
+ // The eignvalues are in ascending order so the third axis correspond to the first vector
+ return {pca[0], pca[1], pca[2]};
+}
+
+double MathOps::computeQuantile(const std::vector<double> &values, double q)
+{
+ std::vector<double> sortedValues(values);
+ std::sort(sortedValues.begin(), sortedValues.end());
+ double pos = q * (sortedValues.size() - 1);
+ size_t index = static_cast<size_t>(pos);
+ if (pos == index)
+ return sortedValues[index];
+ else
+ {
+ double frac = pos - index;
+ return sortedValues[index] * (1 - frac) + sortedValues[index + 1] * frac;
+ }
+}
+
+double MathOps::computeAngle(const std::array<double, 3> &direction, std::array<double, 3> axis)
+{
+ double angle = dot(direction, axis);
+ double normAxis = computeNorm(axis);
+ double normDirection = computeNorm(direction);
+ angle /= normDirection * normAxis;
+ if (fabs(angle) >= 1)
+ return 0;
+ else
+ return acos(angle);
+}
+
+std::vector<double> MathOps::computeAngles(
+ const std::vector<double> &directions, std::array<double, 3> axis)
+{
+ size_t nbDirections = directions.size() / 3;
+ std::vector<double> angles(nbDirections, 0.0);
+ cblas_dgemv(
+ CBLAS_LAYOUT::CblasRowMajor,
+ CBLAS_TRANSPOSE::CblasNoTrans,
+ nbDirections, 3, 1.0,
+ directions.data(), 3,
+ axis.data(), 1,
+ 0.0, angles.data(), 1);
+ double normAxis = computeNorm(axis);
+ std::vector<double> normDirections = computeNorm(directions);
+ for (size_t i = 0; i < nbDirections; ++i)
+ {
+ angles[i] /= normAxis * normDirections[i];
+ if (fabs(angles[i]) >= 1.0)
+ angles[i] = 0.0;
+ angles[i] = acos(angles[i]);
+ }
+ return angles;
+}
+
+double MathOps::computeOrientedAngle(const std::array<double, 3> &normal, const std::array<double, 3> &vector1, const std::array<double, 3> &vector2)
+{
+ double angle = computeAngle(vector1, vector2);
+ if (dot(cross(vector1, vector2), normal) >= 0.0)
+ return angle;
+ else
+ return -angle;
+}
+
+double MathOps::computeVariance(std::vector<double> values)
+{
+ size_t n = values.size();
+ double m = mean(values);
+ double d2 = 0;
+ for (size_t i = 0; i < n; ++i)
+ d2 += pow(values[i] - m, 2);
+ return d2 / (double)n;
+}
+
+std::array<double, 6> MathOps::computeBaseFromNormal(std::array<double, 3> normal)
+{
+ std::array<double, 3> n_normal = normalize(normal);
+ std::array<double, 3> s;
+ std::array<double, 1> u;
+ std::array<double, 9> v;
+ LAPACKE_dgesdd(
+ LAPACK_COL_MAJOR,
+ 'A', 1, 3, n_normal.data(), 1,
+ s.data(),
+ u.data(), 1,
+ v.data(), 3);
+ std::array<double, 3> v1 = {v[3], v[4], v[5]};
+ std::array<double, 3> v2 = {v[6], v[7], v[8]};
+ std::array<double, 3> u1 = normalize(v1);
+ std::array<double, 3> u2;
+ u2.fill(0.0);
+ double innerProd = dot(u1, v2);
+ for (size_t i = 0; i < 3; ++i)
+ u2[i] = v2[i] - innerProd * u1[i];
+ u2 = normalize(u2);
+ double sign = dot(cross(u1, u2), normal);
+ if (sign < 0)
+ {
+ for (size_t i = 0; i < 3; ++i)
+ u1[i] *= -1;
+ }
+ return {u1[0], u1[1], u1[2], u2[0], u2[1], u2[2]};
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __MATHOPS_HXX__
+#define __MATHOPS_HXX__
+
+#include <vector>
+#include <array>
+
+namespace MEDCoupling
+{
+ class MathOps
+ {
+ public:
+ static std::vector<double> lstsq(std::vector<double> &a, const std::vector<double> &b);
+ static std::vector<double> lstsq(
+ std::vector<double> &a,
+ const std::vector<double> &b,
+ int m,
+ int n = 3,
+ int nrhs = 1);
+ static std::vector<double> lstsqRow(
+ std::vector<double> &a,
+ const std::vector<double> &b);
+ static std::array<double, 3> cross(const std::array<double, 3> &a, const std::array<double, 3> &b);
+ static std::array<double, 3> normalize(const std::array<double, 3> &a);
+ static double computeNorm(const std::array<double, 3> &a);
+ static std::vector<double> computeNorm(const std::vector<double> &a);
+ static double dot(const std::array<double, 3> &a, const std::array<double, 3> &b);
+ static std::vector<double> dot(const std::vector<double> &a, const std::array<double, 3> &b);
+ static double mean(const std::vector<double> &values);
+ static std::array<double, 3> meanCoordinates(const std::vector<double> &coordinates);
+ static std::array<double, 9> computeCov(const std::vector<double> &coordinates);
+ static std::array<double, 9> computePCA(const std::vector<double> &coordinates);
+ static std::array<double, 3> computePCAFirstAxis(const std::vector<double> &coordinates);
+ static std::array<double, 3> computePCAThirdAxis(const std::vector<double> &coordinates);
+ static double computeQuantile(
+ const std::vector<double> &values,
+ double q);
+ static double computeAngle(
+ const std::array<double, 3> &direction,
+ std::array<double, 3> axis);
+ static std::vector<double> computeAngles(
+ const std::vector<double> &directions,
+ std::array<double, 3> axis);
+ static double computeOrientedAngle(
+ const std::array<double, 3> &normal,
+ const std::array<double, 3> &vector1,
+ const std::array<double, 3> &vector2);
+ static double computeVariance(std::vector<double> values);
+ static std::array<double, 6> computeBaseFromNormal(std::array<double, 3> normal);
+ };
+};
+
+#endif //__MATHOPS_HXX__
--- /dev/null
+#include "Nodes.hxx"
+
+using namespace MEDCoupling;
+
+Nodes::Nodes(
+ const MEDCouplingUMesh *mesh,
+ const DataArrayInt64 *neighbors,
+ const DataArrayInt64 *neighborsIdx)
+ : mesh(mesh), coords(mesh->getCoords()), neighbors(neighbors), neighborsIdx(neighborsIdx)
+{
+}
+
+Nodes::~Nodes()
+{
+ mesh->decrRef();
+ neighbors->decrRef();
+ neighborsIdx->decrRef();
+}
+
+mcIdType Nodes::getNbNodes() const
+{
+ return coords->getNumberOfTuples();
+}
+
+const std::vector<double> &Nodes::getK1() const
+{
+ return k1;
+}
+
+const std::vector<double> &Nodes::getK2() const
+{
+ return k2;
+}
+
+const std::vector<PrimitiveType> &Nodes::getPrimitiveType() const
+{
+ return primitives;
+}
+
+const std::vector<double> &Nodes::getNormals() const
+{
+ return normals;
+}
+
+const std::vector<double> &Nodes::getWeakDirections() const
+{
+ return weakDirections;
+}
+
+const std::vector<double> &Nodes::getMainDirections() const
+{
+ return mainDirections;
+}
+
+std::array<double, 3> Nodes::getNormal(mcIdType nodeId) const
+{
+ return {normals[3 * nodeId], normals[3 * nodeId + 1], normals[3 * nodeId + 2]};
+}
+
+double Nodes::getK1(mcIdType nodeId) const
+{
+ return k1[nodeId];
+}
+
+double Nodes::getK2(mcIdType nodeId) const
+{
+ return k2[nodeId];
+}
+
+double Nodes::getKdiff0(mcIdType nodeId) const
+{
+ return fabs(k1[nodeId]) > fabs(k2[nodeId]) ? k1[nodeId] : k2[nodeId];
+}
+
+double Nodes::getAdimK1(mcIdType nodeId) const
+{
+ return adimK1[nodeId];
+}
+
+double Nodes::getAdimK2(mcIdType nodeId) const
+{
+ return adimK2[nodeId];
+}
+
+double Nodes::getAdimKdiff0(mcIdType nodeId) const
+{
+ return fabs(adimK1[nodeId]) > fabs(adimK2[nodeId]) ? adimK1[nodeId] : adimK2[nodeId];
+}
+
+const std::vector<mcIdType> Nodes::getNeighbors(mcIdType nodeId) const
+{
+ mcIdType start = neighborsIdx->getIJ(nodeId, 0);
+ mcIdType nbNeighbors = neighborsIdx->getIJ(nodeId + 1, 0) - start;
+
+ std::vector<mcIdType> neighborNodes(nbNeighbors, 0);
+ for (mcIdType i = 0; i < nbNeighbors; i++)
+ neighborNodes[i] = neighbors->getIJ(start + i, 0);
+ return neighborNodes;
+}
+
+PrimitiveType Nodes::getPrimitiveType(mcIdType nodeId) const
+{
+ return primitives[nodeId];
+}
+
+std::array<double, 3> Nodes::getWeakDirection(mcIdType nodeId) const
+{
+ return {weakDirections[3 * nodeId],
+ weakDirections[3 * nodeId + 1],
+ weakDirections[3 * nodeId + 2]};
+}
+
+std::array<double, 3> Nodes::getMainDirection(mcIdType nodeId) const
+{
+ return {mainDirections[3 * nodeId],
+ mainDirections[3 * nodeId + 1],
+ mainDirections[3 * nodeId + 2]};
+}
+
+std::array<double, 3> Nodes::getCoordinates(mcIdType nodeId) const
+{
+ std::array<double, 3> nodeCoords;
+ coords->getTuple(nodeId, nodeCoords.data());
+ return nodeCoords;
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __NODES_HXX__
+#define __NODES_HXX__
+
+#include <vector>
+#include "MEDCouplingUMesh.hxx"
+#include "PrimitiveType.hxx"
+
+namespace MEDCoupling
+{
+ class Nodes
+ {
+ public:
+ friend class NodesBuilder;
+ Nodes(const MEDCouplingUMesh *mesh,
+ const DataArrayInt64 *neighbors,
+ const DataArrayInt64 *neighborsIdx);
+ ~Nodes();
+
+ mcIdType getNbNodes() const;
+ const std::vector<double> &getK1() const;
+ const std::vector<double> &getK2() const;
+ const std::vector<PrimitiveType> &getPrimitiveType() const;
+ const std::vector<double> &getNormals() const;
+ const std::vector<double> &getWeakDirections() const;
+ const std::vector<double> &getMainDirections() const;
+
+ std::array<double, 3> getNormal(mcIdType nodeId) const;
+ double getK1(mcIdType nodeId) const;
+ double getK2(mcIdType nodeId) const;
+ double getKdiff0(mcIdType nodeId) const;
+ double getAdimK1(mcIdType nodeId) const;
+ double getAdimK2(mcIdType nodeId) const;
+ double getAdimKdiff0(mcIdType nodeId) const;
+ const std::vector<mcIdType> getNeighbors(mcIdType nodeId) const;
+ PrimitiveType getPrimitiveType(mcIdType nodeId) const;
+ std::array<double, 3> getWeakDirection(mcIdType nodeId) const;
+ std::array<double, 3> getMainDirection(mcIdType nodeId) const;
+ std::array<double, 3> getCoordinates(mcIdType nodeId) const;
+
+ private:
+ const MEDCouplingUMesh *mesh;
+ const DataArrayDouble *coords;
+ // normals 3 * nbNodes
+ std::vector<double> normals;
+ // neighbors
+ const DataArrayInt64 *neighbors;
+ const DataArrayInt64 *neighborsIdx;
+ // curvature
+ std::vector<double> k1;
+ std::vector<double> k2;
+ std::vector<double> adimK1;
+ std::vector<double> adimK2;
+ std::vector<double> weakDirections;
+ std::vector<double> mainDirections;
+ std::vector<PrimitiveType> primitives;
+ };
+};
+
+#endif //__NODES_HXX__
--- /dev/null
+
+#include "NodesBuilder.hxx"
+#include "Nodes.hxx"
+#include "MathOps.hxx"
+#include "ShapeRecongConstants.hxx"
+
+using namespace MEDCoupling;
+
+NodesBuilder::NodesBuilder(
+ const MEDCouplingUMesh *mesh) : mesh(mesh)
+{
+}
+
+Nodes *NodesBuilder::build()
+{
+ DataArrayInt64 *neighbors;
+ DataArrayInt64 *neighborsIdx;
+ mesh->computeNeighborsOfNodes(neighbors, neighborsIdx);
+ nodes = new Nodes(mesh, neighbors, neighborsIdx);
+ computeNormals();
+ computeCurvatures();
+ return nodes;
+}
+
+void NodesBuilder::computeNormals()
+{
+ mcIdType nbNodes = mesh->getNumberOfNodes();
+ mcIdType nbCells = mesh->getNumberOfCells();
+ std::array<double, 3> cellNormal;
+ cellNormal.fill(0.0);
+ // product of normal AreaByCell
+ std::vector<double> prodNormalAreaByCell(3 * nbCells, 0.0);
+ for (int cellId = 0; cellId < nbCells; cellId++)
+ {
+ std::vector<mcIdType> nodeIds;
+ mesh->getNodeIdsOfCell(cellId, nodeIds);
+ computeCellNormal(nodeIds, cellNormal);
+ prodNormalAreaByCell[3 * cellId + 0] = cellNormal[0] / 2.0;
+ prodNormalAreaByCell[3 * cellId + 1] = cellNormal[1] / 2.0;
+ prodNormalAreaByCell[3 * cellId + 2] = cellNormal[2] / 2.0;
+ }
+ //
+ nodes->normals.resize(3 * nbNodes, 1.0);
+ DataArrayInt64 *revNodal = DataArrayInt64::New();
+ DataArrayInt64 *revNodalIdx = DataArrayInt64::New();
+ mesh->getReverseNodalConnectivity(revNodal, revNodalIdx);
+ for (size_t nodeId = 0; nodeId < (size_t)nbNodes; nodeId++)
+ {
+ int nbCells = revNodalIdx->getIJ(nodeId + 1, 0) -
+ revNodalIdx->getIJ(nodeId, 0);
+ std::vector<mcIdType> cellIds(nbCells, 0);
+ int start = revNodalIdx->getIJ(nodeId, 0);
+ for (size_t i = 0; i < cellIds.size(); ++i)
+ cellIds[i] = revNodal->getIJ(start + i, 0);
+ double normal = 0.0;
+ for (size_t i = 0; i < 3; i++)
+ {
+ nodes->normals[3 * nodeId + i] = 0.0;
+ for (mcIdType j = 0; j < nbCells; j++)
+ {
+ nodes->normals[3 * nodeId + i] +=
+ prodNormalAreaByCell[3 * cellIds[j] + i];
+ }
+ nodes->normals[3 * nodeId + i] /= (double)nbCells;
+ normal += pow(nodes->normals[3 * nodeId + i], 2);
+ }
+ for (size_t i = 0; i < 3; i++)
+ nodes->normals[3 * nodeId + i] /= sqrt(normal);
+ }
+ revNodal->decrRef();
+ revNodalIdx->decrRef();
+}
+
+void NodesBuilder::computeCurvatures(double tol)
+{
+ mcIdType nbNodes = mesh->getNumberOfNodes();
+ nodes->k1.resize(nbNodes);
+ nodes->k2.resize(nbNodes);
+ nodes->adimK1.resize(nbNodes);
+ nodes->adimK2.resize(nbNodes);
+ nodes->primitives.resize(nbNodes);
+ nodes->weakDirections.resize(3 * nbNodes);
+ nodes->mainDirections.resize(3 * nbNodes);
+ for (mcIdType nodeId = 0; nodeId < nbNodes; nodeId++)
+ computeCurvatures(nodeId, tol);
+}
+
+void NodesBuilder::computeCurvatures(mcIdType nodeId, double tol)
+{
+ std::array<double, 3> normal = nodes->getNormal(nodeId);
+ std::vector<mcIdType> neighborIds = nodes->getNeighbors(nodeId);
+ double theta0 = 0.0;
+ double k1 = 0.0;
+ double k2 = 0.0;
+ double adimK1 = 0.0;
+ double adimK2 = 0.0;
+ PrimitiveType primitive = PrimitiveType::Unknown;
+ std::array<double, 3> mainDirection{0.0, 0.0, 0.0};
+ std::array<double, 3> weakDirection{0.0, 0.0, 0.0};
+ if (neighborIds.size() > 0)
+ {
+ std::vector<double> discreteCurvatures = computeDiscreteCurvatures(nodeId, neighborIds);
+ std::vector<double> tangents = computeTangentDirections(nodeId, neighborIds);
+ mcIdType maxCurvatureId = std::distance(
+ discreteCurvatures.begin(),
+ std::max_element(discreteCurvatures.begin(), discreteCurvatures.end()));
+ std::array<double, 3> e1{
+ tangents[3 * maxCurvatureId],
+ tangents[3 * maxCurvatureId + 1],
+ tangents[3 * maxCurvatureId + 2]};
+ std::array<double, 3> e2 = MathOps::normalize(MathOps::cross(e1, normal));
+ std::vector<double> coeffs = computeNormalCurvatureCoefficients(
+ discreteCurvatures, tangents, normal, e1);
+ double a = coeffs[0], b = coeffs[1], c = coeffs[2];
+ double h = (a + c) / 2.0;
+ double kg = a * c - pow(b, 2) / 4.0;
+ k1 = h + sqrt(fabs(pow(h, 2) - kg));
+ k2 = h - sqrt(fabs(pow(h, 2) - kg));
+ if (fabs(k1 - k2) >= tol)
+ theta0 = 0.5 * asin(b / (k1 - k2));
+ std::array<double, 3> direction1{0.0, 0.0, 0.0};
+ std::array<double, 3> direction2{0.0, 0.0, 0.0};
+ for (size_t i = 0; i < 3; ++i)
+ {
+ direction1[i] = cos(theta0) * e1[i] + sin(theta0) * e2[i];
+ direction2[i] = -sin(theta0) * e1[i] + cos(theta0) * e2[i];
+ }
+ double averageDistance = computeAverageDistance(nodeId, neighborIds);
+ double adimK1 = k1 * averageDistance;
+ double adimK2 = k2 * averageDistance;
+ double adimKdiff0, adimKis0;
+ if (fabs(k1) > fabs(k2))
+ {
+ adimKdiff0 = adimK1;
+ adimKis0 = adimK2;
+ mainDirection = direction1;
+ weakDirection = direction2;
+ }
+ else
+ {
+ adimKdiff0 = adimK2;
+ adimKis0 = adimK1;
+ mainDirection = direction2;
+ weakDirection = direction1;
+ }
+ primitive = findPrimitiveType(adimK1, adimK2, adimKdiff0, adimKis0);
+ }
+ // Populate nodes
+ nodes->k1[nodeId] = k1;
+ nodes->k2[nodeId] = k2;
+ nodes->adimK1[nodeId] = adimK1;
+ nodes->adimK2[nodeId] = adimK2;
+ for (size_t i = 0; i < 3; ++i)
+ {
+ nodes->weakDirections[3 * nodeId + i] = weakDirection[i];
+ nodes->mainDirections[3 * nodeId + i] = mainDirection[i];
+ }
+ nodes->primitives[nodeId] = primitive;
+}
+
+PrimitiveType NodesBuilder::findPrimitiveType(double k1, double k2, double kdiff0, double kis0) const
+{
+ if ((fabs(k1 - k2) < EPSILON_PRIMITIVE) && (fabs((k1 + k2) / 2) < EPSILON_PRIMITIVE))
+ return PrimitiveType::Plane;
+ else if ((fabs(k1 - k2) < EPSILON_PRIMITIVE) && (fabs((k1 + k2) / 2) > EPSILON_PRIMITIVE))
+ return PrimitiveType::Sphere;
+ else if ((fabs(kdiff0) > EPSILON_PRIMITIVE) && (fabs(kis0) < EPSILON_PRIMITIVE))
+ return PrimitiveType::Cylinder;
+ else if ((fabs(k1) > EPSILON_PRIMITIVE) && (fabs(k2) > EPSILON_PRIMITIVE))
+ return PrimitiveType::Torus;
+ else
+ return PrimitiveType::Unknown;
+}
+
+PrimitiveType NodesBuilder::findPrimitiveType2(double k1, double k2, double kdiff0, double kis0) const
+{
+ double epsilon2 = pow(EPSILON_PRIMITIVE, 2);
+ double diffCurvature = fabs(k1 - k2);
+ double gaussianCurvature = k1 * k2;
+ double meanCurvature = (k1 + k2) / 2.0;
+ if (fabs(k1) < EPSILON_PRIMITIVE &&
+ fabs(k2) < EPSILON_PRIMITIVE &&
+ gaussianCurvature < epsilon2 &&
+ meanCurvature < EPSILON_PRIMITIVE)
+ return PrimitiveType::Plane;
+ else if (diffCurvature < EPSILON_PRIMITIVE && k1 > EPSILON_PRIMITIVE && k2 > EPSILON_PRIMITIVE)
+ return PrimitiveType::Sphere;
+ else if (
+ (fabs(k1) > EPSILON_PRIMITIVE && fabs(k2) < EPSILON_PRIMITIVE) ||
+ (fabs(k1) < EPSILON_PRIMITIVE && fabs(k2) > EPSILON_PRIMITIVE))
+ return PrimitiveType::Cylinder;
+ else if (
+ std::signbit(k1) != std::signbit(k2) ||
+ (fabs(k1) < EPSILON_PRIMITIVE && fabs(k2) < EPSILON_PRIMITIVE))
+ return PrimitiveType::Torus;
+ else
+ return PrimitiveType::Unknown;
+}
+
+std::vector<double> NodesBuilder::computeNormalCurvatureCoefficients(
+ const std::vector<double> &discreteCurvatures,
+ const std::vector<double> &tangents,
+ const std::array<double, 3> &normal,
+ const std::array<double, 3> &e1) const
+{
+ size_t nbNeighbors = discreteCurvatures.size();
+ std::vector<double> a(3 * nbNeighbors, 0.0);
+ for (size_t i = 0; i < nbNeighbors; ++i)
+ {
+ std::array<double, 3> tangent{tangents[3 * i], tangents[3 * i + 1], tangents[3 * i + 2]};
+ double theta = MathOps::computeOrientedAngle(normal, tangent, e1);
+ double cos_theta = cos(theta);
+ double sin_theta = sin(theta);
+ a[i] = pow(cos_theta, 2);
+ a[nbNeighbors + i] = cos_theta * sin_theta;
+ a[2 * nbNeighbors + i] = pow(sin_theta, 2);
+ }
+ return MathOps::lstsq(a, discreteCurvatures);
+}
+
+void NodesBuilder::computeCellNormal(
+ const std::vector<mcIdType> &nodeIds,
+ std::array<double, 3> &cellNormal) const
+{
+ std::vector<double> point1;
+ std::vector<double> point2;
+ std::vector<double> point3;
+ mesh->getCoordinatesOfNode(nodeIds[0], point1);
+ mesh->getCoordinatesOfNode(nodeIds[1], point2);
+ mesh->getCoordinatesOfNode(nodeIds[2], point3);
+ std::array<double, 3> a;
+ a.fill(3);
+ std::array<double, 3> b;
+ b.fill(3);
+ for (int i = 0; i < 3; i++)
+ {
+ a[i] = point2[i] - point1[i];
+ b[i] = point3[i] - point1[i];
+ }
+ cellNormal[0] = a[1] * b[2] - a[2] * b[1];
+ cellNormal[1] = a[2] * b[0] - a[0] * b[2];
+ cellNormal[2] = a[0] * b[1] - a[1] * b[0];
+}
+
+double NodesBuilder::computeAverageDistance(mcIdType nodeId, const std::vector<mcIdType> &neighborIds) const
+{
+ double distance = 0.0;
+ std::array<double, 3>
+ nodeCoords = nodes->getCoordinates(nodeId);
+ for (size_t i = 0; i < neighborIds.size(); ++i)
+ {
+ std::array<double, 3> neighborCoords = nodes->getCoordinates(neighborIds[i]);
+ double distanceToNeighbor = 0.0;
+ for (size_t j = 0; j < 3; ++j)
+ distanceToNeighbor += pow(neighborCoords[j] - nodeCoords[j], 2);
+ distance += sqrt(distanceToNeighbor);
+ }
+ return distance / (double)neighborIds.size();
+}
+
+std::vector<double> NodesBuilder::computeDiscreteCurvatures(
+ mcIdType nodeId,
+ const std::vector<mcIdType> &neighborIds) const
+{
+ std::vector<double> discreteCurvatures(neighborIds.size(), 0.0);
+ for (size_t i = 0; i < neighborIds.size(); ++i)
+ discreteCurvatures[i] = computeDiscreteCurvature(nodeId, neighborIds[i]);
+ return discreteCurvatures;
+}
+
+double NodesBuilder::computeDiscreteCurvature(
+ mcIdType nodeId,
+ mcIdType neighborId) const
+{
+ double curvature = 0.0;
+ double n = 0.0;
+ for (size_t i = 0; i < 3; i++)
+ {
+ double ni = nodes->coords->getIJ(neighborId, i) - nodes->coords->getIJ(nodeId, i);
+ curvature += ni * (nodes->normals[3 * neighborId + i] - nodes->normals[3 * nodeId + i]);
+ n += ni * ni;
+ }
+ return curvature / n;
+}
+
+std::vector<double> NodesBuilder::computeTangentDirections(
+ mcIdType nodeId,
+ const std::vector<mcIdType> &neighborIds) const
+{
+ size_t nbNeighbors = neighborIds.size();
+ std::vector<double> tangent(3 * nbNeighbors, 0.0);
+ for (size_t i = 0; i < nbNeighbors; ++i)
+ {
+ mcIdType neighborId = neighborIds[i];
+ double s = 0.0;
+ for (size_t j = 0; j < 3; j++)
+ {
+ tangent[3 * i + j] = nodes->coords->getIJ(neighborId, j) - nodes->coords->getIJ(nodeId, j);
+ s += tangent[3 * i + j] * nodes->normals[3 * nodeId + j];
+ }
+ double n = 0.0;
+ for (size_t j = 0; j < 3; j++)
+ {
+ tangent[3 * i + j] -= s * nodes->normals[3 * nodeId + j];
+ n += tangent[3 * i + j] * tangent[3 * i + j];
+ }
+ for (size_t j = 0; j < 3; j++)
+ tangent[3 * i + j] /= sqrt(n);
+ }
+ return tangent;
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __NODECURVATURECALCULATOR_HXX__
+#define __NODECURVATURECALCULATOR_HXX__
+
+#include <string>
+#include <vector>
+#include <array>
+#include "MEDCouplingUMesh.hxx"
+#include "PrimitiveType.hxx"
+
+namespace MEDCoupling
+{
+ class Nodes;
+
+ class NodesBuilder
+ {
+ public:
+ NodesBuilder(const MEDCouplingUMesh *);
+
+ Nodes *build();
+
+ private:
+ void computeNormals();
+ void computeCurvatures(double tol = 0.000001);
+ void computeCurvatures(mcIdType nodeId, double tol);
+ PrimitiveType findPrimitiveType(double k1, double k2, double kdiff0, double kis0) const;
+ PrimitiveType findPrimitiveType2(double k1, double k2, double kdiff0, double kis0) const;
+ std::vector<double> computeNormalCurvatureCoefficients(
+ const std::vector<double> &discreteCurvatures,
+ const std::vector<double> &tangents,
+ const std::array<double, 3> &normal,
+ const std::array<double, 3> &e1) const;
+ void computeCellNormal(const std::vector<mcIdType> &nodeIds, std::array<double, 3> &cellNormal) const;
+ double computeAverageDistance(mcIdType nodeId, const std::vector<mcIdType> &neighborIds) const;
+ std::vector<double> computeDiscreteCurvatures(mcIdType nodeId, const std::vector<mcIdType> &neighborIds) const;
+ double computeDiscreteCurvature(mcIdType nodeId, mcIdType neighborId) const;
+ std::vector<double> computeTangentDirections(mcIdType nodeId, const std::vector<mcIdType> &neighborIds) const;
+
+ const MEDCouplingUMesh *mesh;
+ Nodes *nodes;
+ };
+};
+
+#endif //__NODECURVATURECALCULATOR_HXX__
\ No newline at end of file
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __PRIMITIVETYPE_HXX__
+#define __PRIMITIVETYPE_HXX__
+
+#include <string>
+namespace MEDCoupling
+{
+ enum PrimitiveType
+ {
+ Plane = 0,
+ Sphere = 1,
+ Cylinder = 2,
+ Cone = 3,
+ Torus = 4,
+ Unknown = 5
+ };
+
+ inline std::string convertPrimitiveToString(PrimitiveType type)
+ {
+ std::string typeName = "";
+ switch (type)
+ {
+ case PrimitiveType::Plane:
+ typeName = "Plane";
+ break;
+ case PrimitiveType::Sphere:
+ typeName = "Sphere";
+ break;
+ case PrimitiveType::Cylinder:
+ typeName = "Cylinder";
+ break;
+ case PrimitiveType::Cone:
+ typeName = "Cone";
+ break;
+ case PrimitiveType::Torus:
+ typeName = "Torus";
+ break;
+ case PrimitiveType::Unknown:
+ typeName = "Unknown";
+ break;
+ default:
+ break;
+ }
+ return typeName;
+ };
+
+ inline int convertPrimitiveToInt(PrimitiveType type)
+ {
+ int typeInt = 5;
+ switch (type)
+ {
+ case PrimitiveType::Plane:
+ typeInt = 0;
+ break;
+ case PrimitiveType::Sphere:
+ typeInt = 1;
+ break;
+ case PrimitiveType::Cylinder:
+ typeInt = 2;
+ break;
+ case PrimitiveType::Cone:
+ typeInt = 3;
+ break;
+ case PrimitiveType::Torus:
+ typeInt = 4;
+ break;
+ case PrimitiveType::Unknown:
+ default:
+ break;
+ }
+ return typeInt;
+ };
+};
+
+#endif // __PRIMITIVETYPE_HXX__
\ No newline at end of file
--- /dev/null
+# ShapeRecogn
+A tool leveraging the MEDCoupling library for recognizing canonical shapes in 3D mesh files.
+
+## Table Of Contents
+
+1. [Quickstart](#quickstart)
+2. [Installation](#installation)
+3. [Description](#description)
+
+## Quickstart
+
+### With output file
+
+```python
+>> import ShapeRecogn as sr
+>>
+>> shape_recogn_builder = sr.ShapeRecognMeshBuilder("resources/ShapeRecognCone.med")
+>> shape_recogn = shape_recogn_builder.recognize()
+>> shape_recogn.save("ShapeRecognCone_areas.med")
+```
+
+### Without output file
+
+```python
+>> import ShapeRecogn as sr
+>>
+>> shape_recogn_builder = sr.ShapeRecognMeshBuilder("resources/ShapeRecognCone.med")
+>> shape_recogn = shape_recogn_builder.recognize()
+>> radius_field = shape_recogn.getRadius()
+>> radius_field
+MEDCouplingFieldDouble C++ instance at 0x55f3418c6700. Name : "Radius (Area)".
+Nature of field : NoNature.
+P1 spatial discretization.
+
+Mesh info : MEDCouplingUMesh C++ instance at 0x55f34171a2d0. Name : "Mesh_1". Mesh dimension : 2. Space dimension : 3.
+
+Array info : DataArrayDouble C++ instance at 0x55f3418b40e0. Number of tuples : 957. Number of components : 1.
+[9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, 9.8948383131651862, ... ]
+```
+
+## Installation
+
+### Prerequisites
+
+**LAPACKE Library**: ShapeRecogn requires the [LAPACKE](https://www.netlib.org/lapack/lapacke.html) library for numerical computations. Ensure it is installed on your system.
+
+### CMAKE configuration
+
+Run CMake to configure the build. Make sure to enable ShapeRecogn by setting the `-DMEDCOUPLING_ENABLE_SHAPERECOGN=ON` option.
+
+## Description
+
+ShapeRecogn is a tool that leverages the MEDCoupling library to recognize
+canonical shapes from 3D meshes using triangle cells.
+
+The tool is based on the thesis work of Roseline Bénière *[Reconstruction d’un modèle B-Rep à partir d’un maillage 3D](https://theses.fr/2012MON20034)*, and it recognizes five canonical shapes:
+ - Plane (0)
+ - Sphere (1)
+ - Cylinder (2)
+ - Cone (3)
+ - Torus (4)
+ - Unknown (5) (When the algorithm failed the recognition)
+
+The recognized shapes are divided into areas within the mesh.
+The tool also generates an output file with detailed fields describing the areas and their properties. This makes it easier to analyze or further process the mesh data.
+ - Area Id: To distinguish the areas
+ - Primitive Type (Area): One of the canonical shape with the id describe above
+ - Normal (Area):
+ * Normal of a plane area
+ - Minor Radius (Area)
+ * Minor radius of a torus area
+ - Radius (Area)
+ * Radius of a sphere area
+ * Radius of a cylinder area
+ * Radius of the base of a cone area
+ * Major radius of a torus area
+ - Angle (Area)
+ * Angle between the central axis and the slant of a cone area
+ - Center (Area)
+ * Center of a sphere area
+ * Center of a torus area
+ - Axis (Area)
+ * Central axis of a cylinder area
+ * Central axis of a cone area
+ - Apex (Area)
+ * Apex of a cone area
+
+Some intermediate computation values concerning the nodes are also available:
+ - K1 (Node) and K2 (Node): the *[principal curvatures](https://en.wikipedia.org/wiki/Principal_curvature)*
+ - Primitive Type (Node) : One of the canonical shape with the id describe above. The primitive type is deduced usint K1 and K2 values.
+ - Normal (Node): Normal of the nodes using neighbor nodes
+
+Each field can be retrieved as a `MEDCouplingDoubleField` using the following methods of the `ShapeRecognMesh` class:
+ - `getAreaId()`
+ - `getAreaPrimitiveType()`
+ - `getAreaNormal()`
+ - `getMinorRadius()`
+ - `getRadius()`
+ - `getAngle()`
+ - `getCenter()`
+ - `getAxis()`
+ - `getApex()`
+ - `getNodeK1()`
+ - `getNodeK2()`
+ - `getNodePrimitiveType()`
+ - `getNodeNormal()`
--- /dev/null
+#include "ShapeRecognMesh.hxx"
+
+#include "MEDLoader.hxx"
+
+using namespace MEDCoupling;
+
+ShapeRecognMesh::ShapeRecognMesh()
+ : nodeK1(0), nodeK2(0), nodePrimitiveType(0),
+ nodeNormal(0), areaId(0), areaPrimitiveType(0),
+ areaNormal(0), minorRadius(0), radius(0),
+ angle(0), center(0), axis(0), apex(0)
+{
+}
+ShapeRecognMesh::~ShapeRecognMesh()
+{
+ nodeK1->decrRef();
+ nodeK2->decrRef();
+ nodePrimitiveType->decrRef();
+ nodeNormal->decrRef();
+ areaId->decrRef();
+ areaPrimitiveType->decrRef();
+ areaNormal->decrRef();
+ minorRadius->decrRef();
+ radius->decrRef();
+ angle->decrRef();
+ center->decrRef();
+ axis->decrRef();
+ apex->decrRef();
+}
+
+std::size_t ShapeRecognMesh::getHeapMemorySizeWithoutChildren() const
+{
+ return 0;
+}
+
+std::vector<const BigMemoryObject *> ShapeRecognMesh::getDirectChildrenWithNull() const
+{
+ std::vector<const BigMemoryObject *> ret;
+ ret.push_back(nodeK1);
+ ret.push_back(nodeK2);
+ ret.push_back(nodePrimitiveType);
+ ret.push_back(nodeNormal);
+ ret.push_back(areaId);
+ ret.push_back(areaPrimitiveType);
+ ret.push_back(areaNormal);
+ ret.push_back(minorRadius);
+ ret.push_back(radius);
+ ret.push_back(angle);
+ ret.push_back(center);
+ ret.push_back(axis);
+ ret.push_back(apex);
+ return ret;
+}
+
+ShapeRecognMesh *ShapeRecognMesh::New()
+{
+ return new ShapeRecognMesh;
+}
+
+void ShapeRecognMesh::save(const std::string &outputFile, bool writeFromScratch) const
+{
+ // Nodes
+ // - k1
+ WriteField(outputFile, nodeK1, writeFromScratch);
+ // - k2
+ WriteField(outputFile, nodeK2, false);
+ // - primitive types
+ WriteField(outputFile, nodePrimitiveType, false);
+ // - Normal
+ WriteField(outputFile, nodeNormal, false);
+ // Areas
+ // - Area Id
+ WriteField(outputFile, areaId, false);
+ // - Primitive Types
+ WriteField(outputFile, areaPrimitiveType, false);
+ // - Normal
+ WriteField(outputFile, areaNormal, false);
+ // - Minor Radius
+ WriteField(outputFile, minorRadius, false);
+ // - Radius
+ WriteField(outputFile, radius, false);
+ // - Angle
+ WriteField(outputFile, angle, false);
+ // - Center
+ WriteField(outputFile, center, false);
+ // - Axis
+ WriteField(outputFile, axis, false);
+ // - Apex
+ WriteField(outputFile, apex, false);
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getNodeK1() const
+{
+ return nodeK1;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getNodeK2() const
+{
+ return nodeK2;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getNodePrimitiveType() const
+{
+ return nodePrimitiveType;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getNodeNormal() const
+{
+ return nodeNormal;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getAreaId() const
+{
+ return areaId;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getAreaPrimitiveType() const
+{
+ return areaPrimitiveType;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getAreaNormal() const
+{
+ return areaNormal;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getMinorRadius() const
+{
+ return minorRadius;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getRadius() const
+{
+ return radius;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getAngle() const
+{
+ return angle;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getCenter() const
+{
+ return center;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getAxis() const
+{
+ return axis;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMesh::getApex() const
+{
+ return apex;
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __SHAPERECOGNMESH_HXX__
+#define __SHAPERECOGNMESH_HXX__
+
+#include <string>
+
+#include "MEDCouplingUMesh.hxx"
+#include "MEDCouplingFieldDouble.hxx"
+#include "MEDCouplingRefCountObject.hxx"
+
+namespace MEDCoupling
+{
+ class ShapeRecognMesh : public RefCountObject
+ {
+ friend class ShapeRecognMeshBuilder;
+
+ public:
+ static ShapeRecognMesh *New();
+ std::size_t getHeapMemorySizeWithoutChildren() const;
+ std::vector<const BigMemoryObject *> getDirectChildrenWithNull() const;
+
+ void save(const std::string &outputFile, bool writeFromScratch = true) const;
+
+ // Node properties
+ MEDCoupling::MEDCouplingFieldDouble *getNodeK1() const;
+ MEDCoupling::MEDCouplingFieldDouble *getNodeK2() const;
+ MEDCoupling::MEDCouplingFieldDouble *getNodePrimitiveType() const;
+ MEDCoupling::MEDCouplingFieldDouble *getNodeNormal() const;
+
+ // Area properties
+ MEDCoupling::MEDCouplingFieldDouble *getAreaId() const;
+ MEDCoupling::MEDCouplingFieldDouble *getAreaPrimitiveType() const;
+ MEDCoupling::MEDCouplingFieldDouble *getAreaNormal() const;
+ MEDCoupling::MEDCouplingFieldDouble *getMinorRadius() const;
+ MEDCoupling::MEDCouplingFieldDouble *getRadius() const;
+ MEDCoupling::MEDCouplingFieldDouble *getAngle() const;
+ MEDCoupling::MEDCouplingFieldDouble *getCenter() const;
+ MEDCoupling::MEDCouplingFieldDouble *getAxis() const;
+ MEDCoupling::MEDCouplingFieldDouble *getApex() const;
+
+ protected:
+ ShapeRecognMesh();
+ ~ShapeRecognMesh();
+
+ private:
+ MEDCoupling::MEDCouplingFieldDouble *nodeK1;
+ MEDCoupling::MEDCouplingFieldDouble *nodeK2;
+ MEDCoupling::MEDCouplingFieldDouble *nodePrimitiveType;
+ MEDCoupling::MEDCouplingFieldDouble *nodeNormal;
+ MEDCoupling::MEDCouplingFieldDouble *areaId;
+ MEDCoupling::MEDCouplingFieldDouble *areaPrimitiveType;
+ MEDCoupling::MEDCouplingFieldDouble *areaNormal;
+ MEDCoupling::MEDCouplingFieldDouble *minorRadius;
+ MEDCoupling::MEDCouplingFieldDouble *radius;
+ MEDCoupling::MEDCouplingFieldDouble *angle;
+ MEDCoupling::MEDCouplingFieldDouble *center;
+ MEDCoupling::MEDCouplingFieldDouble *axis;
+ MEDCoupling::MEDCouplingFieldDouble *apex;
+ };
+};
+
+#endif // __SHAPERECOGNMESH_HXX__
--- /dev/null
+#include "ShapeRecognMeshBuilder.hxx"
+
+#include "NodesBuilder.hxx"
+#include "AreasBuilder.hxx"
+#include "MEDLoader.hxx"
+#include "ShapeRecognMesh.hxx"
+#include "MEDCouplingFieldDouble.hxx"
+
+using namespace MEDCoupling;
+
+ShapeRecognMeshBuilder::ShapeRecognMeshBuilder(const std::string &fileName, int meshDimRelToMax)
+{
+ mesh = ReadUMeshFromFile(fileName, meshDimRelToMax);
+ if (mesh->getMeshDimension() != 2)
+ throw INTERP_KERNEL::Exception("Expect a mesh with a dimension equal to 2");
+ if (mesh->getNumberOfCellsWithType(INTERP_KERNEL::NORM_TRI3) != mesh->getNumberOfCells())
+ throw INTERP_KERNEL::Exception("Expect a mesh containing exclusively triangular cells");
+}
+
+ShapeRecognMeshBuilder::~ShapeRecognMeshBuilder()
+{
+ if (areas != nullptr)
+ delete areas;
+ if (nodes != nullptr)
+ delete nodes;
+ mesh->decrRef();
+}
+
+ShapeRecognMesh *ShapeRecognMeshBuilder::recognize()
+{
+ mesh->incrRef();
+ NodesBuilder nodesBuilder(mesh);
+ nodes = nodesBuilder.build();
+ AreasBuilder areasBuilder(nodes);
+ areasBuilder.build();
+ areas = areasBuilder.getAreas();
+ MCAuto<ShapeRecognMesh> recognMesh = ShapeRecognMesh::New();
+ recognMesh->nodeK1 = buildNodeK1();
+ recognMesh->nodeK2 = buildNodeK2();
+ recognMesh->nodePrimitiveType = buildNodePrimitiveType();
+ recognMesh->nodeNormal = buildNodeNormal();
+ recognMesh->areaId = buildAreaId();
+ recognMesh->areaPrimitiveType = buildAreaPrimitiveType();
+ recognMesh->areaNormal = buildAreaNormal();
+ recognMesh->minorRadius = buildMinorRadius();
+ recognMesh->radius = buildRadius();
+ recognMesh->angle = buildAngle();
+ recognMesh->center = buildCenter();
+ recognMesh->axis = buildAxis();
+ recognMesh->apex = buildApex();
+ return recognMesh.retn();
+}
+
+const Nodes *ShapeRecognMeshBuilder::getNodes() const
+{
+ return nodes;
+}
+
+const Areas *ShapeRecognMeshBuilder::getAreas() const
+{
+ return areas;
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildNodeK1() const
+{
+ if (nodes == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ return buildField("K1 (Node)", 1, nodes->getK1());
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildNodeK2() const
+{
+ if (nodes == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ return buildField("K2 (Node)", 1, nodes->getK2());
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildNodePrimitiveType() const
+{
+ if (nodes == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ return buildField("Primitive Type (Node)", 1, nodes->getPrimitiveType());
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildNodeNormal() const
+{
+ if (nodes == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ return buildField("Normal (Node)", 3, nodes->getNormals());
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildAreaId() const
+{
+ if (areas == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ return buildField("Area Id", 1, areas->getAreaIdByNodes());
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildAreaPrimitiveType() const
+{
+ if (areas == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ double *values = buildAreaArray([](Areas *areas, mcIdType areaId) -> double
+ { return (double)areas->getPrimitiveType(areaId); });
+ return buildField("Primitive Type (Area)", 1, values);
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildAreaNormal() const
+{
+ if (areas == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ double *values = buildArea3DArray([](Areas *areas, mcIdType areaId) -> const std::array<double, 3> &
+ { return areas->getNormal(areaId); });
+ return buildField("Normal (Area)", 3, values);
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildMinorRadius() const
+{
+ if (areas == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ double *values = buildAreaArray([](Areas *areas, mcIdType areaId) -> double
+ { return areas->getMinorRadius(areaId); });
+ return buildField("Minor Radius (Area)", 1, values);
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildRadius() const
+{
+ if (areas == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ double *values = buildAreaArray([](Areas *areas, mcIdType areaId) -> double
+ { return areas->getRadius(areaId); });
+ return buildField("Radius (Area)", 1, values);
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildAngle() const
+{
+ if (areas == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ double *values = buildAreaArray([](Areas *areas, mcIdType areaId) -> double
+ { return areas->getAngle(areaId); });
+ return buildField("Angle (Area)", 1, values);
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildCenter() const
+{
+ if (areas == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ double *values = buildArea3DArray([](Areas *areas, mcIdType areaId) -> const std::array<double, 3> &
+ { return areas->getCenter(areaId); });
+ return buildField("Center (Area)", 3, values);
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildAxis() const
+{
+ if (areas == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ double *values = buildArea3DArray([](Areas *areas, mcIdType areaId) -> const std::array<double, 3> &
+ { return areas->getAxis(areaId); });
+ return buildField("Axis (Area)", 3, values);
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildApex() const
+{
+ if (areas == nullptr)
+ throw INTERP_KERNEL::Exception("recognize must be called before building any fields");
+ double *values = buildArea3DArray([](Areas *areas, mcIdType areaId) -> const std::array<double, 3> &
+ { return areas->getApex(areaId); });
+ return buildField("Apex (Area)", 3, values);
+}
+
+template <typename T>
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildField(
+ const std::string &name,
+ size_t nbOfCompo,
+ const std::vector<T> &values) const
+{
+ DataArrayDouble *data = DataArrayDouble::New();
+ data->setName(name);
+ data->alloc(nodes->getNbNodes(), nbOfCompo);
+ std::copy(values.begin(), values.end(), data->getPointer());
+ data->declareAsNew();
+ return buildField(name, nbOfCompo, data);
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildField(
+ const std::string &name,
+ size_t nbOfCompo,
+ double *values) const
+{
+ DataArrayDouble *data = DataArrayDouble::New();
+ data->setName(name);
+ data->useArray(
+ values,
+ true,
+ MEDCoupling::DeallocType::CPP_DEALLOC,
+ nodes->getNbNodes(),
+ nbOfCompo);
+ return buildField(name, nbOfCompo, data);
+}
+
+MEDCouplingFieldDouble *ShapeRecognMeshBuilder::buildField(
+ const std::string &name,
+ size_t nbOfCompo,
+ DataArrayDouble *data) const
+{
+ MEDCouplingFieldDouble *field = MEDCouplingFieldDouble::New(ON_NODES);
+ field->setName(name);
+ field->setMesh(mesh);
+ if (nbOfCompo == 3)
+ {
+ data->setInfoOnComponent(0, "X");
+ data->setInfoOnComponent(1, "Y");
+ data->setInfoOnComponent(2, "Z");
+ }
+ field->setArray(data);
+ data->decrRef();
+ return field;
+}
+
+double *ShapeRecognMeshBuilder::buildArea3DArray(
+ const std::array<double, 3> &(*areaFunc)(Areas *, mcIdType)) const
+{
+ double *values = new double[3 * nodes->getNbNodes()];
+ const std::vector<mcIdType> &areaIdByNodes = areas->getAreaIdByNodes();
+ for (size_t nodeId = 0; nodeId < areaIdByNodes.size(); ++nodeId)
+ {
+ mcIdType areaId = areaIdByNodes[nodeId];
+ if (areaId != -1)
+ {
+ const std::array<double, 3> &areaValues = areaFunc(areas, areaId);
+ values[3 * nodeId] = areaValues[0];
+ values[3 * nodeId + 1] = areaValues[1];
+ values[3 * nodeId + 2] = areaValues[2];
+ }
+ }
+ return values;
+}
+
+double *ShapeRecognMeshBuilder::buildAreaArray(double (*areaFunc)(Areas *, mcIdType)) const
+{
+ const std::vector<mcIdType> &areaIdByNodes = areas->getAreaIdByNodes();
+ double *values = new double[nodes->getNbNodes()];
+ for (size_t nodeId = 0; nodeId < areaIdByNodes.size(); ++nodeId)
+ {
+ mcIdType areaId = areaIdByNodes[nodeId];
+ if (areaId != -1)
+ values[nodeId] = areaFunc(areas, areaId);
+ }
+ return values;
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __SHAPERECOGNMESHBUILDER_HXX__
+#define __SHAPERECOGNMESHBUILDER_HXX__
+
+#include <string>
+
+#include "MEDCouplingUMesh.hxx"
+#include "MEDCouplingFieldDouble.hxx"
+
+namespace MEDCoupling
+{
+ class Nodes;
+ class Areas;
+ class ShapeRecognMesh;
+
+ class ShapeRecognMeshBuilder
+ {
+ public:
+ ShapeRecognMeshBuilder(const std::string &fileName, int meshDimRelToMax = 0);
+ ~ShapeRecognMeshBuilder();
+
+ const Nodes *getNodes() const;
+ const Areas *getAreas() const;
+
+ ShapeRecognMesh *recognize();
+
+ private:
+ // Node properties
+ MEDCoupling::MEDCouplingFieldDouble *buildNodeK1() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildNodeK2() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildNodePrimitiveType() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildNodeNormal() const;
+
+ // Area properties
+ MEDCoupling::MEDCouplingFieldDouble *buildAreaId() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildAreaPrimitiveType() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildAreaNormal() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildMinorRadius() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildRadius() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildAngle() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildCenter() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildAxis() const;
+ MEDCoupling::MEDCouplingFieldDouble *buildApex() const;
+
+ template <typename T>
+ MEDCouplingFieldDouble *buildField(
+ const std::string &name,
+ size_t nbOfCompo,
+ const std::vector<T> &values) const;
+ MEDCouplingFieldDouble *buildField(
+ const std::string &name,
+ size_t nbOfCompo,
+ double *values) const;
+ MEDCouplingFieldDouble *buildField(
+ const std::string &name,
+ size_t nbOfCompo,
+ DataArrayDouble *values) const;
+ double *buildArea3DArray(const std::array<double, 3> &(*areaFunc)(Areas *, mcIdType)) const;
+ double *buildAreaArray(double (*areaFunc)(Areas *, mcIdType)) const;
+
+ const MEDCouplingUMesh *mesh;
+ Nodes *nodes = nullptr;
+ Areas *areas = nullptr;
+ };
+};
+
+#endif // __SHAPERECOGNMESHBUILDER_HXX__
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __SHAPERECOGNCONSTANTS_HXX__
+#define __SHAPERECOGNCONSTANTS_HXX__
+
+namespace MEDCoupling
+{
+ // Nodes
+ constexpr double EPSILON_PRIMITIVE = 0.005;
+ // Areas
+ // - Match
+ constexpr double TOL_MATCH_CYLINDER = 0.05;
+ constexpr double TOL_MATCH_SPHERE = 0.05;
+ // - Relative distance
+ constexpr double DELTA_PLANE = 0.05;
+ constexpr double DELTA_SPHERE = 0.05;
+ constexpr double DELTA_CYLINDER = 0.05;
+ constexpr double DELTA_CONE = 0.05;
+ // - Invalid Zones
+ constexpr double THETA_MAX_CYLINDER = 0.05;
+ // - Thresholds
+ constexpr int THRESHOLD_MIN_NB_NODES = 5;
+ constexpr int THRESHOLD_MAX_NB_AREAS = 500;
+};
+
+#endif //__SHAPERECOGNCONSTANTS_HXX__
\ No newline at end of file
--- /dev/null
+# Copyright (C) 2012-2024 CEA, EDF
+#
+# 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 2.1 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, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+
+FIND_PACKAGE(SWIG REQUIRED)
+INCLUDE(${SWIG_USE_FILE})
+
+ADD_DEFINITIONS(${PYTHON_DEFINITIONS})
+
+SET_SOURCE_FILES_PROPERTIES(ShapeRecogn.i PROPERTIES CPLUSPLUS ON)
+IF ("${PYTHON_VERSION_MAJOR}" STREQUAL "3")
+ SET_SOURCE_FILES_PROPERTIES(ShapeRecogn.i PROPERTIES SWIG_FLAGS "-py3")
+ELSE()
+ SET_SOURCE_FILES_PROPERTIES(ShapeRecogn.i PROPERTIES SWIG_DEFINITIONS "-shadow")
+ENDIF()
+
+SET(SWIG_MODULE_ShapeRecogn_EXTRA_FLAGS "")
+IF(MEDCOUPLING_USE_64BIT_IDS)
+ STRING(APPEND SWIG_MODULE_ShapeRecogn_EXTRA_FLAGS "-DMEDCOUPLING_USE_64BIT_IDS")
+ENDIF(MEDCOUPLING_USE_64BIT_IDS)
+
+SET (ShapeRecogn_SWIG_DPYS_FILES
+ ShapeRecogn.i
+)
+
+INCLUDE_DIRECTORIES(
+ ${PYTHON_INCLUDE_DIRS}
+ ${MEDFILE_INCLUDE_DIRS}
+ ${HDF5_INCLUDE_DIRS}
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ ${CMAKE_CURRENT_SOURCE_DIR}/..
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../MEDLoader
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../MEDCoupling_Swig
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../MEDCoupling
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../INTERP_KERNEL
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../INTERP_KERNEL/Bases
+)
+
+SET (SWIG_MODULE_ShapeRecogn_EXTRA_DEPS ${ShapeRecogn_SWIG_DPYS_FILES}
+ ${shaperecogn_HEADERS_HXX}
+ ${medcoupling_HEADERS_HXX} ${medcoupling_HEADERS_TXX}
+ ${interpkernel_HEADERS_HXX} ${interpkernel_HEADERS_TXX})
+
+IF(WIN32)
+ SET_PROPERTY(SOURCE ShapeRecogn.i PROPERTY COMPILE_DEFINITIONS WIN32)
+ENDIF()
+
+IF(${CMAKE_VERSION} VERSION_LESS "3.8.0")
+ SWIG_ADD_MODULE(ShapeRecogn python ShapeRecogn.i)
+ELSE()
+ SWIG_ADD_LIBRARY(ShapeRecogn LANGUAGE python SOURCES ShapeRecogn.i)
+ENDIF()
+
+SWIG_LINK_LIBRARIES(ShapeRecogn ${PYTHON_LIBRARIES} ${PLATFORM_LIBS} shaperecogn medloader medcouplingcpp)
+SWIG_CHECK_GENERATION(ShapeRecogn)
+IF(WIN32)
+ SET_TARGET_PROPERTIES(_ShapeRecogn PROPERTIES DEBUG_OUTPUT_NAME _ShapeRecogn_d)
+ # To increase the size of the .obj file on Windows because ShapeRecognPYTHON_wrap.cxx, generated by SWIG, is too big
+ TARGET_COMPILE_OPTIONS(_ShapeRecogn PRIVATE /bigobj)
+ENDIF(WIN32)
+
+INSTALL(TARGETS _ShapeRecogn DESTINATION ${MEDCOUPLING_INSTALL_PYTHON})
+INSTALL(FILES ShapeRecogn.i DESTINATION ${MEDCOUPLING_INSTALL_HEADERS})
+
+SALOME_INSTALL_SCRIPTS(
+ ${CMAKE_CURRENT_BINARY_DIR}/ShapeRecogn.py
+ ${MEDCOUPLING_INSTALL_PYTHON}
+ EXTRA_DPYS "${SWIG_MODULE_ShapeRecogn_REAL_NAME}")
--- /dev/null
+%module ShapeRecogn
+
+%include "std_string.i"
+%include "MEDCouplingCommon.i"
+
+%{
+#include "ShapeRecognMesh.hxx"
+#include "ShapeRecognMeshBuilder.hxx"
+using namespace MEDCoupling;
+%}
+
+%ignore getAreas() const;
+%ignore getNodes() const;
+%include "ShapeRecognMesh.hxx"
+%include "ShapeRecognMeshBuilder.hxx"
--- /dev/null
+# Copyright (C) 2012-2024 CEA, EDF
+#
+# 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 2.1 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, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+
+INCLUDE_DIRECTORIES(
+ ${CPPUNIT_INCLUDE_DIRS}
+ ${HDF5_INCLUDE_DIRS}
+ ${MEDFILE_INCLUDE_DIRS}
+ ${CMAKE_CURRENT_SOURCE_DIR}/..
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../MEDLoader
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../MEDCoupling
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../INTERP_KERNEL
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../INTERP_KERNEL/Bases
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../INTERP_KERNELTest # For common CppUnitTest.hxx file and TestIKUtils.hxx
+ )
+
+SET(TestShapeRecogn_SOURCES
+ TestShapeRecogn.cxx
+ MathOpsTest.cxx
+ PlaneTest.cxx
+ CylinderTest.cxx
+ ConeTest.cxx
+ SphereTest.cxx
+ TorusTest.cxx
+)
+
+SALOME_ACCUMULATE_ENVIRONMENT(MEDCOUPLING_RESOURCE_DIR "${CMAKE_BINARY_DIR}/resources")
+SALOME_GENERATE_TESTS_ENVIRONMENT(tests_env)
+
+ADD_EXECUTABLE(TestShapeRecogn ${TestShapeRecogn_SOURCES})
+TARGET_LINK_LIBRARIES(TestShapeRecogn shaperecogn InterpKernelTestUtils ${CPPUNIT_LIBRARIES} ${PLATFORM_LIBS})
+
+INSTALL(TARGETS TestShapeRecogn DESTINATION ${MEDCOUPLING_INSTALL_BINS})
+
+SET(BASE_TESTS TestShapeRecogn)
+
+FOREACH(test ${BASE_TESTS})
+ ADD_TEST(NAME ${test}
+ COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MCTestLauncher.py ${CMAKE_CURRENT_BINARY_DIR}/${test})
+ SET_TESTS_PROPERTIES(${test} PROPERTIES ENVIRONMENT "${tests_env}")
+ENDFOREACH()
+
+# Application tests
+
+SET(TEST_INSTALL_DIRECTORY ${MEDCOUPLING_INSTALL_TESTS}/ShapeRecogn)
+INSTALL(TARGETS TestShapeRecogn DESTINATION ${TEST_INSTALL_DIRECTORY})
+
+INSTALL(FILES CTestTestfileInstall.cmake
+ DESTINATION ${TEST_INSTALL_DIRECTORY}
+ RENAME CTestTestfile.cmake)
--- /dev/null
+# Copyright (C) 2015-2024 CEA, EDF
+#
+# 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 2.1 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, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+
+SET(TEST_NAMES
+ TestShapeRecogn
+)
+
+FOREACH(tfile ${TEST_NAMES})
+ SET(TEST_NAME ${COMPONENT_NAME}_${tfile})
+ ADD_TEST(${TEST_NAME} python3 MCTestLauncher.py ${tfile})
+ SET_TESTS_PROPERTIES(${TEST_NAME} PROPERTIES
+ LABELS "${COMPONENT_NAME}"
+ TIMEOUT ${TIMEOUT}
+ )
+ENDFOREACH()
--- /dev/null
+#include "ConeTest.hxx"
+
+#include "ShapeRecognMeshBuilder.hxx"
+#include "Areas.hxx"
+#include "MathOps.hxx"
+#include "TestInterpKernelUtils.hxx" // getResourceFile()
+
+using namespace MEDCoupling;
+
+void ConeTest::setUp()
+{
+ std::string file = INTERP_TEST::getResourceFile("ShapeRecognCone.med", 3);
+ srMesh = new ShapeRecognMeshBuilder(file);
+ srMesh->recognize();
+ areas = srMesh->getAreas();
+}
+
+void ConeTest::tearDown()
+{
+ if (srMesh != 0)
+ delete srMesh;
+ areas = 0;
+}
+
+void ConeTest::testNumberOfAreas()
+{
+ CPPUNIT_ASSERT_EQUAL(3, (int)areas->getNumberOfAreas());
+}
+
+void ConeTest::testComputePlaneProperties()
+{
+ const Nodes *nodes = srMesh->getNodes();
+ Areas areas(nodes);
+ mcIdType areaId = areas.addArea();
+ std::vector<mcIdType> nodeIds{560, 561, 562, 563, 564};
+ // Check the coordinates
+ std::vector<double> coordsRef{
+ -0.13274028,
+ -1.19035036,
+ 7.0,
+ 0.62579096,
+ -0.97231524,
+ 7.0,
+ -0.75942852,
+ -0.99273623,
+ 7.0,
+ 0.76454326,
+ 0.85588873,
+ 7.0,
+ 1.12546475,
+ -0.04924272,
+ 7.0,
+ };
+ for (size_t i = 0; i < nodeIds.size(); ++i)
+ {
+ mcIdType nodeId = nodeIds[i];
+ areas.addNode(areaId, nodeId);
+ std::array<double, 3> coords = nodes->getCoordinates(nodeId);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i], coords[0], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 1], coords[1], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 2], coords[2], 1E-6);
+ }
+ areas.computeProperties(areaId);
+ // Check normal
+ std::array<double, 3> normal = areas.getNormal(areaId);
+ std::array<double, 3> normalRef{0.0, 0.0, 1.0};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(normalRef[i], normal[i], 1E-6);
+}
+
+void ConeTest::testComputeCylinderProperties()
+{
+ const Nodes *nodes = srMesh->getNodes();
+ Areas areas(nodes);
+ mcIdType areaId = areas.addArea();
+ // check coordinates
+ std::vector<double> coordsRef{
+ 8.58402821e+00,
+ -2.10248053e-15,
+ 1.41597179e+00,
+ 7.87604232e+00,
+ -1.92907400e-15,
+ 2.12395768e+00,
+ 7.16805642e+00,
+ -1.75566747e-15,
+ 2.83194358e+00,
+ };
+ // Add nodes
+ std::vector<mcIdType> nodeIds{26, 27, 28};
+ for (size_t i = 0; i < nodeIds.size(); ++i)
+ {
+ mcIdType nodeId = nodeIds[i];
+ areas.addNode(areaId, nodeId);
+ std::array<double, 3> coords = nodes->getCoordinates(nodeId);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i], coords[0], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 1], coords[1], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 2], coords[2], 1E-6);
+ }
+ areas.computeProperties(areaId);
+ // Check radius
+ double radius = areas.getRadius(areaId);
+ double radiusRef = 11.016809459292505;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(radiusRef, radius, 1E-6);
+ // Check axis
+ std::array<double, 3> axis = areas.getAxis(areaId);
+ std::array<double, 3> axisRef{-2.85185021e-02, 2.47000552e-05, 9.99593264e-01};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(axisRef[i], axis[i], 1E-6);
+ // Check axis point
+ std::array<double, 3> axisPoint = areas.getAxisPoint(areaId);
+ std::array<double, 3> axisPointRef{9.55163389e-02, -2.27195412e-05, -5.67562579e+00};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(axisPointRef[i], axisPoint[i], 1E-6);
+}
+
+void ConeTest::testComputeConeProperties()
+{
+ std::srand(0);
+ const Nodes *nodes = srMesh->getNodes();
+ Areas areas(nodes);
+ mcIdType areaId = areas.addArea(PrimitiveType::Cone);
+ std::vector<double> coordsRef{
+ 8.58402821e+00,
+ -2.10248053e-15,
+ 1.41597179e+00,
+ 7.87604232e+00,
+ -1.92907400e-15,
+ 2.12395768e+00,
+ 7.16805642e+00,
+ -1.75566747e-15,
+ 2.83194358e+00,
+ 6.46007053e+00,
+ -1.58226094e-15,
+ 3.53992947e+00};
+ // Add nodes
+ std::vector<mcIdType> nodeIds{26, 27, 28, 29};
+ for (size_t i = 0; i < nodeIds.size(); ++i)
+ {
+ mcIdType nodeId = nodeIds[i];
+ areas.addNode(areaId, nodeId);
+ std::array<double, 3> coords = nodes->getCoordinates(nodeId);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i], coords[0], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 1], coords[1], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(coordsRef[3 * i + 2], coords[2], 1E-6);
+ }
+ areas.computeProperties(areaId);
+ // Radius
+ double radiusRef = 10.546213172989718;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(radiusRef, areas.getRadius(areaId), 1E-6);
+ // Angle
+ double angleRef = 0.7567719849399294;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(angleRef, areas.getAngle(areaId), 1E-6);
+ // Axis
+ std::array<double, 3> axisRef{-2.99093478e-02, 1.48776702e-05, 9.99552615e-01};
+ std::array<double, 3> axis = areas.getAxis(areaId);
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(axisRef[i], axis[i], 1E-6);
+ // Axis Point
+ std::array<double, 3> axisPointRef{
+ 7.43242419e-02, -1.70396559e-05, -4.98890939e+00};
+ std::array<double, 3> axisPoint = areas.getAxisPoint(areaId);
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(axisPointRef[i], axisPoint[i], 1E-6);
+ // Apex
+ std::array<double, 3> apexRef{
+ -3.87870047e-01, 1.98282729e-04, 1.03928162e+01};
+ std::array<double, 3> apex = areas.getApex(areaId);
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(apexRef[i], apex[i], 1E-6);
+}
+
+void ConeTest::testFirstArea()
+{
+ // primitive type
+ CPPUNIT_ASSERT_EQUAL(PrimitiveType::Plane, areas->getPrimitiveType(0));
+ // node ids
+ std::vector<mcIdType> nodeIdsRef{
+ 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559,
+ 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570,
+ 571};
+ std::vector<mcIdType> nodeIds = areas->getNodeIds(0);
+ std::sort(nodeIds.begin(), nodeIds.end());
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size());
+ for (size_t i = 0; i < nodeIds.size(); ++i)
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]);
+ // normal
+ std::array<double, 3> normal = areas->getNormal(0);
+ std::array<double, 3> normalRef{0.0, 0.0, 1.0};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(normalRef[i], normal[i], 1E-6);
+}
+
+void ConeTest::testSecondArea()
+{
+ // primitive type
+ CPPUNIT_ASSERT_EQUAL(PrimitiveType::Plane, areas->getPrimitiveType(0));
+ // node ids
+ std::vector<mcIdType> nodeIdsRef = {
+ 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583,
+ 584, 585, 586, 587, 588, 589, 590, 591, 592, 593, 594, 595,
+ 596, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606, 607,
+ 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619,
+ 620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631,
+ 632, 633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643,
+ 644, 645, 646, 647, 648, 649, 650, 651, 652, 653, 654, 655,
+ 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667,
+ 668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679,
+ 680, 681, 682, 683, 684, 685, 686, 687, 688, 689, 690, 691,
+ 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703,
+ 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715,
+ 716, 717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727,
+ 728, 729, 730, 731, 732, 733, 734, 735, 736, 737, 738, 739,
+ 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751,
+ 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762, 763,
+ 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775,
+ 776, 777, 778, 779, 780, 781, 782, 783, 784, 785, 786, 787,
+ 788, 789, 790, 791, 792, 793, 794, 795, 796, 797, 798, 799,
+ 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 810, 811,
+ 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823,
+ 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835,
+ 836, 837, 838, 839, 840, 841, 842, 843, 844, 845, 846, 847,
+ 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 858, 859,
+ 860, 861, 862, 863, 864, 865, 866, 867, 868, 869, 870, 871,
+ 872, 873, 874, 875, 876, 877, 878, 879, 880, 881, 882, 883,
+ 884, 885, 886, 887, 888, 889, 890, 891, 892, 893, 894, 895,
+ 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907,
+ 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919,
+ 920, 921, 922, 923, 924, 925, 926, 927, 928, 929, 930, 931,
+ 932, 933, 934, 935, 936, 937, 938, 939, 940, 941, 942, 943,
+ 944, 945, 946, 947, 948, 949, 950, 951, 952, 953, 954, 955,
+ 956};
+ std::vector<mcIdType> nodeIds = areas->getNodeIds(1);
+ std::sort(nodeIds.begin(), nodeIds.end());
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size());
+ for (size_t i = 0; i < nodeIds.size(); ++i)
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]);
+ // normal
+ std::array<double, 3> normal = areas->getNormal(1);
+ std::array<double, 3> normalRef = {0.0, 0.0, -1.0};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(normalRef[i], normal[i], 1E-6);
+}
+
+void ConeTest::testThirdArea()
+{
+ // primitive type
+ CPPUNIT_ASSERT_EQUAL(PrimitiveType::Cone, areas->getPrimitiveType(2));
+ // node ids
+ std::vector<mcIdType> nodeIdsRef{
+ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
+ 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
+ 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43,
+ 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57,
+ 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,
+ 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85,
+ 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,
+ 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110,
+ 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121,
+ 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132,
+ 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
+ 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154,
+ 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165,
+ 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176,
+ 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187,
+ 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198,
+ 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209,
+ 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220,
+ 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231,
+ 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242,
+ 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253,
+ 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264,
+ 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275,
+ 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286,
+ 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297,
+ 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308,
+ 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319,
+ 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330,
+ 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341,
+ 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352,
+ 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363,
+ 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374,
+ 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385,
+ 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396,
+ 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407,
+ 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418,
+ 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429,
+ 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440,
+ 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451,
+ 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462,
+ 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473,
+ 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484,
+ 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495,
+ 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506,
+ 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517,
+ 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528,
+ 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539,
+ 540, 541, 542, 543, 544, 545, 546, 547, 548};
+ std::vector<mcIdType> nodeIds = areas->getNodeIds(2);
+ std::sort(nodeIds.begin(), nodeIds.end());
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size());
+ for (size_t i = 0; i < nodeIds.size(); ++i)
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]);
+ // radius
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(10.0, areas->getRadius(2), 5E-1);
+ // angle
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(M_PI_4, areas->getAngle(2), 1E-2);
+ // axis
+ std::array<double, 3> axis = areas->getAxis(2);
+ std::array<double, 3> axisRef{0.0, 0.0, 1.0};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(axisRef[i], axis[i], 1E-2);
+ // // axis point
+ // std::array<double, 3> axisPoint = areas->getAxisPoint(2);
+ // std::array<double, 3> axisPointRef{
+ // -6.50039462e-03, -5.76665576e-05, -3.99198765e+00};
+ // for (size_t i = 0; i < 3; ++i)
+ // CPPUNIT_ASSERT_DOUBLES_EQUAL(
+ // axisPointRef[i], axisPoint[i], 1E-2);
+ // apex
+ std::array<double, 3> apex = areas->getApex(2);
+ std::array<double, 3> apexRef{0.0, 0.0, 10.0};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(apexRef[i], apex[i], 1E-1);
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __CONETEST_HXX__
+#define __CONETEST_HXX__
+
+#include <cppunit/TestFixture.h>
+#include <cppunit/extensions/HelperMacros.h>
+
+namespace MEDCoupling
+{
+ class ShapeRecognMeshBuilder;
+ class Areas;
+
+ class ConeTest : public CppUnit::TestFixture
+ {
+ CPPUNIT_TEST_SUITE(ConeTest);
+ CPPUNIT_TEST(testNumberOfAreas);
+ CPPUNIT_TEST(testComputePlaneProperties);
+ CPPUNIT_TEST(testComputeCylinderProperties);
+ CPPUNIT_TEST(testComputeConeProperties);
+ CPPUNIT_TEST(testFirstArea);
+ CPPUNIT_TEST(testSecondArea);
+ CPPUNIT_TEST(testThirdArea);
+ CPPUNIT_TEST_SUITE_END();
+
+ public:
+ void setUp() override;
+ void tearDown() override;
+
+ void testComputePlaneProperties();
+ void testComputeCylinderProperties();
+ void testComputeConeProperties();
+
+ void testNumberOfAreas();
+ void testFirstArea();
+ void testSecondArea();
+ void testThirdArea();
+
+ private:
+ ShapeRecognMeshBuilder *srMesh = 0;
+ const Areas *areas;
+ };
+};
+
+#endif // __CONETEST_HXX__
--- /dev/null
+#include "CylinderTest.hxx"
+
+#include "ShapeRecognMeshBuilder.hxx"
+#include "Areas.hxx"
+#include "MathOps.hxx"
+#include "TestInterpKernelUtils.hxx" // getResourceFile()
+
+using namespace MEDCoupling;
+
+void CylinderTest::setUp()
+{
+ std::string file = INTERP_TEST::getResourceFile("ShapeRecognCylinder.med", 3);
+ srMesh = new ShapeRecognMeshBuilder(file);
+ srMesh->recognize();
+ areas = srMesh->getAreas();
+}
+
+void CylinderTest::tearDown()
+{
+ if (srMesh != 0)
+ delete srMesh;
+ areas = 0;
+}
+
+void CylinderTest::testNumberOfAreas()
+{
+ CPPUNIT_ASSERT_EQUAL(3, (int)areas->getNumberOfAreas());
+}
+
+void CylinderTest::testFirstArea()
+{
+ // primitive type
+ CPPUNIT_ASSERT_EQUAL(PrimitiveType::Cylinder, areas->getPrimitiveType(0));
+ // node ids
+ std::vector<mcIdType> nodeIdsRef{
+ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
+ 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
+ 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46,
+ 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61,
+ 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
+ 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91,
+ 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104,
+ 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
+ 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126,
+ 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137,
+ 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148,
+ 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159,
+ 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170,
+ 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181,
+ 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192,
+ 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203,
+ 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214,
+ 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225,
+ 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236,
+ 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247,
+ 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258,
+ 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269,
+ 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280,
+ 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291,
+ 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302,
+ 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313,
+ 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324,
+ 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335,
+ 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346,
+ 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357,
+ 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368,
+ 369};
+ std::vector<mcIdType> nodeIds = areas->getNodeIds(0);
+ std::sort(nodeIds.begin(), nodeIds.end());
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size());
+ for (size_t i = 0; i < nodeIds.size(); ++i)
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]);
+ // radius
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(4.993494657779537, areas->getRadius(0), 1E-2);
+ // axis
+ std::array<double, 3> axis = areas->getAxis(0);
+ std::array<double, 3> axisRef{7.66631075e-04, -1.59966800e-04, 9.99999693e-01};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(axisRef[i], axis[i], 1E-4);
+ // axis point
+ std::array<double, 3> axisPoint = areas->getAxisPoint(0);
+ std::array<double, 3> axisPointRef{3.78927537e-03, -2.03409415e-03, 5.03355802e+00};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(
+ axisPointRef[i], axisPoint[i], 1E-2);
+}
+
+void CylinderTest::testSecondArea()
+{
+ // primitive type
+ CPPUNIT_ASSERT_EQUAL(PrimitiveType::Plane, areas->getPrimitiveType(1));
+ // node ids
+ std::vector<mcIdType> nodeIdsRef{
+ 370, 371, 372, 373, 374, 375, 376, 377, 378, 379,
+ 380, 381, 382, 383, 384, 385, 386, 387, 388, 389,
+ 390, 391, 392, 393, 394, 395, 396, 397, 398, 399,
+ 400, 401, 402, 403, 404, 405, 406, 407, 408, 409,
+ 410, 411, 412, 413, 414, 415, 416, 417, 418, 419,
+ 420, 421, 422, 423, 424, 425, 426, 427, 428, 429,
+ 430, 431, 432, 433, 434, 435, 436, 437, 438, 439,
+ 440, 441, 442};
+ std::vector<mcIdType> nodeIds = areas->getNodeIds(1);
+ std::sort(nodeIds.begin(), nodeIds.end());
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size());
+ for (size_t i = 0; i < nodeIds.size(); ++i)
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]);
+ // normal
+ std::array<double, 3> normal = areas->getNormal(1);
+ std::array<double, 3> normalRef{0.0, 0.0, 1.0};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(normalRef[i], normal[i], 1E-6);
+}
+
+void CylinderTest::testThirdArea()
+{
+ // primitive type
+ CPPUNIT_ASSERT_EQUAL(PrimitiveType::Plane, areas->getPrimitiveType(2));
+ // node ids
+ std::vector<mcIdType> nodeIdsRef{
+ 443, 444, 445, 446, 447, 448, 449, 450, 451, 452,
+ 453, 454, 455, 456, 457, 458, 459, 460, 461, 462,
+ 463, 464, 465, 466, 467, 468, 469, 470, 471, 472,
+ 473, 474, 475, 476, 477, 478, 479, 480, 481, 482,
+ 483, 484, 485, 486, 487, 488, 489, 490, 491, 492,
+ 493, 494, 495, 496, 497, 498, 499};
+ std::vector<mcIdType> nodeIds = areas->getNodeIds(2);
+ std::sort(nodeIds.begin(), nodeIds.end());
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef.size(), nodeIds.size());
+ for (size_t i = 0; i < nodeIds.size(); ++i)
+ CPPUNIT_ASSERT_EQUAL(nodeIdsRef[i], nodeIds[i]);
+ // normal
+ std::array<double, 3> normal = areas->getNormal(2);
+ std::array<double, 3> normalRef{0.0, 0.0, -1.0};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(normalRef[i], normal[i], 1E-6);
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __CYLINDERTEST_HXX__
+#define __CYLINDERTEST_HXX__
+
+#include <cppunit/TestFixture.h>
+#include <cppunit/extensions/HelperMacros.h>
+
+namespace MEDCoupling
+{
+ class ShapeRecognMeshBuilder;
+ class Areas;
+
+ class CylinderTest : public CppUnit::TestFixture
+ {
+ CPPUNIT_TEST_SUITE(CylinderTest);
+ CPPUNIT_TEST(testNumberOfAreas);
+ CPPUNIT_TEST(testFirstArea);
+ CPPUNIT_TEST(testSecondArea);
+ CPPUNIT_TEST(testThirdArea);
+ CPPUNIT_TEST_SUITE_END();
+
+ public:
+ void setUp() override;
+ void tearDown() override;
+
+ void testNumberOfAreas();
+ void testFirstArea();
+ void testSecondArea();
+ void testThirdArea();
+
+ private:
+ ShapeRecognMeshBuilder *srMesh = 0;
+ const Areas *areas;
+ };
+};
+
+#endif // __CYLINDERTEST_HXX__
--- /dev/null
+#include "MathOpsTest.hxx"
+#include "MathOps.hxx"
+
+using namespace MEDCoupling;
+
+void MathOpsTest::testLstsq()
+{
+ std::vector<double> a = {
+ 0.69473263, 0.83318004, 0.60822673, 0.59243878, 0.82872553,
+ 0.84048546, 0.95698819, 0.02171218, 0.27683381, 0.20628928,
+ 0.80895323, 0.4207767, 0.37468575, 0.86258204, 0.42571846};
+ std::vector<double> b = {
+ 0.91167508, 0.95878824, 0.7234827, 0.51753917, 0.18603306};
+ std::vector<double> x = MathOps::lstsq(a, b);
+ std::array<double, 3> xRef = {0.35719095, 0.5134345, 0.26039343};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(xRef[i], x[i], 1E-6);
+}
+
+void MathOpsTest::testLstsq2()
+{
+ std::vector<double> a = {
+ 0.4564562, 0.3517006,
+ 0.28928215, 0.72309086,
+ 0.05944836, 0.56024464};
+ std::vector<double> b = {0.98902712, 0.46791812};
+ std::vector<double> x = MathOps::lstsq(a, b);
+ std::array<double, 3> xRef = {2.10752524, 0.2636243, -0.82807416};
+ for (size_t i = 0; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(xRef[i], x[i], 1E-6);
+}
+
+void MathOpsTest::testLstsqBig()
+{
+ std::vector<double> a(5000000 * 3, 1.0);
+ std::vector<double> b(5000000, 1.0);
+ std::vector<double> x = MathOps::lstsq(a, b);
+ std::array<double, 3> xRef = {1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0};
+ for (size_t i = 2; i < 3; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(xRef[i], x[i], 1E-6);
+}
+
+void MathOpsTest::testComputeCov()
+{
+ std::vector<double> coordinates{
+ 1, 2, -3, -8, 6, -3, 9, 2, -1, -10, -11, -120};
+ std::array<double, 9>
+ covMatrix = MathOps::computeCov(coordinates);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(76.666666, covMatrix[0], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(26.666666, covMatrix[1], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(319.333333, covMatrix[2], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(26.666666, covMatrix[3], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(54.916666, covMatrix[4], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(420.75, covMatrix[5], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(319.333333, covMatrix[6], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(420.75, covMatrix[7], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(3462.25, covMatrix[8], 1E-6);
+}
+
+void MathOpsTest::testComputePCAFirstAxis()
+{
+ std::vector<double> coordinates{
+ 1, 2, -3, -8, 6, -3, 9, 2, -1, -10, -11, -120};
+ std::array<double, 3>
+ axis = MathOps::computePCAFirstAxis(coordinates);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09198798, axis[0], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.11994164, axis[1], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.9885101, axis[2], 1E-6);
+}
+
+void MathOpsTest::testComputeAngles()
+{
+ std::vector<double> directions{
+ 1, 2, -3, -8, 6, -3, 9, 2, -1, -10, -11, -120};
+ std::array<double, 3> axis{1, 2, 3};
+ std::vector<double> angles = MathOps::computeAngles(
+ directions, axis);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.86054803, angles[0], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.69914333, angles[1], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.27845478, angles[2], 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(2.61880376, angles[3], 1E-6);
+}
+
+void MathOpsTest::testComputeBaseFromNormal()
+{
+ std::array<double, 3> normal = {1.0, 2.0, 3.0};
+ std::array<double, 6> base = MathOps::computeBaseFromNormal(normal);
+ std::array<double, 6> baseRef = {
+ -0.53452248, 0.77454192, -0.33818712,
+ -0.80178373, -0.33818712, 0.49271932};
+ for (size_t i = 0; i < 6; ++i)
+ {
+ std::ostringstream message;
+ message << "Mismatch at index " << i;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(message.str().c_str(), baseRef[i], base[i], 1E-6);
+ }
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __MATHOPSTEST_HXX__
+#define __MATHOPSTEST_HXX__
+
+#include <cppunit/TestFixture.h>
+#include <cppunit/extensions/HelperMacros.h>
+
+namespace MEDCoupling
+{
+ class MathOpsTest : public CppUnit::TestFixture
+ {
+ CPPUNIT_TEST_SUITE(MathOpsTest);
+ CPPUNIT_TEST(testLstsq);
+ CPPUNIT_TEST(testLstsq2);
+ CPPUNIT_TEST(testLstsqBig);
+ CPPUNIT_TEST(testComputeCov);
+ CPPUNIT_TEST(testComputePCAFirstAxis);
+ CPPUNIT_TEST(testComputeAngles);
+ CPPUNIT_TEST(testComputeBaseFromNormal);
+ CPPUNIT_TEST_SUITE_END();
+
+ public:
+ static void testLstsq();
+ static void testLstsq2();
+ static void testLstsqBig();
+ static void testComputeCov();
+ static void testComputePCAFirstAxis();
+ static void testComputeAngles();
+ static void testComputeBaseFromNormal();
+ };
+};
+
+#endif // __MATHOPSTEST_HXX__
--- /dev/null
+#include "PlaneTest.hxx"
+
+#include "ShapeRecognMeshBuilder.hxx"
+#include "Areas.hxx"
+#include "MathOps.hxx"
+#include "TestInterpKernelUtils.hxx" // getResourceFile()
+
+using namespace MEDCoupling;
+
+void PlaneTest::setUp()
+{
+ std::string file = INTERP_TEST::getResourceFile("ShapeRecognPlane.med", 3);
+ srMesh = new ShapeRecognMeshBuilder(file);
+ srMesh->recognize();
+ areas = srMesh->getAreas();
+}
+
+void PlaneTest::tearDown()
+{
+ if (srMesh != 0)
+ delete srMesh;
+ areas = 0;
+}
+
+void PlaneTest::testArea()
+{
+ CPPUNIT_ASSERT_EQUAL(36, (int)areas->getNumberOfNodes(0));
+ CPPUNIT_ASSERT_EQUAL(1, (int)areas->getNumberOfAreas());
+ // Normal
+ std::array<double, 3> normal = areas->getNormal(0);
+ std::array<double, 3> normalRef = {0.781525, 0.310606, -0.541056};
+ std::array<double, 3> affinePoint = areas->getAffinePoint(0);
+ double proportion0 = normal[0] / normalRef[0];
+ double proportion1 = normal[1] / normalRef[1];
+ double proportion2 = normal[2] / normalRef[2];
+ double proportion3 = MathOps::dot(normal, affinePoint) / MathOps::dot(normalRef, affinePoint);
+ // Check proportions between the normal vectors of the two planes
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(proportion0, proportion1, 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(proportion1, proportion2, 1E-6);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(proportion2, proportion3, 1E-6);
+ // Check the angle
+ double angle = MathOps::computeAngle(normal, normalRef);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, angle, 1E-6);
+}
\ No newline at end of file
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __PLANETEST_HXX__
+#define __PLANETEST_HXX__
+
+#include <cppunit/TestFixture.h>
+#include <cppunit/extensions/HelperMacros.h>
+
+namespace MEDCoupling
+{
+ class ShapeRecognMeshBuilder;
+ class Areas;
+
+ class PlaneTest : public CppUnit::TestFixture
+ {
+ CPPUNIT_TEST_SUITE(PlaneTest);
+ CPPUNIT_TEST(testArea);
+ CPPUNIT_TEST_SUITE_END();
+
+ public:
+ void setUp() override;
+ void tearDown() override;
+
+ void testArea();
+
+ private:
+ ShapeRecognMeshBuilder *srMesh = 0;
+ const Areas *areas;
+ };
+};
+
+#endif // __PLANETEST_HXX__
--- /dev/null
+#include "SphereTest.hxx"
+
+#include "ShapeRecognMeshBuilder.hxx"
+#include "Areas.hxx"
+#include "MathOps.hxx"
+#include "TestInterpKernelUtils.hxx" // getResourceFile()
+
+using namespace MEDCoupling;
+
+void SphereTest::setUp()
+{
+ std::string file = INTERP_TEST::getResourceFile("ShapeRecognSphere.med", 3);
+ srMesh = new ShapeRecognMeshBuilder(file);
+ srMesh->recognize();
+ areas = srMesh->getAreas();
+}
+
+void SphereTest::tearDown()
+{
+ if (srMesh != 0)
+ delete srMesh;
+ areas = 0;
+}
+
+void SphereTest::testArea()
+{
+ CPPUNIT_ASSERT_EQUAL(1, (int)areas->getNumberOfAreas());
+ // 8 double nodes so 147 - 6 nodes
+ CPPUNIT_ASSERT_EQUAL(141, (int)areas->getNumberOfNodes(0));
+ CPPUNIT_ASSERT_EQUAL(PrimitiveType::Sphere, areas->getPrimitiveType(0));
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, areas->getRadius(0), 1E-2);
+ std::array<double, 3> centerRef = {5.3, -6.7, -9.02};
+ std::array<double, 3> center = areas->getCenter(0);
+ for (size_t j = 0; j < 3; ++j)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(centerRef[j], center[j], 1E-2);
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __SPHERETEST_HXX__
+#define __SPHERETEST_HXX__
+
+#include <cppunit/TestFixture.h>
+#include <cppunit/extensions/HelperMacros.h>
+
+namespace MEDCoupling
+{
+ class ShapeRecognMeshBuilder;
+ class Areas;
+
+ class SphereTest : public CppUnit::TestFixture
+ {
+ CPPUNIT_TEST_SUITE(SphereTest);
+ CPPUNIT_TEST(testArea);
+ CPPUNIT_TEST_SUITE_END();
+
+ public:
+ void setUp() override;
+ void tearDown() override;
+
+ void testArea();
+
+ private:
+ ShapeRecognMeshBuilder *srMesh = 0;
+ const Areas *areas;
+ };
+};
+
+#endif // __SPHERETEST_HXX__
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#include "MathOpsTest.hxx"
+#include "PlaneTest.hxx"
+#include "CylinderTest.hxx"
+#include "ConeTest.hxx"
+#include "SphereTest.hxx"
+#include "TorusTest.hxx"
+
+CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::MathOpsTest);
+CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::PlaneTest);
+CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::CylinderTest);
+CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::ConeTest);
+CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::SphereTest);
+CPPUNIT_TEST_SUITE_REGISTRATION(MEDCoupling::TorusTest);
+
+#include "BasicMainTest.hxx"
--- /dev/null
+#include "TorusTest.hxx"
+
+#include "ShapeRecognMeshBuilder.hxx"
+#include "Areas.hxx"
+#include "MathOps.hxx"
+#include "TestInterpKernelUtils.hxx" // getResourceFile()
+
+using namespace MEDCoupling;
+
+void TorusTest::setUp()
+{
+ std::string file = INTERP_TEST::getResourceFile("ShapeRecognTorus.med", 3);
+ srMesh = new ShapeRecognMeshBuilder(file);
+ srMesh->recognize();
+ areas = srMesh->getAreas();
+}
+
+void TorusTest::tearDown()
+{
+ if (srMesh != 0)
+ delete srMesh;
+ areas = 0;
+}
+
+void TorusTest::testArea()
+{
+ CPPUNIT_ASSERT_EQUAL(275, (int)srMesh->getNodes()->getNbNodes());
+ CPPUNIT_ASSERT_EQUAL(1, (int)areas->getNumberOfAreas());
+ CPPUNIT_ASSERT_EQUAL(PrimitiveType::Torus, areas->getPrimitiveType(0));
+ // Some nodes are unknown
+ CPPUNIT_ASSERT_EQUAL(272, (int)areas->getNumberOfNodes(0));
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.843297, areas->getMinorRadius(0), 1E-2);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.156428, areas->getRadius(0), 1E-1);
+ std::array<double, 3> centerRef = {7.687022, -3.726887, -9.02};
+ std::array<double, 3> center = areas->getCenter(0);
+ for (size_t j = 0; j < 3; ++j)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(centerRef[j], center[j], 1E-2);
+}
--- /dev/null
+// Copyright (C) 2007-2024 CEA, EDF
+//
+// 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 2.1 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, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __TORUSTEST_HXX__
+#define __TORUSTEST_HXX__
+
+#include <cppunit/TestFixture.h>
+#include <cppunit/extensions/HelperMacros.h>
+
+namespace MEDCoupling
+{
+ class ShapeRecognMeshBuilder;
+ class Areas;
+
+ class TorusTest : public CppUnit::TestFixture
+ {
+ CPPUNIT_TEST_SUITE(TorusTest);
+ CPPUNIT_TEST(testArea);
+ CPPUNIT_TEST_SUITE_END();
+
+ public:
+ void setUp() override;
+ void tearDown() override;
+
+ void testArea();
+
+ private:
+ ShapeRecognMeshBuilder *srMesh = 0;
+ const Areas *areas;
+ };
+};
+
+#endif // __TORUSTEST_HXX__
