diff --git a/src/levelset/levelSetObject.cpp b/src/levelset/levelSetObject.cpp index b8c9953300..6ecae15602 100644 --- a/src/levelset/levelSetObject.cpp +++ b/src/levelset/levelSetObject.cpp @@ -1402,6 +1402,48 @@ LevelSetIntersectionStatus LevelSetObject::isCellIntersected(long id, LevelSetIn return _isCellIntersected(id, distance, mode); } +/*! + * Check if the specified interface intersects the zero-levelset iso-surface. + * + * The iso-surface is considered planar. Calculate the intersection + * as a segment defined by two points and the part of the interface belonging to the subspace + * of negative level-set. If the interface is not intersected, an empty polygon is returned. + * The interface belonging to the possitive level-set subspace is returned if the + * "positivePart" argument is true. + * + * @param[in] id is the interface index + * @param[in] positivePart when true the part of the interface occupying the possitive levelset area is returned + * @param[in] signedLevelSet controls if signed levelset function will be used + * @param[out] intersection If intersects, return the cutting segment. In case of a liner interface + * the intersection degenerates to a point which coincides with the begining and the ending + * of the segment + * @param[out] polygon is the polygon constructed by the interface intersection placed on the + * opposite subspace of the one pointed by the plane's normal + * @return the level set status + */ +LevelSetIntersectionStatus LevelSetObject::isInterfaceIntersected(long id, bool positivePart, + std::array, 2> *intersection, + std::vector> *polygon) const +{ + double distanceTolerance = m_kernel->getDistanceTolerance(); + return _isInterfaceIntersected(id, positivePart, true, distanceTolerance, intersection, polygon); +} + +/*! + * Check if the specified interface intersects the zero-levelset iso-surface. + * + * The iso-surface is considered planar. + * + * @param[in] id is the interface index + * @return the level set status + */ +LevelSetIntersectionStatus LevelSetObject::isInterfaceIntersected(long id) const +{ + std::array, 2> intersection; + std::vector> polygon; + return isInterfaceIntersected(id, false, &intersection, &polygon); +} + /*! * Check if the specified cell intersects the zero-levelset iso-surface. * @@ -1511,6 +1553,50 @@ LevelSetIntersectionStatus LevelSetObject::_isCellIntersected(long id, double di } +/*! + * Check if the specified interface intersects the zero-levelset iso-surface. + * + * The iso-surface is considered planar. Calculate the intersection + * as a segment defined by two points and the part of the interface belonging to the subspace + * of negative level-set. If the interface is not intersected, an empty polygon is returned. + * The interface belonging to the possitive level-set subspace is returned if the + * "positivePart" argument is true. + * + * @param[in] id is the interface index + * @param[in] positivePart when true the part of the interface occupying the possitive levelset area is returned + * @param[in] signedLevelSet controls if signed levelset function will be used + * @param[in] tolerance is the tolerance used for distance comparisons + * @param[out] intersection If intersects, return the cutting segment. In case of a liner interface + * the intersection degenerates to a point which coincides with the begining and the ending + * of the segment + * @param[out] polygon is the polygon constructed by the interface intersection placed on the + * opposite subspace of the one pointed by the plane's normal + * @return the level set status + */ +LevelSetIntersectionStatus LevelSetObject::_isInterfaceIntersected(long id, bool positivePart, + bool signedLevelSet, double tolerance, + std::array, 2> *intersection, + std::vector> *polygon) const +{ + BITPIT_UNUSED(signedLevelSet); + BITPIT_UNUSED(tolerance); + + const Interface &interface = m_kernel->getMesh()->getInterface(id); + std::array centroid = m_kernel->getMesh()->evalElementCentroid(interface); + + std::array root = evalProjectionPoint(centroid); + std::array normal = evalGradient(centroid, true); + if (positivePart) { + normal *= -1.; + } + + if( m_kernel->getMesh()->intersectInterfacePlane(id, root, normal, intersection, polygon) ){ + return LevelSetIntersectionStatus::TRUE; + } else { + return LevelSetIntersectionStatus::FALSE; + } +} + /*! * Evaluate levelset sign at the specified cell. * @param[in] id cell id diff --git a/src/levelset/levelSetObject.hpp b/src/levelset/levelSetObject.hpp index f28a555fc5..bc0f543c56 100644 --- a/src/levelset/levelSetObject.hpp +++ b/src/levelset/levelSetObject.hpp @@ -90,8 +90,10 @@ friend class LevelSetProxyObject; virtual bool isCellInNarrowBand(long id) const; virtual bool isInNarrowBand(const std::array &point) const; - BITPIT_DEPRECATED(LevelSetIntersectionStatus intersectSurface(long, LevelSetIntersectionMode=LevelSetIntersectionMode::FAST_FUZZY) const); - LevelSetIntersectionStatus isCellIntersected(long, LevelSetIntersectionMode=LevelSetIntersectionMode::FAST_FUZZY) const; + BITPIT_DEPRECATED(LevelSetIntersectionStatus intersectSurface(long id, LevelSetIntersectionMode=LevelSetIntersectionMode::FAST_FUZZY) const); + LevelSetIntersectionStatus isCellIntersected(long id, LevelSetIntersectionMode=LevelSetIntersectionMode::FAST_FUZZY) const; + LevelSetIntersectionStatus isInterfaceIntersected(long id, bool positivePart, std::array, 2> *intersection, std::vector> *polygon) const; + LevelSetIntersectionStatus isInterfaceIntersected(long id) const; virtual short evalCellSign(long id) const; virtual double evalCellValue(long id, bool signedLevelSet) const; @@ -181,6 +183,7 @@ friend class LevelSetProxyObject; virtual void restore(std::istream &); virtual LevelSetIntersectionStatus _isCellIntersected(long, double distance, LevelSetIntersectionMode=LevelSetIntersectionMode::FAST_FUZZY) const; + virtual LevelSetIntersectionStatus _isInterfaceIntersected(long id, bool positivePart, bool signedLevelSet, double tolerance, std::array, 2> *intersection, std::vector> *polygon) const; virtual short _evalCellSign(long id) const = 0; virtual double _evalCellValue(long id, bool signedLevelSet) const = 0; diff --git a/src/levelset/levelSetSegmentationObject.cpp b/src/levelset/levelSetSegmentationObject.cpp index 8d7057b078..dffa3df48f 100644 --- a/src/levelset/levelSetSegmentationObject.cpp +++ b/src/levelset/levelSetSegmentationObject.cpp @@ -1817,6 +1817,93 @@ LevelSetIntersectionStatus LevelSetSegmentationBaseObject::_isCellIntersected(lo return LevelSetObject::_isCellIntersected(id, distance, mode); } +/*! + * Check if the specified interface intersects the zero-levelset iso-surface. + * + * The iso-surface is considered planar. Calculate the intersection + * as a segment defined by two points and the part of the interface belonging to the subspace + * of negative level-set. If the interface is not intersected, an empty polygon is returned. + * The interface belonging to the possitive level-set subspace is returned if the + * "positivePart" argument is true. + * + * The intersection is computed by considering the zero-levelset as a curved surface. + * + * @param[in] id is the interface index + * @param[in] positivePart when true the part of the interface occupying the positive levelset area is returned + * @param[in] signedLevelSet controls if signed levelset function will be used + * @param[in] tolerance is the tolerance used for distance comparisons + * @param[out] intersection If intersects, return the cutting segment. In case of a liner interface + * the intersection degenerates to a point which coincides with the begining and the ending + * of the segment + * @param[out] polygon is the polygon constructed by the interface intersection placed on the + * opposite subspace of the one pointed by the plane's normal + * @return the level set status + */ +LevelSetIntersectionStatus LevelSetSegmentationBaseObject::_isInterfaceIntersected(long id, bool positivePart, + bool signedLevelSet, double tolerance, + std::array, 2> *intersection, + std::vector> *polygon) const +{ + // Evaluate projection of interface centroid on surface + const Interface &interface = m_kernel->getMesh()->getInterface(id); + std::array interfaceCentroid = m_kernel->getMesh()->evalElementCentroid(interface); + + std::array projectionPoint; + std::array projectionNormal; + evalProjection(interfaceCentroid, signedLevelSet, &projectionPoint, &projectionNormal); + if (positivePart) { + projectionNormal *= -1.; + } + + // Early return if projection point is not the closest point of the + // plane (defined by the projection point and projection normal) to + // the interface centroid + long support = evalSupport(interfaceCentroid); + const SurfUnstructured &surface = evalSurface(interfaceCentroid); + LevelSetSegmentationSurfaceInfo::SegmentConstIterator segmentItr = surface.getCellConstIterator(support); + + const Cell &segment = *segmentItr; + ElementType segmentType = segment.getType(); + if ((segmentType == ElementType::VERTEX) + || (m_kernel->getMesh()->getDimension() == 3 && segmentType == ElementType::LINE)) { + return LevelSetIntersectionStatus::FALSE; + } + + // Detect the intersection by imposing a Newton algorithm + std::array intersectionCentroid = interfaceCentroid; + std::array intersectionCentroidNew; + + double residual = tolerance + 1.0; + int iter = 0; + int iterMax = 5; + bool isIntersected = true; + while (residual > tolerance && iter < iterMax && isIntersected) { + ++ iter; + + // Eval intersection between arbitrary 2D surface and planar interface + isIntersected = m_kernel->getMesh()->intersectInterfacePlane(id, projectionPoint, projectionNormal, intersection, polygon); + + intersectionCentroidNew = 0.5 * ((*intersection)[0] + (*intersection)[1]); + + evalProjection(intersectionCentroidNew, signedLevelSet, &projectionPoint, &projectionNormal); + if (positivePart) { + projectionNormal *= -1.; + } + + residual = norm2(intersectionCentroid - intersectionCentroidNew); + intersectionCentroid = intersectionCentroidNew; + } + + // Eval intersection corresponding to the latest projection point and normal + if (isIntersected) { + isIntersected = m_kernel->getMesh()->intersectInterfacePlane(id, projectionPoint, projectionNormal, intersection, polygon); + if (isIntersected) { + return LevelSetIntersectionStatus::TRUE; + } + } + return LevelSetIntersectionStatus::FALSE; +} + /*! * Evaluate the normal of the surface at the segment closest to the specified cell. * @@ -3163,6 +3250,41 @@ void LevelSetSegmentationObject::_evalProjection(const std::array &poi } } +/*! + * Check if the specified interface intersects the zero-levelset iso-surface. + * + * The iso-surface is considered planar. Calculate the intersection + * as a segment defined by two points and the part of the interface belonging to the subspace + * of negative level-set. If the interface is not intersected, an empty polygon is returned. + * The interface belonging to the possitive level-set subspace is returned if the + * "positivePart" argument is true. + * + * Depending on the smoothing order, the intersection is computed by considering the + * zero-levelset as a plane or a curved surface. + * + * @param[in] id is the interface index + * @param[in] positivePart when true the part of the interface occupying the positive levelset area is returned + * @param[in] signedLevelSet controls if signed levelset function will be used + * @param[in] tolerance is the tolerance used for distance comparisons + * @param[out] intersection If intersects, return the cutting segment. In case of a liner interface + * the intersection degenerates to a point which coincides with the begining and the ending + * of the segment + * @param[out] polygon is the polygon constructed by the interface intersection placed on the + * opposite subspace of the one pointed by the plane's normal + * @return the level set status + */ +LevelSetIntersectionStatus LevelSetSegmentationObject::_isInterfaceIntersected(long id, bool positivePart, + bool signedLevelSet, double tolerance, + std::array, 2> *intersection, + std::vector> *polygon) const +{ + // Early return if low order smoothing is chosen + if (m_surfaceInfo->getSurfaceSmoothing() == LevelSetSurfaceSmoothing::LOW_ORDER) { + return LevelSetObject::_isInterfaceIntersected(id, positivePart, signedLevelSet, tolerance, intersection, polygon); + } + return LevelSetSegmentationBaseObject::_isInterfaceIntersected(id, positivePart, signedLevelSet, tolerance, intersection, polygon); +} + /*! * Get the smallest characteristic size within the triangulation * This function is only provided for guarantee backwards compatibility with older versions. diff --git a/src/levelset/levelSetSegmentationObject.hpp b/src/levelset/levelSetSegmentationObject.hpp index 3414582cd0..47d10559ae 100644 --- a/src/levelset/levelSetSegmentationObject.hpp +++ b/src/levelset/levelSetSegmentationObject.hpp @@ -165,6 +165,7 @@ class LevelSetSegmentationBaseObject : public LevelSetObject { using LevelSetObject::fillFieldCellCache; LevelSetIntersectionStatus _isCellIntersected(long, double distance, LevelSetIntersectionMode=LevelSetIntersectionMode::FAST_FUZZY) const override; + virtual LevelSetIntersectionStatus _isInterfaceIntersected(long id, bool positivePart, bool signedLevelSet, double tolerance, std::array, 2> *intersection, std::vector> *polygon) const override; virtual const SurfUnstructured & _evalCellSurface(long id) const = 0; virtual int _evalCellPart(long id) const; @@ -232,6 +233,8 @@ class LevelSetSegmentationObject : public LevelSetSegmentationBaseObject { long _evalSupport(const std::array &point, double searchRadius) const override; void _evalProjection(const std::array &point, bool signedLevelSet, std::array *projectionPoint, std::array *projectionNormal) const override; + LevelSetIntersectionStatus _isInterfaceIntersected(long id, bool positivePart, bool signedLevelSet, double tolerance, std::array, 2> *intersection, std::vector> *polygon) const override; + private: std::unique_ptr m_surfaceInfo; diff --git a/test/integration_tests/levelset/test_levelset_00001.cpp b/test/integration_tests/levelset/test_levelset_00001.cpp index ec407bd143..951f0d067d 100644 --- a/test/integration_tests/levelset/test_levelset_00001.cpp +++ b/test/integration_tests/levelset/test_levelset_00001.cpp @@ -178,7 +178,8 @@ int subtest_001() /*! * Subtest 002 * -* Testing projection on surface feature of a 2D levelset on a Cartesian mesh in default memory mode. +* Testing projection and interface intersection on surface feature of a 2D levelset on +* a Cartesian mesh in default memory mode. */ int subtest_002() { @@ -214,9 +215,10 @@ int subtest_002() delta = meshMax -meshMin ; bitpit::VolCartesian mesh( 1, dimensions, meshMin, delta, nc); - mesh.update() ; + mesh.switchMemoryMode(bitpit::VolCartesian::MEMORY_NORMAL); mesh.initializeAdjacencies() ; mesh.initializeInterfaces() ; + mesh.update() ; // Compute level set in narrow band std::chrono::time_point start, end; @@ -246,10 +248,6 @@ int subtest_002() levelSetSegmentation->evalProjection(point, true, &projectionPoint, &projectionNormal); - end = std::chrono::system_clock::now(); - elapsed_seconds = std::chrono::duration_cast(end-start).count(); - bitpit::log::cout() << "elapsed time: " << elapsed_seconds << " ms" << std::endl; - std::array projectionPointTarget = {0.0021934049109738, -0.008062212041176, 0.0}; std::array projectionNormalTarget = {-0.87574759289519, -0.48276925496379, 0.0}; @@ -257,8 +255,42 @@ int subtest_002() double normalDeviation = norm2(projectionNormal - projectionNormalTarget); double distanceTolerance = mesh.getTol(); - return !(bitpit::utils::DoubleFloatingEqual()(pointDeviation, 0., distanceTolerance, distanceTolerance) - ||bitpit::utils::DoubleFloatingEqual()(normalDeviation, 0., distanceTolerance, distanceTolerance)); + bool projectionFound = (bitpit::utils::DoubleFloatingEqual()(pointDeviation, 0., distanceTolerance, distanceTolerance) + && bitpit::utils::DoubleFloatingEqual()(normalDeviation, 0., distanceTolerance, distanceTolerance)); + + // Compute intersected interface + int intrId = 744; + + std::array, 2> intersection; + std::vector> polygon; + bool intersectionFound = (levelSetSegmentation->isInterfaceIntersected(intrId, true, &intersection, &polygon) == bitpit::LevelSetIntersectionStatus::TRUE); + + std::array intersectionMean = 0.5 * (intersection[0] + intersection[1]); + + int size = polygon.size(); + std::unique_ptr connectivity(new long[size + 1]); + connectivity[0] = size; + for (int i = 0; i < size; ++i) { + connectivity[i + 1] = i; + } + + bitpit::Element element(0, bitpit::ElementType::POLYGON, std::move(connectivity)); + std::array polygonMean = element.evalCentroid(polygon.data()); + + std::array intersectionMeanTarget = {0.234654275, -0.059034230965178, 0.0}; + std::array polygonMeanTarget = {0.234654275, -0.059897734876829, 0.0}; + + double intersectionDeviation = norm2(intersectionMean - intersectionMeanTarget); + double polygonDeviation = norm2(polygonMean - polygonMeanTarget); + + intersectionFound = intersectionFound && (bitpit::utils::DoubleFloatingEqual()(intersectionDeviation, 0., distanceTolerance, distanceTolerance) + && bitpit::utils::DoubleFloatingEqual()(polygonDeviation, 0., distanceTolerance, distanceTolerance)); + + end = std::chrono::system_clock::now(); + elapsed_seconds = std::chrono::duration_cast(end-start).count(); + bitpit::log::cout() << "elapsed time: " << elapsed_seconds << " ms" << std::endl; + + return !(projectionFound && intersectionFound); } /*! diff --git a/test/integration_tests/levelset/test_levelset_00002.cpp b/test/integration_tests/levelset/test_levelset_00002.cpp index d322c62705..fc2d489cda 100644 --- a/test/integration_tests/levelset/test_levelset_00002.cpp +++ b/test/integration_tests/levelset/test_levelset_00002.cpp @@ -483,8 +483,8 @@ int subtest_004() /*! * Subtest 005 * -* Testing projection on surface feature of a 3D levelset on a Cartesian mesh in default memory mode. -* Surface consist of triangular elements. +* Testing projection and interface intersection on surface feature of a 3D levelset on a Cartesian mesh +* in default memory mode. Surface consist of triangular elements. */ int subtest_005() { @@ -515,6 +515,7 @@ int subtest_005() std::unique_ptr mesh = generateCartesianMesh(*segmentation); mesh->switchMemoryMode(bitpit::VolCartesian::MEMORY_NORMAL); mesh->initializeAdjacencies(); + mesh->initializeInterfaces() ; mesh->update(); // Initialize levelset @@ -558,8 +559,42 @@ int subtest_005() double normalDeviation = norm2(projectionNormal - projectionNormalTarget); double distanceTolerance = mesh->getTol(); - return !(bitpit::utils::DoubleFloatingEqual()(pointDeviation, 0., distanceTolerance, distanceTolerance) - ||bitpit::utils::DoubleFloatingEqual()(normalDeviation, 0., distanceTolerance, distanceTolerance)); + bool projectionFound = (bitpit::utils::DoubleFloatingEqual()(pointDeviation, 0., distanceTolerance, distanceTolerance) + && bitpit::utils::DoubleFloatingEqual()(normalDeviation, 0., distanceTolerance, distanceTolerance)); + + // Compute intersected interface + int intrId = 71736; + + std::array, 2> intersection; + std::vector> polygon; + bool intersectionFound = (levelSetSegmentation->isInterfaceIntersected(intrId, true, &intersection, &polygon) == bitpit::LevelSetIntersectionStatus::TRUE); + + std::array intersectionMean = 0.5 * (intersection[0] + intersection[1]); + + int size = polygon.size(); + std::unique_ptr connectivity(new long[size + 1]); + connectivity[0] = size; + for (int i = 0; i < size; ++i) { + connectivity[i + 1] = i; + } + + bitpit::Element element(0, bitpit::ElementType::POLYGON, std::move(connectivity)); + std::array polygonMean = element.evalCentroid(polygon.data()); + + std::array intersectionMeanTarget = {-0.0534969303756953, -0.0926593981683255, -0.99584790177713}; + std::array polygonMeanTarget = {-0.0534969303756953, -0.0926593981683255, -1.00417395088856}; + + double intersectionDeviation = norm2(intersectionMean - intersectionMeanTarget); + double polygonDeviation = norm2(polygonMean - polygonMeanTarget); + + intersectionFound = intersectionFound && (bitpit::utils::DoubleFloatingEqual()(intersectionDeviation, 0., distanceTolerance, distanceTolerance) + && bitpit::utils::DoubleFloatingEqual()(polygonDeviation, 0., distanceTolerance, distanceTolerance)); + + end = std::chrono::system_clock::now(); + elapsed_seconds = std::chrono::duration_cast(end-start).count(); + bitpit::log::cout() << "elapsed time: " << elapsed_seconds << " ms" << std::endl; + + return !(projectionFound && intersectionFound); } /*!