From facc192d5acb93a2a725fc1ebca54c78e9549e79 Mon Sep 17 00:00:00 2001 From: Andrea Iob Date: Fri, 30 Jun 2023 13:27:27 +0200 Subject: [PATCH 1/4] levelset: fix sign propagation for empty objects --- src/levelset/levelSetSignPropagator.cpp | 121 +++++++++++++++--------- 1 file changed, 74 insertions(+), 47 deletions(-) diff --git a/src/levelset/levelSetSignPropagator.cpp b/src/levelset/levelSetSignPropagator.cpp index a912771a9a..036ff374c4 100644 --- a/src/levelset/levelSetSignPropagator.cpp +++ b/src/levelset/levelSetSignPropagator.cpp @@ -406,11 +406,7 @@ void LevelSetSignPropagator::initializePropagation(const LevelSetObjectInterface // flagged as external. const LevelSetBoundedObject *boundedObject = dynamic_cast(object); - m_nExternal = 0; if (boundedObject) { - // Get tolerance for distance comparison - double distanceTolerance = m_mesh->getTol(); - // Evaluate the bounding box of the object // // The current process may only have the portion of the object needed for @@ -424,63 +420,94 @@ void LevelSetSignPropagator::initializePropagation(const LevelSetObjectInterface boundedObject->getBoundingBox(objectBoxMin, objectBoxMax); #endif - // Check if the patch intersects the bounding box of the object - std::array patchBoxMin; - std::array patchBoxMax; - m_mesh->getBoundingBox(patchBoxMin, patchBoxMax); + // Check if the bounding box of the object is empty + bool isObjectBoxEmpty = false; + for (int i = 0; i < 3; ++i) { + if (objectBoxMax[i] < objectBoxMin[i]) { + isObjectBoxEmpty = true; + break; + } + } - bool isPatchIntersected = CGElem::intersectBoxBox(patchBoxMin, patchBoxMax, objectBoxMin, objectBoxMax, 3, distanceTolerance); + // Identify cell outside the bounding box of the object + if (!isObjectBoxEmpty) { + // Get tolerance for distance comparison + double distanceTolerance = m_mesh->getTol(); - // Detect external cells - VolumeKernel::CellConstIterator cellBegin = m_mesh->cellConstBegin(); - VolumeKernel::CellConstIterator cellEnd = m_mesh->cellConstEnd(); - for (VolumeKernel::CellConstIterator itr = cellBegin; itr != cellEnd; ++itr) { - // Cells inside the narrowband cannot be external - long cellId = itr.getId(); - if (object->isInNarrowBand(cellId)) { - continue; - } + // Check if the patch intersects the bounding box of the object + std::array patchBoxMin; + std::array patchBoxMax; + m_mesh->getBoundingBox(patchBoxMin, patchBoxMax); - // Check if the centroid is inside the bounding box - // - // Cells with the centroid inside the bounding box of the object - // cannot be external cells - if (isPatchIntersected) { - double geometricTolerance = m_mesh->getTol(); - std::array cellCentroid = object->getKernel()->computeCellCentroid(cellId); - - bool isCentroidInternal = true; - for (int i = 0; i < 3; ++i) { - if (cellCentroid[i] < objectBoxMin[i] - geometricTolerance || cellCentroid[i] > objectBoxMin[i] + geometricTolerance) { - isCentroidInternal = false; - break; + bool isPatchIntersected = CGElem::intersectBoxBox(patchBoxMin, patchBoxMax, objectBoxMin, objectBoxMax, 3, distanceTolerance); + + // Detect external cells + VolumeKernel::CellConstIterator cellBegin = m_mesh->cellConstBegin(); + VolumeKernel::CellConstIterator cellEnd = m_mesh->cellConstEnd(); + + m_nExternal = 0; + for (VolumeKernel::CellConstIterator itr = cellBegin; itr != cellEnd; ++itr) { + // Cells inside the narrowband cannot be external + long cellId = itr.getId(); + if (object->isInNarrowBand(cellId)) { + continue; + } + + // Check if the centroid is inside the bounding box + // + // Cells with the centroid inside the bounding box of the object + // cannot be external cells + if (isPatchIntersected) { + double geometricTolerance = m_mesh->getTol(); + std::array cellCentroid = object->getKernel()->computeCellCentroid(cellId); + + bool isCentroidInternal = true; + for (int i = 0; i < 3; ++i) { + if (cellCentroid[i] < objectBoxMin[i] - geometricTolerance || cellCentroid[i] > objectBoxMin[i] + geometricTolerance) { + isCentroidInternal = false; + break; + } + } + + if (isCentroidInternal) { + continue; } } - if (isCentroidInternal) { + // Check if the cell is inside the bounding box of the object + std::array cellBoxMin; + std::array cellBoxMax; + m_mesh->evalCellBoundingBox(cellId, &cellBoxMin, &cellBoxMax); + + bool isCellIntersected = CGElem::intersectBoxBox(cellBoxMin, cellBoxMax, objectBoxMin, objectBoxMax, 3, distanceTolerance); + if (isCellIntersected) { continue; } - } - // Check if the cell is inside the bounding box of the object - std::array cellBoxMin; - std::array cellBoxMax; - m_mesh->evalCellBoundingBox(cellId, &cellBoxMin, &cellBoxMax); - - bool isCellIntersected = CGElem::intersectBoxBox(cellBoxMin, cellBoxMax, objectBoxMin, objectBoxMax, 3, distanceTolerance); - if (isCellIntersected) { - continue; + std::size_t cellRawId = itr.getRawIndex(); + m_propagationStates.rawAt(cellRawId) = STATE_EXTERNAL; + ++m_nExternal; } - std::size_t cellRawId = itr.getRawIndex(); - m_propagationStates.rawAt(cellRawId) = STATE_EXTERNAL; - ++m_nExternal; + // Initialize the sign of the external region + m_externalSign = LevelSetSignStorage::SIGN_UNDEFINED; + } else { + // All cells are external + m_nExternal = m_mesh->getCellCount(); + + // Initialize the sign of the external region + m_externalSign = LevelSetSignStorage::SIGN_POSITIVE; } + + // Update the number of cells waiting for sing propagation m_nWaiting -= m_nExternal; - } + } else { + // It's not possible to identify + m_nExternal = 0; - // Initialize the sign of the external region - m_externalSign = LevelSetSignStorage::SIGN_UNDEFINED; + // Initialize the sign of the external region + m_externalSign = LevelSetSignStorage::SIGN_UNDEFINED; + } } /*! From 6ade2165fc05dd1f78a0ffec5bdf2c3891da2652 Mon Sep 17 00:00:00 2001 From: Andrea Iob Date: Fri, 30 Jun 2023 20:03:57 +0200 Subject: [PATCH 2/4] levelset: generalize generation of immutable objects --- src/levelset/levelSetImmutableObject.hpp | 4 +-- src/levelset/levelSetImmutableObject.tpp | 37 ++++++++++++------------ src/levelset/levelSetObjectFactory.tpp | 12 ++++---- 3 files changed, 27 insertions(+), 26 deletions(-) diff --git a/src/levelset/levelSetImmutableObject.hpp b/src/levelset/levelSetImmutableObject.hpp index d5dca2db34..1b8b2e0152 100644 --- a/src/levelset/levelSetImmutableObject.hpp +++ b/src/levelset/levelSetImmutableObject.hpp @@ -26,7 +26,6 @@ #define __BITPIT_LEVELSET_IMMUTABLE_OBJECT_HPP__ #include "levelSetObject.hpp" -#include "levelSetBooleanObject.hpp" #include "levelSetCachedObject.hpp" #include "levelSetSignedObject.hpp" @@ -47,8 +46,9 @@ class LevelSetImmutableObject : public LevelSetObject, public LevelSetImmutableO { public: + template + LevelSetImmutableObject(object_t *object); LevelSetImmutableObject(LevelSetCachedObject *object); - LevelSetImmutableObject(LevelSetBooleanObject *object); LevelSetImmutableObject * clone() const override; diff --git a/src/levelset/levelSetImmutableObject.tpp b/src/levelset/levelSetImmutableObject.tpp index 3a213d4c93..d57232cd36 100644 --- a/src/levelset/levelSetImmutableObject.tpp +++ b/src/levelset/levelSetImmutableObject.tpp @@ -52,31 +52,17 @@ LevelSetImmutableObject::LevelSetImmutableObject(int id) /*! * Constructor. * - * \param[in,out] source is an object whose contents will be acquired by the - * immutable object. Since source contents will be moved into the object, on - * output the source will be in an undefined state - */ -template -LevelSetImmutableObject::LevelSetImmutableObject(LevelSetCachedObject *source) - : LevelSetObject(*source), - LevelSetCachedObjectInterface(std::move(*source)), - LevelSetSignedObjectInterface(std::move(*source)) -{ -} - -/*! - * Constructor. - * - * \param[in,out] source is an object whose contents will be acquired by the +* \param[in,out] source is an object whose contents will be acquired by the * immutable object. Since source contents will be moved into the object, on * output the source will be in an undefined state */ template -LevelSetImmutableObject::LevelSetImmutableObject(LevelSetBooleanObject *source) +template +LevelSetImmutableObject::LevelSetImmutableObject(object_t *source) : LevelSetObject(*source) { // Mesh information - const VolumeKernel *mesh = (const_cast(source))->getKernel()->getMesh(); + const VolumeKernel *mesh = (const_cast(source))->getKernel()->getMesh(); VolumeKernel::CellConstIterator cellBegin = mesh->cellConstBegin(); VolumeKernel::CellConstIterator cellEnd = mesh->cellConstEnd(); @@ -116,6 +102,21 @@ LevelSetImmutableObject::LevelSetImmutableObject(LevelSetBo } } +/*! + * Constructor. + * + * \param[in,out] source is an object whose contents will be acquired by the + * immutable object. Since source contents will be moved into the object, on + * output the source will be in an undefined state + */ +template +LevelSetImmutableObject::LevelSetImmutableObject(LevelSetCachedObject *source) + : LevelSetObject(*source), + LevelSetCachedObjectInterface(std::move(*source)), + LevelSetSignedObjectInterface(std::move(*source)) +{ +} + /*! * Get LevelSetInfo of cell * @param[in] i cell idex diff --git a/src/levelset/levelSetObjectFactory.tpp b/src/levelset/levelSetObjectFactory.tpp index d9a75a81fd..6212b593df 100644 --- a/src/levelset/levelSetObjectFactory.tpp +++ b/src/levelset/levelSetObjectFactory.tpp @@ -87,18 +87,18 @@ std::unique_ptr LevelSetObjectFactory::_createImmutableObject(co switch (storageType) { case LevelSetStorageType::SPARSE: - if (LevelSetBooleanObject *booleanSource = dynamic_cast(source)) { - return std::unique_ptr(new LevelSetImmutableObject>(booleanSource)); - } else if (LevelSetCachedObject> *cachedSource = dynamic_cast> *>(source)) { + if (LevelSetCachedObject> *cachedSource = dynamic_cast> *>(source)) { return std::unique_ptr(new LevelSetImmutableObject>(cachedSource)); + } else { + return std::unique_ptr(new LevelSetImmutableObject>(source)); } break; case LevelSetStorageType::DENSE: - if (LevelSetBooleanObject *booleanSource = dynamic_cast(source)) { - return std::unique_ptr(new LevelSetImmutableObject>(booleanSource)); - } else if (LevelSetCachedObject> *cachedSource = dynamic_cast> *>(source)) { + if (LevelSetCachedObject> *cachedSource = dynamic_cast> *>(source)) { return std::unique_ptr(new LevelSetImmutableObject>(cachedSource)); + } else { + return std::unique_ptr(new LevelSetImmutableObject>(source)); } break; From 8df087f38e57ce96d08de6a836c834be2e19c7dc Mon Sep 17 00:00:00 2001 From: Andrea Iob Date: Fri, 30 Jun 2023 20:05:07 +0200 Subject: [PATCH 3/4] levelset: fix typo in output message --- src/levelset/levelSet.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/levelset/levelSet.cpp b/src/levelset/levelSet.cpp index e9df2f4630..9171f0ebd5 100644 --- a/src/levelset/levelSet.cpp +++ b/src/levelset/levelSet.cpp @@ -259,7 +259,7 @@ int LevelSet::addObject( LevelSetBooleanOperation operation, const std::vector &list, int id ) { - assert(m_kernel && " levelset: setMesh must be called befor adding a mask object "); + assert(m_kernel && " levelset: setMesh must be called before adding a mask object "); std::unique_ptr object = LevelSetObjectFactory::createMaskObject &, const VolumeKernel &>( m_kernel.get(), m_storageType, id, list, *m_kernel->getMesh() ) ; @@ -277,7 +277,7 @@ int LevelSet::addObject( const std::unordered_set &list, int id ) { */ int LevelSet::addObject( const std::vector &list, long refInterface, bool invert, int id ) { - assert(m_kernel && " levelset: setMesh must be called befor adding a mask object "); + assert(m_kernel && " levelset: setMesh must be called before adding a mask object "); std::unique_ptr object = LevelSetObjectFactory::createMaskObject &, long &, bool &, const VolumeKernel &>( m_kernel.get(), m_storageType, id, list, refInterface, invert, *m_kernel->getMesh() ) ; From d1f0a22f57c41f0cf330938e73620485648aaa14 Mon Sep 17 00:00:00 2001 From: Andrea Iob Date: Fri, 30 Jun 2023 20:04:50 +0200 Subject: [PATCH 4/4] levelset: allow to evaluate the complement of an object --- src/levelset/bitpit_levelset.hpp | 1 + src/levelset/levelSet.cpp | 18 + src/levelset/levelSet.hpp | 2 + src/levelset/levelSetComplementObject.cpp | 229 +++++++ src/levelset/levelSetComplementObject.hpp | 74 +++ src/levelset/levelSetObjectFactory.hpp | 4 + src/levelset/levelSetObjectFactory.tpp | 16 + .../integration_tests/levelset/CMakeLists.txt | 1 + .../levelset/test_levelset_00009.cpp | 600 ++++++++++++++++++ 9 files changed, 945 insertions(+) create mode 100644 src/levelset/levelSetComplementObject.cpp create mode 100644 src/levelset/levelSetComplementObject.hpp create mode 100644 test/integration_tests/levelset/test_levelset_00009.cpp diff --git a/src/levelset/bitpit_levelset.hpp b/src/levelset/bitpit_levelset.hpp index 7c7978f32d..cea7d313ca 100644 --- a/src/levelset/bitpit_levelset.hpp +++ b/src/levelset/bitpit_levelset.hpp @@ -46,6 +46,7 @@ #include "levelSetSegmentationObject.hpp" #include "levelSetSignedObject.hpp" #include "levelSetBooleanObject.hpp" +#include "levelSetComplementObject.hpp" #include "levelSetMaskObject.hpp" #include "levelSet.hpp" diff --git a/src/levelset/levelSet.cpp b/src/levelset/levelSet.cpp index 9171f0ebd5..de2e86b660 100644 --- a/src/levelset/levelSet.cpp +++ b/src/levelset/levelSet.cpp @@ -142,6 +142,24 @@ void LevelSet::setMesh( VolumeKernel* mesh ) { } +/*! + * Adds the complement of the specified object. + * Objects can be added to the levelset only after setting the mesh. + * @param[in] sourceId id of source object + * @param[in] id id to be assigned to object. In case default value is passed the insertion order will be used as identifier + * @return identifier of new object + */ +int LevelSet::addObjectComplement( int sourceId, int id ) { + + assert(m_kernel && " levelset: setMesh must be called before adding a mask object "); + + LevelSetObject *sourceObject = m_objects.at(sourceId).get() ; + + std::unique_ptr object = LevelSetObjectFactory::createComplementObject( m_kernel.get(), m_storageType, id, sourceObject ) ; + + return registerObject(std::move(object)); +}; + /*! * Adds a segmentation object * Objects can be added to the levelset only after setting the mesh. diff --git a/src/levelset/levelSet.hpp b/src/levelset/levelSet.hpp index d147d674f8..ea908d3c09 100644 --- a/src/levelset/levelSet.hpp +++ b/src/levelset/levelSet.hpp @@ -88,6 +88,8 @@ class LevelSet{ void setMesh( VolumeKernel* ) ; + int addObjectComplement( int, int id=levelSetDefaults::OBJECT ) ; + int addObject( std::unique_ptr &&, double, int id = levelSetDefaults::OBJECT ) ; int addObject( SurfaceKernel *, double, int id = levelSetDefaults::OBJECT ) ; int addObject( std::unique_ptr &&, double, int id = levelSetDefaults::OBJECT ) ; diff --git a/src/levelset/levelSetComplementObject.cpp b/src/levelset/levelSetComplementObject.cpp new file mode 100644 index 0000000000..67658c31a1 --- /dev/null +++ b/src/levelset/levelSetComplementObject.cpp @@ -0,0 +1,229 @@ +/*---------------------------------------------------------------------------*\ + * + * bitpit + * + * Copyright (C) 2015-2021 OPTIMAD engineering Srl + * + * ------------------------------------------------------------------------- + * License + * This file is part of bitpit. + * + * bitpit is free software: you can redistribute it and/or modify it + * under the terms of the GNU Lesser General Public License v3 (LGPL) + * as published by the Free Software Foundation. + * + * bitpit is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + * License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with bitpit. If not, see . + * +\*---------------------------------------------------------------------------*/ + +# include + +# include "bitpit_IO.hpp" +# include "bitpit_common.hpp" + +# include "levelSetObject.hpp" +# include "levelSetProxyObject.hpp" +# include "levelSetComplementObject.hpp" + +namespace bitpit { + +/*! + * \class LevelSetComplementObject + * \ingroup levelset + * + * \brief Class that allows to evaluate the complement of a LevelSetObjects. + */ + +/*! + * Constructor. + * + * \param[in] id identifier of object + * \param[in] source pointer to source object + */ +LevelSetComplementObject::LevelSetComplementObject(int id, const LevelSetObject *source) + : LevelSetProxyObject(id), + m_sourceObject(source) +{ +} + +/*! + * Copy constructor. + * + * Assigns same id to new object; + * \param[in] other object to be copied + */ +LevelSetComplementObject::LevelSetComplementObject(const LevelSetComplementObject &other) + : LevelSetProxyObject(other), + m_sourceObject(other.m_sourceObject) +{ +} + +/*! + * Returns LevelSetInfo. + * + * \param[in] id cell id + * \return LevelSetInfo +*/ +LevelSetInfo LevelSetComplementObject::getLevelSetInfo(long id) const +{ + return LevelSetInfo(getValue(id), getGradient(id)); +} + +/*! + * Get the levelset value. + * + * \param[in] id cell id + * \return levelset value in cell + */ +double LevelSetComplementObject::getLS(long id) const +{ + return getValue(id); +} + +/*! + * Get the levelset value. + * + * \param[in] id cell id + * \return levelset value in cell + */ +double LevelSetComplementObject::getValue(long id) const +{ + return (- m_sourceObject->getValue(id)); +} + +/*! + * Get the levelset gradient. + * + * \param[in] id cell id + * \return levelset gradient in cell + */ +std::array LevelSetComplementObject::getGradient(long id) const +{ + return (- 1. * m_sourceObject->getGradient(id)); +} + +/*! + * Computes the LevelSetInfo in a point. + * + * \param[in] coords point coordinates + * \return LevelSetInfo +*/ +LevelSetInfo LevelSetComplementObject::computeLevelSetInfo(const std::array &coords) const +{ + LevelSetInfo levelSetInfo = m_sourceObject->computeLevelSetInfo(coords); + levelSetInfo.value *= -1.; + levelSetInfo.gradient *= -1.; + + return levelSetInfo; +} + +/*! + * Replace a source object. + * + * \param[in] current current source object + * \param[in] updated updated source object + */ +void LevelSetComplementObject::replaceSourceObject(const LevelSetObject *current, const LevelSetObject *updated) +{ + if (current != m_sourceObject) { + throw std::runtime_error("Unable to find the source that should be replaced."); + } + + m_sourceObject = updated; +} + +/*! + * Clones the object. + * + * \return pointer to cloned object + */ +LevelSetComplementObject* LevelSetComplementObject::clone() const +{ + return new LevelSetComplementObject(*this); +} + +/*! + * Gets the surface normal at the projection point. + * + * \param[in] id cell index + * \return closest part + */ +std::array LevelSetComplementObject::getNormal(long id) const +{ + return (-1. * m_sourceObject->getNormal(id)); +} + +/*! + * Gets the closest part index. + * + * \param[in] id cell index + * \return closest part + */ +int LevelSetComplementObject::getPart(long id) const +{ + return m_sourceObject->getPart(id); +} + +/*! + * Get surface feature size. + * + * \param[in] id cell index + * \return characteristic size + */ +double LevelSetComplementObject::getSurfaceFeatureSize(long id) const +{ + return m_sourceObject->getSurfaceFeatureSize(id); +} + +/*! + * Get the smallest surface feature size. + * + * \return characteristic size + */ +double LevelSetComplementObject::getMinSurfaceFeatureSize() const +{ + return m_sourceObject->getMinSurfaceFeatureSize(); +} + +/*! + * Get the largest surface feature size. + * + * \return characteristic size + */ +double LevelSetComplementObject::getMaxSurfaceFeatureSize() const +{ + return m_sourceObject->getMaxSurfaceFeatureSize(); +} + +/*! + * Get the object that defines the levelset information for the specified cell. + * + * \param[in] id cell index + * \return The object that defines the levelset information for the specified + * cell. + */ +const LevelSetObject * LevelSetComplementObject::getReferenceObject(long id) const +{ + BITPIT_UNUSED(id); + + return m_sourceObject; +} + +/*! + * Get all objects that compose the boolean object. + * + * \return pointers to all primary objects involved in the definition of the + * boolean object levelset information. + */ +std::vector LevelSetComplementObject::getSourceObjects() const +{ + return std::vector(1, m_sourceObject); +} + +} diff --git a/src/levelset/levelSetComplementObject.hpp b/src/levelset/levelSetComplementObject.hpp new file mode 100644 index 0000000000..7d2fdad72a --- /dev/null +++ b/src/levelset/levelSetComplementObject.hpp @@ -0,0 +1,74 @@ +/*---------------------------------------------------------------------------*\ + * + * bitpit + * + * Copyright (C) 2015-2021 OPTIMAD engineering Srl + * + * ------------------------------------------------------------------------- + * License + * This file is part of bitpit. + * + * bitpit is free software: you can redistribute it and/or modify it + * under the terms of the GNU Lesser General Public License v3 (LGPL) + * as published by the Free Software Foundation. + * + * bitpit is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + * License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with bitpit. If not, see . + * +\*---------------------------------------------------------------------------*/ + +# ifndef __BITPIT_LEVELSET_COMPLEMENT_OBJECT_HPP__ +# define __BITPIT_LEVELSET_COMPLEMENT_OBJECT_HPP__ + +# include + +# include "levelSetCommon.hpp" +# include "levelSetProxyObject.hpp" + +namespace bitpit{ + +class LevelSetObject ; +class LevelSetProxyObject ; + +class LevelSetComplementObject: public LevelSetProxyObject { + + private: + const LevelSetObject* m_sourceObject; /**< Pointers to source object */ + + protected: + void replaceSourceObject(const LevelSetObject *current, const LevelSetObject *updated) override ; + + public: + LevelSetComplementObject(int id, const LevelSetObject*); + LevelSetComplementObject(const LevelSetComplementObject &); + + LevelSetComplementObject* clone() const override; + + LevelSetInfo getLevelSetInfo(long ) const override; + double getLS(long ) const override; + double getValue(long ) const override; + std::array getGradient(long ) const override; + + std::array getNormal(long ) const override; + int getPart(long ) const override; + double getSurfaceFeatureSize(long ) const override; + + LevelSetInfo computeLevelSetInfo(const std::array &) const override; + + double getMinSurfaceFeatureSize() const override; + double getMaxSurfaceFeatureSize() const override; + + const LevelSetObject * getReferenceObject( long ) const override; + + std::vector getSourceObjects() const override; + +}; + +} + +#endif diff --git a/src/levelset/levelSetObjectFactory.hpp b/src/levelset/levelSetObjectFactory.hpp index 8465e0a3aa..5f38c2f69f 100644 --- a/src/levelset/levelSetObjectFactory.hpp +++ b/src/levelset/levelSetObjectFactory.hpp @@ -27,6 +27,7 @@ #include "levelSetBooleanObject.hpp" #include "levelSetCartesianKernel.hpp" +#include "levelSetComplementObject.hpp" #include "levelSetImmutableObject.hpp" #include "levelSetMaskObject.hpp" #include "levelSetOctreeKernel.hpp" @@ -40,6 +41,9 @@ class LevelSetObjectFactory { template static std::unique_ptr createBooleanObject(const LevelSetKernel *kernel, LevelSetStorageType storageType, Args&&... args); + template + static std::unique_ptr createComplementObject(const LevelSetKernel *kernel, LevelSetStorageType storageType, Args&&... args); + template class narrow_band_cache_t> static std::unique_ptr createImmutableObject(const LevelSetKernel *kernel, LevelSetStorageType storageType, LevelSetObject *source); diff --git a/src/levelset/levelSetObjectFactory.tpp b/src/levelset/levelSetObjectFactory.tpp index 6212b593df..d0ab9e26fe 100644 --- a/src/levelset/levelSetObjectFactory.tpp +++ b/src/levelset/levelSetObjectFactory.tpp @@ -50,6 +50,22 @@ std::unique_ptr LevelSetObjectFactory::createBooleanObject(const return std::unique_ptr(new LevelSetBooleanObject(std::forward(args)...)); } +/*! + * Create a new complement object for the specified kernel. + * + * \param kernel is the kernel + * \param storageType is the storage type + * \param args are the arguments that will be forwarded to the constructor + */ +template +std::unique_ptr LevelSetObjectFactory::createComplementObject(const LevelSetKernel *kernel, LevelSetStorageType storageType, Args&&... args) +{ + BITPIT_UNUSED(kernel); + BITPIT_UNUSED(storageType); + + return std::unique_ptr(new LevelSetComplementObject(std::forward(args)...)); +} + /*! * Create a new immutable object for the specified kernel. * diff --git a/test/integration_tests/levelset/CMakeLists.txt b/test/integration_tests/levelset/CMakeLists.txt index df6b85c496..4627d0cb1f 100644 --- a/test/integration_tests/levelset/CMakeLists.txt +++ b/test/integration_tests/levelset/CMakeLists.txt @@ -35,6 +35,7 @@ list(APPEND TESTS "test_levelset_00005") list(APPEND TESTS "test_levelset_00006") list(APPEND TESTS "test_levelset_00007") list(APPEND TESTS "test_levelset_00008") +list(APPEND TESTS "test_levelset_00009") if (BITPIT_ENABLE_MPI) list(APPEND TESTS "test_levelset_parallel_00001:3") list(APPEND TESTS "test_levelset_parallel_00002:3") diff --git a/test/integration_tests/levelset/test_levelset_00009.cpp b/test/integration_tests/levelset/test_levelset_00009.cpp new file mode 100644 index 0000000000..45b8073423 --- /dev/null +++ b/test/integration_tests/levelset/test_levelset_00009.cpp @@ -0,0 +1,600 @@ +/*---------------------------------------------------------------------------*\ + * + * bitpit + * + * Copyright (C) 2015-2021 OPTIMAD engineering Srl + * + * ------------------------------------------------------------------------- + * License + * This file is part of bitpit. + * + * bitpit is free software: you can redistribute it and/or modify it + * under the terms of the GNU Lesser General Public License v3 (LGPL) + * as published by the Free Software Foundation. + * + * bitpit is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public + * License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with bitpit. If not, see . + * +\*---------------------------------------------------------------------------*/ + +//Standard Template Library +# include +# include + +#if BITPIT_ENABLE_MPI==1 +# include +#endif + +// bitpit +# include "bitpit_surfunstructured.hpp" +# include "bitpit_volcartesian.hpp" +# include "bitpit_voloctree.hpp" +# include "bitpit_volunstructured.hpp" +# include "bitpit_levelset.hpp" + +using namespace bitpit; + +/*! + * Generate segmentation. + * + * \result The generated segmentation. + */ +std::unique_ptr generateSegmentation() +{ + // Input geometry +#if BITPIT_ENABLE_MPI + std::unique_ptr segmentation( new SurfUnstructured(2, MPI_COMM_NULL) ); +#else + std::unique_ptr segmentation( new SurfUnstructured(2) ); +#endif + + segmentation->importSTL("./data/cube.stl", true); + + segmentation->deleteCoincidentVertices() ; + segmentation->initializeAdjacencies() ; + + segmentation->getVTK().setName("geometry_002") ; + segmentation->write() ; + + return segmentation; +} + +/*! + * Generate the Cartesian mesh. + * + * \result The generated Cartesian mesh. + */ +std::unique_ptr generateCartesianMesh(const SurfUnstructured &segmentation) +{ + int dimensions = 3; + + std::array meshMin, meshMax; + segmentation.getBoundingBox(meshMin, meshMax); + + std::array delta = meshMax - meshMin; + meshMin -= 0.1 * delta; + meshMax += 0.1 * delta; + delta = meshMax - meshMin; + + std::array nc = {{64, 64, 64}}; + + std::unique_ptr mesh(new VolCartesian(dimensions, meshMin, delta, nc)); + + return mesh; +} + +/*! + * Generate the Octree mesh. + * + * \result The generated Octree mesh. + */ +std::unique_ptr generateOctreeMesh(const SurfUnstructured &segmentation) +{ + int dimensions = 3; + + std::array segmentationMin; + std::array segmentationMax; + segmentation.getBoundingBox(segmentationMin, segmentationMax); + + std::array delta = segmentationMax - segmentationMin; + segmentationMin -= 0.1 * delta; + segmentationMax += 0.1 * delta; + delta = segmentationMax - segmentationMin; + + std::array origin = segmentationMin; + + double length = 0.; + for (int i = 0; i < 3; ++i) { + length = std::max(length, segmentationMax[i] - segmentationMin[i]); + }; + + double dh = length / 64; + +#if BITPIT_ENABLE_MPI + std::unique_ptr mesh(new VolOctree(dimensions, origin, length, dh, MPI_COMM_NULL)); +#else + std::unique_ptr mesh(new VolOctree(dimensions, origin, length, dh)); +#endif + + return mesh; +} + +/*! + * Generate the Unstructured mesh. + * + * \result The generated Unstructured mesh. + */ +std::unique_ptr generateUnstructuredMesh(const SurfUnstructured &segmentation) +{ + int dimensions = 3; + + std::array segmentationMin; + std::array segmentationMax; + segmentation.getBoundingBox(segmentationMin, segmentationMax); + + std::array delta = segmentationMax - segmentationMin; + segmentationMin -= 0.1 * delta; + segmentationMax += 0.1 * delta; + delta = segmentationMax - segmentationMin; + + std::array nCells = {{64, 64, 64}}; + std::array nVertices = {{nCells[0] + 1, nCells[1] + 1, nCells[2] + 1}}; + + // Create patch +#if BITPIT_ENABLE_MPI + std::unique_ptr mesh(new VolUnstructured(dimensions, MPI_COMM_NULL)); +#else + std::unique_ptr mesh(new VolUnstructured(dimensions)); +#endif + mesh->setVertexAutoIndexing(false); + + // Create vertices + for (int i = 0; i < nVertices[0]; ++i) { + double x = segmentationMin[0] + delta[0] / nCells[0] * i; + for (int j = 0; j < nVertices[1]; ++j) { + double y = segmentationMin[1] + delta[1] / nCells[1] * j; + for (int k = 0; k < nVertices[2]; ++k) { + double z = segmentationMin[2] + delta[2] / nCells[2] * k; + + long vertexId = i + nVertices[0] * j + nVertices[0] * nVertices[1] * k; + + mesh->addVertex({{x, y, z}}, vertexId); + + } + } + } + + // Create cells + std::unordered_set customCellIds; + customCellIds.insert(171039); + customCellIds.insert(187359); + + int cellConnectSize = ReferenceElementInfo::MAX_ELEM_VERTICES; + std::vector cellConnect(cellConnectSize); + for (int i = 0; i < nCells[0]; ++i) { + for (int j = 0; j < nCells[1]; ++j) { + for (int k = 0; k < nCells[2]; ++k) { + long cellId = i + nCells[0] * j + nCells[0] * nCells[1] * k; + if (customCellIds.count(cellId) != 0) { + continue; + } + + cellConnect[0] = i + nVertices[0] * j + nVertices[0] * nVertices[1] * k; + cellConnect[1] = (i + 1) + nVertices[0] * j + nVertices[0] * nVertices[1] * k; + cellConnect[2] = (i + 1) + nVertices[0] * (j + 1) + nVertices[0] * nVertices[1] * k; + cellConnect[3] = i + nVertices[0] * (j + 1) + nVertices[0] * nVertices[1] * k; + cellConnect[4] = i + nVertices[0] * j + nVertices[0] * nVertices[1] * (k + 1); + cellConnect[5] = (i + 1) + nVertices[0] * j + nVertices[0] * nVertices[1] * (k + 1); + cellConnect[6] = (i + 1) + nVertices[0] * (j + 1) + nVertices[0] * nVertices[1] * (k + 1); + cellConnect[7] = i + nVertices[0] * (j + 1) + nVertices[0] * nVertices[1] * (k + 1); + + mesh->addCell(ElementType::HEXAHEDRON, cellConnect, cellId); + } + } + } + + cellConnect[0] = 176377; + cellConnect[1] = 176442; + cellConnect[2] = 180602; + cellConnect[3] = 180667; + cellConnect[4] = 176376; + cellConnect[5] = 176441; + cellConnect[6] = 180601; + cellConnect[7] = 180666 ; + mesh->addCell(ElementType::VOXEL, cellConnect); + + std::vector faceStream(31); + faceStream[ 0] = 6; + faceStream[ 1] = 4; + faceStream[ 2] = 197437; + faceStream[ 3] = 193212; + faceStream[ 4] = 193277; + faceStream[ 5] = 197502; + faceStream[ 6] = 4; + faceStream[ 7] = 193212; + faceStream[ 8] = 193211; + faceStream[ 9] = 193276; + faceStream[10] = 193277; + faceStream[11] = 4; + faceStream[12] = 197436; + faceStream[13] = 197437; + faceStream[14] = 197502; + faceStream[15] = 197501; + faceStream[16] = 4; + faceStream[17] = 197502; + faceStream[18] = 193277; + faceStream[19] = 193276; + faceStream[20] = 197501; + faceStream[21] = 4; + faceStream[22] = 197436; + faceStream[23] = 193211; + faceStream[24] = 193212; + faceStream[25] = 197437; + faceStream[26] = 4; + faceStream[27] = 193211; + faceStream[28] = 197436; + faceStream[29] = 197501; + faceStream[30] = 193276; + mesh->addCell(ElementType::POLYHEDRON, faceStream); + + return mesh; +} + +/*! +* Subtest 001 +* +* Testing evaluation of complement object of a 3D levelset on a Cartesian mesh in default memory mode. +*/ +int subtest_001() +{ + log::cout() << std::endl; + log::cout() << "Testing three-dimensional levelset on a Cartesian mesh in default memory mode" << std::endl; + + // Input geometry + log::cout() << " - Loading geometry" << std::endl; + + std::unique_ptr segmentation = generateSegmentation(); + segmentation->getVTK().setName("geometry_002"); + segmentation->write(); + + log::cout() << "n. vertex: " << segmentation->getVertexCount() << std::endl; + log::cout() << "n. simplex: " << segmentation->getCellCount() << std::endl; + + // Create the mesh + log::cout() << " - Setting mesh" << std::endl; + + std::unique_ptr mesh = generateCartesianMesh(*segmentation); + mesh->switchMemoryMode(VolCartesian::MEMORY_NORMAL); + mesh->initializeAdjacencies(); + mesh->update(); + + // Initialize levelset + log::cout() << " - Initializing levelset" << std::endl; + + LevelSet levelset ; + levelset.setPropagateSign(true); + levelset.setMesh(mesh.get()); + + int objectId = levelset.addObject(segmentation.get(), BITPIT_PI); + int complementId = levelset.addObjectComplement(objectId); + + // Compute levelset + log::cout() << " - Evaluating levelset" << std::endl; + + std::chrono::time_point start = std::chrono::system_clock::now(); + + levelset.compute( ) ; + + std::chrono::time_point end = std::chrono::system_clock::now(); + int elapsed_seconds = std::chrono::duration_cast(end-start).count(); + log::cout() << "elapsed time: " << elapsed_seconds << " ms" << std::endl; + + // Write output + log::cout() << " - Writing output" << std::endl; + + levelset.getObject(objectId).enableVTKOutput(LevelSetWriteField::VALUE); + levelset.getObject(complementId).enableVTKOutput(LevelSetWriteField::VALUE); + + mesh->getVTK().setName("levelset_009_cartesian_normal"); + mesh->write(); + + // Check if the levelset is correct + log::cout() << " - Checking results" << std::endl; + + std::vector testCells{0, 111462, 128486, 191398}; + for (long cellId : testCells) { + double objectLevelset = levelset.getObject(objectId).getValue(cellId); + double complementLevelset = levelset.getObject(complementId).getValue(cellId); + + log::cout() << " - Ojbect levelset on cell " << cellId << " = " << objectLevelset << std::endl; + log::cout() << " - Complement levelset on cell " << cellId << " = " << complementLevelset << std::endl; + if (!utils::DoubleFloatingEqual()(objectLevelset + complementLevelset, 0.)) { + log::cout() << " Expected complement levelset on cell " << cellId << " = " << (- objectLevelset) << std::endl; + throw std::runtime_error("Error in evaluation of levelset of complement object."); + } + } + + return 0; +} + +/*! +* Subtest 002 +* +* Testing evaluation of complement object of a 3D levelset on a Cartesian mesh in light memory mode. +*/ +int subtest_002() +{ + log::cout() << std::endl; + log::cout() << "Testing three-dimensional levelset on a Cartesian mesh in light memory mode" << std::endl; + + // Input geometry + log::cout() << " - Loading geometry" << std::endl; + + std::unique_ptr segmentation = generateSegmentation(); + segmentation->getVTK().setName("geometry_002"); + segmentation->write(); + + log::cout() << "n. vertex: " << segmentation->getVertexCount() << std::endl; + log::cout() << "n. simplex: " << segmentation->getCellCount() << std::endl; + + // Create the mesh + log::cout() << " - Setting mesh" << std::endl; + + std::unique_ptr mesh = generateCartesianMesh(*segmentation); + mesh->switchMemoryMode(VolCartesian::MEMORY_LIGHT); + + // Initialize levelset + log::cout() << " - Initializing levelset" << std::endl; + + LevelSet levelset ; + levelset.setPropagateSign(false); + levelset.setMesh(mesh.get()); + + int objectId = levelset.addObject(segmentation.get(), BITPIT_PI); + int complementId = levelset.addObjectComplement(objectId); + + // Compute levelset + log::cout() << " - Evaluating levelset" << std::endl; + + std::chrono::time_point start = std::chrono::system_clock::now(); + + levelset.compute( ) ; + + std::chrono::time_point end = std::chrono::system_clock::now(); + int elapsed_seconds = std::chrono::duration_cast(end-start).count(); + log::cout() << "elapsed time: " << elapsed_seconds << " ms" << std::endl; + + // Write output + log::cout() << " - Writing output" << std::endl; + + levelset.getObject(objectId).enableVTKOutput(LevelSetWriteField::VALUE); + levelset.getObject(complementId).enableVTKOutput(LevelSetWriteField::VALUE); + + mesh->switchMemoryMode(VolCartesian::MEMORY_NORMAL); + mesh->getVTK().setName("levelset_009_cartesian_light"); + mesh->write(); + + // Check if the levelset is correct + log::cout() << " - Checking results" << std::endl; + + std::vector testCells{0, 111462, 128486, 191398}; + for (long cellId : testCells) { + double objectLevelset = levelset.getObject(objectId).getValue(cellId); + double complementLevelset = levelset.getObject(complementId).getValue(cellId); + + log::cout() << " - Ojbect levelset on cell " << cellId << " = " << objectLevelset << std::endl; + log::cout() << " - Complement levelset on cell " << cellId << " = " << complementLevelset << std::endl; + if (!utils::DoubleFloatingEqual()(objectLevelset + complementLevelset, 0.)) { + log::cout() << " Expected complement levelset on cell " << cellId << " = " << (- objectLevelset) << std::endl; + throw std::runtime_error("Error in evaluation of levelset of complement object."); + } + } + + return 0; +} + +/*! +* Subtest 003 +* +* Testing evaluation of complement object of a 3D levelset on an Octreee mesh. +*/ +int subtest_003() +{ + log::cout() << std::endl; + log::cout() << "Testing three-dimensional levelset on an Octree mesh" << std::endl; + + // Input geometry + log::cout() << " - Loading geometry" << std::endl; + + std::unique_ptr segmentation = generateSegmentation(); + segmentation->getVTK().setName("geometry_002"); + segmentation->write(); + + log::cout() << "n. vertex: " << segmentation->getVertexCount() << std::endl; + log::cout() << "n. simplex: " << segmentation->getCellCount() << std::endl; + + // Create the mesh + log::cout() << " - Setting mesh" << std::endl; + + std::unique_ptr mesh = generateOctreeMesh(*segmentation); + mesh->initializeAdjacencies(); + mesh->update(); + + // Initialize levelset + log::cout() << " - Initializing levelset" << std::endl; + + LevelSet levelset ; + levelset.setPropagateSign(true); + levelset.setMesh(mesh.get()); + + int objectId = levelset.addObject(segmentation.get(), BITPIT_PI); + int complementId = levelset.addObjectComplement(objectId); + + // Compute levelset + log::cout() << " - Evaluating levelset" << std::endl; + + std::chrono::time_point start = std::chrono::system_clock::now(); + + levelset.compute( ) ; + + std::chrono::time_point end = std::chrono::system_clock::now(); + int elapsed_seconds = std::chrono::duration_cast(end-start).count(); + log::cout() << "elapsed time: " << elapsed_seconds << " ms" << std::endl; + + // Write output + log::cout() << " - Writing output" << std::endl; + + levelset.getObject(objectId).enableVTKOutput(LevelSetWriteField::VALUE); + levelset.getObject(complementId).enableVTKOutput(LevelSetWriteField::VALUE); + + mesh->getVTK().setName("levelset_009_octree"); + mesh->write(); + + // Check if the levelset is correct + log::cout() << " - Checking results" << std::endl; + + std::vector testCells{0, 136783, 145227, 204399}; + for (long cellId : testCells) { + double objectLevelset = levelset.getObject(objectId).getValue(cellId); + double complementLevelset = levelset.getObject(complementId).getValue(cellId); + + log::cout() << " - Ojbect levelset on cell " << cellId << " = " << objectLevelset << std::endl; + log::cout() << " - Complement levelset on cell " << cellId << " = " << complementLevelset << std::endl; + if (!utils::DoubleFloatingEqual()(objectLevelset + complementLevelset, 0.)) { + log::cout() << " Expected complement levelset on cell " << cellId << " = " << (- objectLevelset) << std::endl; + throw std::runtime_error("Error in evaluation of levelset of complement object."); + } + } + + return 0; +} + +/*! +* Subtest 004 +* +* Testing evaluation of complement object of a 3D levelset on an Unstructured mesh. +*/ +int subtest_004() +{ + log::cout() << std::endl; + log::cout() << "Testing three-dimensional levelset on an Unstructured mesh" << std::endl; + + // Input geometry + log::cout() << " - Loading geometry" << std::endl; + + std::unique_ptr segmentation = generateSegmentation(); + segmentation->getVTK().setName("geometry_002"); + segmentation->write(); + + log::cout() << "n. vertex: " << segmentation->getVertexCount() << std::endl; + log::cout() << "n. simplex: " << segmentation->getCellCount() << std::endl; + + // Create the mesh + log::cout() << " - Setting mesh" << std::endl; + + std::unique_ptr mesh = generateUnstructuredMesh(*segmentation); + mesh->initializeAdjacencies(); + mesh->update(); + + // Initialize levelset + log::cout() << " - Initializing levelset" << std::endl; + + LevelSet levelset ; + levelset.setPropagateSign(true); + levelset.setMesh(mesh.get()); + + int objectId = levelset.addObject(segmentation.get(), BITPIT_PI); + int complementId = levelset.addObjectComplement(objectId); + + // Compute levelset + log::cout() << " - Evaluating levelset" << std::endl; + + std::chrono::time_point start = std::chrono::system_clock::now(); + + levelset.compute( ) ; + + std::chrono::time_point end = std::chrono::system_clock::now(); + int elapsed_seconds = std::chrono::duration_cast(end-start).count(); + log::cout() << "elapsed time: " << elapsed_seconds << " ms" << std::endl; + + // Write output + log::cout() << " - Writing output" << std::endl; + + levelset.getObject(objectId).enableVTKOutput(LevelSetWriteField::VALUE); + levelset.getObject(complementId).enableVTKOutput(LevelSetWriteField::VALUE); + + mesh->getVTK().setName("levelset_009_unstructured"); + mesh->write(); + + // Check if the levelset is correct + log::cout() << " - Checking results" << std::endl; + + std::vector testCells{0, 87014, 187302, 145958}; + for (long cellId : testCells) { + double objectLevelset = levelset.getObject(objectId).getValue(cellId); + double complementLevelset = levelset.getObject(complementId).getValue(cellId); + + log::cout() << " - Ojbect levelset on cell " << cellId << " = " << objectLevelset << std::endl; + log::cout() << " - Complement levelset on cell " << cellId << " = " << complementLevelset << std::endl; + if (!utils::DoubleFloatingEqual()(objectLevelset + complementLevelset, 0.)) { + log::cout() << " Expected complement levelset on cell " << cellId << " = " << (- objectLevelset) << std::endl; + throw std::runtime_error("Error in evaluation of levelset of complement object."); + } + } + + return 0; +} + +/*! +* Main program. +*/ +int main(int argc, char *argv[]) +{ +#if BITPIT_ENABLE_MPI==1 + MPI_Init(&argc,&argv); +#else + BITPIT_UNUSED(argc); + BITPIT_UNUSED(argv); +#endif + + // Initialize the logger + log::manager().initialize(log::MODE_COMBINE); + + // Run the subtests + log::cout() << "Testing basic levelset features" << std::endl; + + int status; + try { + status = subtest_001(); + if (status != 0) { + return status; + } + + status = subtest_002(); + if (status != 0) { + return status; + } + + status = subtest_003(); + if (status != 0) { + return status; + } + + status = subtest_004(); + if (status != 0) { + return status; + } + } catch (const std::exception &exception) { + log::cout() << exception.what(); + exit(1); + } + +#if BITPIT_ENABLE_MPI==1 + MPI_Finalize(); +#endif +}