From d159be5ca988fb3ff50cd17e09066ff4376053ac Mon Sep 17 00:00:00 2001 From: Konstantinos Date: Fri, 12 Apr 2024 14:32:09 +0200 Subject: [PATCH] levelset: intersect interface with continuous surface --- src/levelset/levelSetObject.cpp | 25 +++++ src/levelset/levelSetObject.hpp | 1 + src/levelset/levelSetSegmentationObject.cpp | 112 ++++++++++++++++++++ src/levelset/levelSetSegmentationObject.hpp | 5 + 4 files changed, 143 insertions(+) diff --git a/src/levelset/levelSetObject.cpp b/src/levelset/levelSetObject.cpp index 1dc59ff149..13f8a93124 100644 --- a/src/levelset/levelSetObject.cpp +++ b/src/levelset/levelSetObject.cpp @@ -1364,6 +1364,31 @@ LevelSetIntersectionStatus LevelSetObject::intersectCellSurface(long id, LevelSe return _intersectCellSurface(id, distance, mode); } +/*! + * Check if the specified cell interface intersects the zero-levelset iso-surface. + * + * The iso-surface is considered planar. + * + * @param[in] id cell id + * @param[in] intrLocalId is the local index of the interface for the given cell + * @param[in] tolerance is the tolerance used for distance comparisons + * @return the level set status + */ +LevelSetIntersectionStatus LevelSetObject::_intersectInterfaceSurface(long id, int intrLocalId, double tolerance) const +{ + const Interface &interface = m_kernel->getMesh()->getCell(id).getInterface(intrLocalId); + std::array centroid = m_kernel->getMesh()->evalElementCentroid(interface); + + std::array root = evalProjectionPoint(centroid); + std::array normal = evalGradient(centroid, true); + + if( m_kernel->intersectInterfacePlane(id, intrLocalId, root, normal, tolerance, nullptr, nullptr, nullptr) ){ + return LevelSetIntersectionStatus::TRUE; + } else { + return LevelSetIntersectionStatus::FALSE; + } +} + /*! * Check if the specified cell intersects the zero-levelset iso-surface. * diff --git a/src/levelset/levelSetObject.hpp b/src/levelset/levelSetObject.hpp index 1f54d62565..662ebe9a54 100644 --- a/src/levelset/levelSetObject.hpp +++ b/src/levelset/levelSetObject.hpp @@ -179,6 +179,7 @@ friend class LevelSetProxyObject; virtual void dump(std::ostream &); virtual void restore(std::istream &); + virtual LevelSetIntersectionStatus _intersectInterfaceSurface(long id, int intrLocalId, double tolerance) const; virtual LevelSetIntersectionStatus _intersectCellSurface(long, double distance, LevelSetIntersectionMode=LevelSetIntersectionMode::FAST_FUZZY) const; virtual short _evalCellSign(long id) const = 0; diff --git a/src/levelset/levelSetSegmentationObject.cpp b/src/levelset/levelSetSegmentationObject.cpp index 96bc70e4ec..ba8603655f 100644 --- a/src/levelset/levelSetSegmentationObject.cpp +++ b/src/levelset/levelSetSegmentationObject.cpp @@ -1589,6 +1589,100 @@ int LevelSetSegmentationBaseObject::evalCellPart(long id) const return part; } +/*! + * Check if the specified cell interface intersects the zero-levelset iso-surface. + * If an intersection is detected, the centroid and area of the intersected interface is computed. + * + * @param[in] id cell id + * @param[in] intrLocalId is the local index of the interface for the given cell + * @param[in] tolerance is the tolerance used for distance comparisons + * @param[in] signedLevelSet controls if signed levelset function will be used + * @param[in] tolerance is the tolerance used for distance comparisons + * @param[out] centroid is the centroid the intersected face computed as the arithmetic avretage of the face vertiecs + * @param[out] area is the area of the intersected face + * @return indicator regarding intersection + */ +LevelSetIntersectionStatus LevelSetSegmentationBaseObject::_evalIntersectedInterfaceData(long id, int intrLocalId, bool signedLevelSet, double tolerance, std::array *centroid, double *area) const +{ + // Evaluate projection of interface centroid on surface + long intrGlobalId = m_kernel->getMesh()->getCell(id).getInterfaces()[intrLocalId]; + const Interface &interface = m_kernel->getMesh()->getInterface(intrGlobalId); + std::array interfaceCentroid = m_kernel->getMesh()->evalElementCentroid(interface); + + std::array projectionPoint; + std::array projectionNormal; + evalProjection(interfaceCentroid, signedLevelSet, &projectionPoint, &projectionNormal); + + // 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 = evalCellSupport(id); + const SurfUnstructured &surface = evalCellSurface(id); + LevelSetSegmentationSurfaceInfo::SegmentConstIterator segmentItr = surface.getCellConstIterator(support); + const Cell &segment = *segmentItr; + ElementType segmentType = segment.getType(); + if (segmentType == ElementType::VERTEX) { + return LevelSetIntersectionStatus::FALSE; + } else if (m_kernel->getMesh()->getDimension() == 3 && segmentType == ElementType::LINE) { + return LevelSetIntersectionStatus::FALSE; + } + + // Eval intersection between arbitrary 2D surface and planar interface + std::array root = projectionPoint; + std::array normal = projectionNormal; + std::array intersection = interfaceCentroid; + std::array intersectionNew; + double residual = tolerance + 1.0; + int iter = 0; + int iterMax = 5; + bool isIntersected = true; + while (residual > tolerance && iter < iterMax && isIntersected) { + ++ iter; + + isIntersected = m_kernel->intersectInterfacePlane(id, intrLocalId, root, normal, tolerance, nullptr, nullptr, &intersectionNew); + + evalProjection(intersectionNew, signedLevelSet, &root, &normal); + + residual = norm2(intersection - intersectionNew); + intersection = intersectionNew; + } + + // Eval centroid and area of the intersected interface + if (isIntersected) { + isIntersected = m_kernel->intersectInterfacePlane(id, intrLocalId, root, normal, tolerance, centroid, area, nullptr); + if (isIntersected) { + return LevelSetIntersectionStatus::TRUE; + } + } else { + double phi = dotProduct(interfaceCentroid - projectionPoint, projectionNormal); + if (utils::DoubleFloatingGreater()(phi, 0.0, tolerance)) { + (*area) = m_kernel->getMesh()->evalInterfaceArea(intrGlobalId); + return LevelSetIntersectionStatus::FALSE; + } + } + (*area) = 0.0; + return LevelSetIntersectionStatus::FALSE; +} + +/*! + * Check if the specified cell interface intersects the zero-levelset iso-surface. + * If an intersection is detected, the centroid and area of the intersected interface is computed. + * + * @param[in] id cell id + * @param[in] intrLocalId is the local index of the interface for the given cell + * @param[in] tolerance is the tolerance used for distance comparisons + * @param[in] signedLevelSet controls if signed levelset function will be used + * @param[out] centroid is the centroid the intersected face computed as the arithmetic avretage of the face vertiecs + * @param[out] area is the area of the intersected face + * @return indicator regarding intersection + */ +LevelSetIntersectionStatus LevelSetSegmentationBaseObject::evalIntersectedInterfaceData(long id, int intrLocalId, bool signedLevelSet, std::array *centroid, double *area) const +{ + double distanceTolerance = m_kernel->getDistanceTolerance(); + + return _evalIntersectedInterfaceData(id, intrLocalId, signedLevelSet, distanceTolerance, centroid, area); +} + /*! * Check if cell intersects the surface. * @@ -2988,6 +3082,24 @@ void LevelSetSegmentationObject::_evalProjection(const std::array &poi } } +/*! + * Check if the specified cell interface intersects the zero-levelset iso-surface. + * + * @param[in] id cell id + * @param[in] intrLocalId is the local index of the interface for the given cell + * @param[in] tolerance is the tolerance used for distance comparisons + * @return indicator regarding intersection + */ +LevelSetIntersectionStatus LevelSetSegmentationObject::_intersectInterfaceSurface(long id, int intrLocalId, double tolerance) const +{ + // Early return if low order smoothing is chosen + if (m_surfaceInfo->getSurfaceSmoothing() == LevelSetSurfaceSmoothing::LOW_ORDER) { + return LevelSetObject::_intersectInterfaceSurface(id, intrLocalId, tolerance); + } + + return _evalIntersectedInterfaceData(id, intrLocalId, false, tolerance, nullptr, nullptr); +} + /*! * 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 3389bfabb4..27fa255d4b 100644 --- a/src/levelset/levelSetSegmentationObject.hpp +++ b/src/levelset/levelSetSegmentationObject.hpp @@ -146,6 +146,8 @@ class LevelSetSegmentationBaseObject : public LevelSetObject { std::array evalNormal(const std::array &point, bool signedLevelSet) const; void evalProjection(const std::array &point, bool signedLevelSet, std::array *projectionPoint, std::array *projectionNormal) const; + LevelSetIntersectionStatus evalIntersectedInterfaceData(long id, int intrLocalId, bool signedLevelSet, std::array *centroid, double *area) const; + BITPIT_DEPRECATED(int getPart(long cellId) const); BITPIT_DEPRECATED(std::array getNormal(long cellId) const); BITPIT_DEPRECATED(long getSupport(long cellId) const); @@ -165,6 +167,7 @@ class LevelSetSegmentationBaseObject : public LevelSetObject { using LevelSetObject::fillFieldCellCache; LevelSetIntersectionStatus _intersectCellSurface(long, double distance, LevelSetIntersectionMode=LevelSetIntersectionMode::FAST_FUZZY) const override; + LevelSetIntersectionStatus _evalIntersectedInterfaceData(long id, int intrLocalId, bool signedLevelSet, double tolerance, std::array *centroid, double *area) const; virtual const SurfUnstructured & _evalCellSurface(long id) const = 0; virtual int _evalCellPart(long id) const; @@ -233,6 +236,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 _intersectInterfaceSurface(long id, int intrLocalId, double tolerance) const override; + private: std::unique_ptr m_surfaceInfo;