levelset: surface smoothing

[ksamouchos]

This pull request introduces a new tool which smooths out a discretized surface (polygon or polyhedral). The smoothed surface is continuous and is equipped with a continuous normal. Also, it passes form the discretized geometry's vertices and respects the normal vectors corresponding to each vertex. The new level set computes the projection point and normal based on the new continuous surface.

This pull request completes task IMM-770.

[andrea-iob]

Since it's still a draft, I've only made some random comments without actually testing the changes.

[andrea-iob]

The code is failing to compile on windows because it's using M_PI. There is no such constant on windows, you should use BITPIT_PI/BITPIT_PI_2/BITPIT_PI_4.

M_PI is defined only if _USE_MATH_DEFINES is set, however this macro should be set before including cmath. If an external library includes cmath without defining _USE_MATH_DEFINES, then M_PI will not be available. To avoid potential problems, we have decided to create our own pi constant.

[ksamouchos]

Branch update:

• I've corrected the surface smoothing on polygonal segments
• I've added a levelset test case to control the surface smoothing over polygonal segments
• I've made the correction proposed by Andrea
• I've rebased the branch on the current master

[ksamouchos]

The code is failing to compile on windows because it's using M_PI. There is no such constant on windows, you should use BITPIT_PI/BITPIT_PI_2/BITPIT_PI_4.

M_PI is defined only if _USE_MATH_DEFINES is set, however this macro should be set before including cmath. If an external library includes cmath without defining _USE_MATH_DEFINES, then M_PI will not be available. To avoid potential problems, we have decided to create our own pi constant.

I didn't know it. Thank you for mentioning it. Now I'm using the BITPI_PI variable.

[ksamouchos]

I changed function

void setSurface(const SurfUnstructured *surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing, bool force = false);

to allow the user to modify the surface smoothing order of object LevelSetSegmentationObject after its creation.

[ksamouchos]

Rebased to the current master

[ksamouchos]

I added a detailed documentation for functions "evalHighOrderProjectionOnLine" and "evalHighOrderProjectionOnTriangle".

[ksamouchos]

I'm answering to the question of @andrea-iob about why I removed the assert from function https://github.com/optimad/bitpit/blob/4040cdd2779366d7605d15582ef7e8f5f031e17c/src/levelset/levelSetSegmentationObject.cpp#L2754

Let's have a look at one of the new integration tests I've added which is not working without removing the assert. The test is the following: https://github.com/optimad/bitpit/blob/4040cdd2779366d7605d15582ef7e8f5f031e17c/test/integration_tests/levelset/test_levelset_00001.cpp#L183

In this case the size of the narrow band is zero. Thus, function "fillCellGeometricNarrowBandLocationCache" finds valoid cell support only for intersected cells. For the rest returns an UNKNOWN location: https://github.com/optimad/bitpit/blob/4040cdd2779366d7605d15582ef7e8f5f031e17c/src/levelset/levelSetSegmentationObject.cpp#L2760

However, when a smoothed geometry is used, there are cells like the one with id=2824 which are intersected by the discretized geometry but not by the smoothed geometry. These cells have a valid cell support, so the function continues its process by checking if the location is NARROW_BAND_INTERSECTED or NARROW_BAND_DISTANCE in line https://github.com/optimad/bitpit/blob/4040cdd2779366d7605d15582ef7e8f5f031e17c/src/levelset/levelSetSegmentationObject.cpp#L2774

Cells like the one mentioned above don't belong to any of these two categories meaning that the assert that follows will throw an error. By removing the assert, the code still functions normally, because the location of these cells is declared later in function "fillCartesianCellZoneCache". Either these cells are neighboring the narrow band and detected in line https://github.com/optimad/bitpit/blob/4040cdd2779366d7605d15582ef7e8f5f031e17c/src/levelset/levelSetSegmentationObject.cpp#L2713 or belong to the bulk and are detected in line https://github.com/optimad/bitpit/blob/4040cdd2779366d7605d15582ef7e8f5f031e17c/src/levelset/levelSetSegmentationObject.cpp#L2726 More specifically, cell 2824 belongs to the second category.

Consequently, the removal of the assert is necessary and does not harm the code.

[andrea-iob]

Consequently, the removal of the assert is necessary and does not harm the code.

Ok.

Let's avoid caching unknown cell locations (like we already do in LevelSetObject::fillCellGeometricNarrowBandLocationCache) and add a comment explaining that, in the high order case, it's fine to have cells whose cell location will be unknown:

    // Update the cell location cache
    //
    // First we need to check if the cell intersects the surface, and only if it doesn't we
    // should check if its distance is lower than the narrow band size.
    //
    // When high order smoothing is active, there may be cells whose support is within the search
    // radius, but they are not intersected and their distance is less than the narrow band size.
    // These cells are not geometrically inside the narrow band, they are neighbours of cells
    // geometrically inside the narrow band and as such it's up to the caller of this function to
    // identify their cell location.
    LevelSetCellLocation cellLocation = LevelSetCellLocation::UNKNOWN;
    if (_intersectSurface(id, cellUnsigendValue, CELL_LOCATION_INTERSECTION_MODE) == LevelSetIntersectionStatus::TRUE) {
        cellLocation = LevelSetCellLocation::NARROW_BAND_INTERSECTED;
    } else if (cellUnsigendValue <= m_narrowBandSize) {
        cellLocation = LevelSetCellLocation::NARROW_BAND_DISTANCE;
    }
    assert((getSurfaceSmoothing() == LevelSetSurfaceSmoothing::HIGH_ORDER) || (cellLocation != LevelSetCellLocation::UNKNOWN));

    if (cellLocation != LevelSetCellLocation::UNKNOWN) {
        CellCacheCollection::ValueCache<char> *locationCache = getCellCache<char>(m_cellLocationCacheId);
        locationCache->insertEntry(id, static_cast<char>(cellLocation));
    }

[andrea-iob]

If you address the last two comments the pull request can be merged.

[ksamouchos]

Consequently, the removal of the assert is necessary and does not harm the code.

Ok.

Let's avoid caching unknown cell locations (like we already do in LevelSetObject::fillCellGeometricNarrowBandLocationCache) and add a comment explaining that, in the high order case, it's fine to have cells whose cell location will be unknown:

    // Update the cell location cache
    //
    // First we need to check if the cell intersects the surface, and only if it doesn't we
    // should check if its distance is lower than the narrow band size.
    //
    // When high order smoothing is active, there may be cells whose support is within the search
    // radius, but they are not intersected and their distance is less than the narrow band size.
    // These cells are not geometrically inside the narrow band, they are neighbours of cells
    // geometrically inside the narrow band and as such it's up to the caller of this function to
    // identify their cell location.
    LevelSetCellLocation cellLocation = LevelSetCellLocation::UNKNOWN;
    if (_intersectSurface(id, cellUnsigendValue, CELL_LOCATION_INTERSECTION_MODE) == LevelSetIntersectionStatus::TRUE) {
        cellLocation = LevelSetCellLocation::NARROW_BAND_INTERSECTED;
    } else if (cellUnsigendValue <= m_narrowBandSize) {
        cellLocation = LevelSetCellLocation::NARROW_BAND_DISTANCE;
    }
    assert((getSurfaceSmoothing() == LevelSetSurfaceSmoothing::HIGH_ORDER) || (cellLocation != LevelSetCellLocation::UNKNOWN));

    if (cellLocation != LevelSetCellLocation::UNKNOWN) {
        CellCacheCollection::ValueCache<char> *locationCache = getCellCache<char>(m_cellLocationCacheId);
        locationCache->insertEntry(id, static_cast<char>(cellLocation));
    }

Done

[ksamouchos]

Should I squash the commits before merging?

[andrea-iob]

Should I squash the commits before merging?

Keep the commits separate.