Skip to content

Commit

Permalink
levelset: intersect interface with continuous surface
Browse files Browse the repository at this point in the history
  • Loading branch information
Konstantinos committed Jul 24, 2024
1 parent e7bde60 commit 380f09b
Show file tree
Hide file tree
Showing 6 changed files with 286 additions and 14 deletions.
79 changes: 79 additions & 0 deletions src/levelset/levelSetObject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1402,6 +1402,47 @@ 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
* "invert" argument is false.
*
* @param[in] id is the interface index
* @param[in] invert when false 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 invert,
std::array<std::array<double, 3>, 2> *intersection,
std::vector<std::array<double, 3>> *polygon) const
{
return _isInterfaceIntersected(id, invert, 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<std::array<double, 3>, 2> intersection;
std::vector<std::array<double, 3>> polygon;
return isInterfaceIntersected(id, false, &intersection, &polygon);
}

/*!
* Check if the specified cell intersects the zero-levelset iso-surface.
*
Expand Down Expand Up @@ -1511,6 +1552,44 @@ 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
* "invert" argument is false.
*
* @param[in] id is the interface index
* @param[in] invert when false the part of the interface occupying the possitive levelset area is returned
* @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 invert,
std::array<std::array<double, 3>, 2> *intersection,
std::vector<std::array<double, 3>> *polygon) const
{
const Interface &interface = m_kernel->getMesh()->getInterface(id);
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 (!invert) {
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
Expand Down
7 changes: 5 additions & 2 deletions src/levelset/levelSetObject.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,8 +90,10 @@ friend class LevelSetProxyObject;
virtual bool isCellInNarrowBand(long id) const;
virtual bool isInNarrowBand(const std::array<double,3> &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 invert, std::array<std::array<double, 3>, 2> *intersection, std::vector<std::array<double, 3>> *polygon) const;
LevelSetIntersectionStatus isInterfaceIntersected(long id) const;

virtual short evalCellSign(long id) const;
virtual double evalCellValue(long id, bool signedLevelSet) const;
Expand Down Expand Up @@ -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, std::array<std::array<double, 3>, 2> *intersection, std::vector<std::array<double, 3>> *polygon) const;

virtual short _evalCellSign(long id) const = 0;
virtual double _evalCellValue(long id, bool signedLevelSet) const = 0;
Expand Down
120 changes: 120 additions & 0 deletions src/levelset/levelSetSegmentationObject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1817,6 +1817,94 @@ 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
* "invert" argument is false.
*
* The intersection is computed by considering the zero-levelset as a curved surface.
*
* @param[in] id is the interface index
* @param[in] invert when false the part of the interface occupying the positive levelset area is returned
* @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 invert,
std::array<std::array<double, 3>, 2> *intersection,
std::vector<std::array<double, 3>> *polygon) const
{

// Get the tolerance used for distance comparisons
double tolerance = getKernel()->getDistanceTolerance();

// Evaluate projection of interface centroid on surface
const Interface &interface = m_kernel->getMesh()->getInterface(id);
std::array<double,3> interfaceCentroid = m_kernel->getMesh()->evalElementCentroid(interface);

std::array<double,3> projectionPoint;
std::array<double,3> projectionNormal;
evalProjection(interfaceCentroid, true, &projectionPoint, &projectionNormal);
if (!invert) {
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<double, 3> intersectionCentroid = interfaceCentroid;
std::array<double, 3> 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, true, &projectionPoint, &projectionNormal);
if (!invert) {
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.
*
Expand Down Expand Up @@ -3187,6 +3275,38 @@ void LevelSetSegmentationObject::_evalProjection(const std::array<double,3> &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
* "invert" argument is false.
*
* 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] invert when false the part of the interface occupying the positive levelset area is returned
* @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 invert,
std::array<std::array<double, 3>, 2> *intersection,
std::vector<std::array<double, 3>> *polygon) const
{
// Early return if low order smoothing is chosen
if (m_surfaceInfo->getSurfaceSmoothing() == LevelSetSurfaceSmoothing::LOW_ORDER) {
return LevelSetObject::_isInterfaceIntersected(id, invert, intersection, polygon);
}
return LevelSetSegmentationBaseObject::_isInterfaceIntersected(id, invert, intersection, polygon);
}

/*!
* Get the smallest characteristic size within the triangulation
* This function is only provided for guarantee backwards compatibility with older versions.
Expand Down
3 changes: 3 additions & 0 deletions src/levelset/levelSetSegmentationObject.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 invert, std::array<std::array<double, 3>, 2> *intersection, std::vector<std::array<double, 3>> *polygon) const override;

virtual const SurfUnstructured & _evalCellSurface(long id) const = 0;
virtual int _evalCellPart(long id) const;
Expand Down Expand Up @@ -232,6 +233,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 _isInterfaceIntersected(long id, bool invert, std::array<std::array<double, 3>, 2> *intersection, std::vector<std::array<double, 3>> *polygon) const override;

private:
std::unique_ptr<LevelSetSegmentationSurfaceInfo> m_surfaceInfo;

Expand Down
48 changes: 40 additions & 8 deletions test/integration_tests/levelset/test_levelset_00001.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
{
Expand Down Expand Up @@ -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<std::chrono::system_clock> start, end;
Expand Down Expand Up @@ -246,19 +248,49 @@ int subtest_002()

levelSetSegmentation->evalProjection(point, true, &projectionPoint, &projectionNormal);

end = std::chrono::system_clock::now();
elapsed_seconds = std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count();
bitpit::log::cout() << "elapsed time: " << elapsed_seconds << " ms" << std::endl;

std::array<double, 3> projectionPointTarget = {0.0021934049109738, -0.008062212041176, 0.0};
std::array<double, 3> projectionNormalTarget = {-0.87574759289519, -0.48276925496379, 0.0};

double pointDeviation = norm2(projectionPoint - projectionPointTarget);
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<std::array<double, 3>, 2> intersection;
std::vector<std::array<double, 3>> polygon;
bool intersectionFound = (levelSetSegmentation->isInterfaceIntersected(intrId, false, &intersection, &polygon) == bitpit::LevelSetIntersectionStatus::TRUE);

std::array<double, 3> intersectionMean = 0.5 * (intersection[0] + intersection[1]);

int size = polygon.size();
std::unique_ptr<long[]> 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<double, 3> polygonMean = element.evalCentroid(polygon.data());

std::array<double, 3> intersectionMeanTarget = {0.234654275, -0.059034230965178, 0.0};
std::array<double, 3> 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<std::chrono::milliseconds>(end-start).count();
bitpit::log::cout() << "elapsed time: " << elapsed_seconds << " ms" << std::endl;

return !(projectionFound && intersectionFound);
}

/*!
Expand Down
Loading

0 comments on commit 380f09b

Please sign in to comment.