levelset: intersect interface with continuous surface

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<double,3> centroid = m_kernel->getMesh()->evalElementCentroid(interface);
+
+    std::array<double,3> root = evalProjectionPoint(centroid);
+    std::array<double,3> 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.
  *

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;

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<double, 3> *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<double,3> interfaceCentroid = m_kernel->getMesh()->evalElementCentroid(interface);
+
+    std::array<double,3> projectionPoint;
+    std::array<double,3> 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<double,3>  root         = projectionPoint;
+    std::array<double,3>  normal       = projectionNormal;
+    std::array<double, 3> intersection = interfaceCentroid;
+    std::array<double, 3> 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<double, 3> *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<double,3> &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.

src/levelset/levelSetSegmentationObject.hpp

@@ -146,6 +146,8 @@ class LevelSetSegmentationBaseObject : public LevelSetObject {
     std::array<double,3> evalNormal(const std::array<double,3> &point, bool signedLevelSet) const;
     void evalProjection(const std::array<double,3> &point, bool signedLevelSet, std::array<double, 3> *projectionPoint, std::array<double, 3> *projectionNormal) const;
 
+    LevelSetIntersectionStatus evalIntersectedInterfaceData(long id, int intrLocalId, bool signedLevelSet, std::array<double, 3> *centroid, double *area) const;
+
     BITPIT_DEPRECATED(int getPart(long cellId) const);
     BITPIT_DEPRECATED(std::array<double BITPIT_COMMA 3> 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<double, 3> *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<double,3> &point, double searchRadius) const override;
     void _evalProjection(const std::array<double,3> &point, bool signedLevelSet, std::array<double, 3> *projectionPoint, std::array<double, 3> *projectionNormal) const override;
 
+    LevelSetIntersectionStatus _intersectInterfaceSurface(long id, int intrLocalId, double tolerance) const override;
+
 private:
     std::unique_ptr<LevelSetSegmentationSurfaceInfo> m_surfaceInfo;