From 829228e90433c33c54ccc103fbf5dfff21cc70fa Mon Sep 17 00:00:00 2001 From: Felix S Date: Fri, 26 Apr 2024 19:34:46 +0200 Subject: [PATCH 1/6] Decouple cell set from psDomain, psVTKWriter and psSmartPointer --- include/viennaps/cellSet/csDenseCellSet.hpp | 51 +++++++++---------- include/viennaps/cellSet/csMeanFreePath.hpp | 28 +++++++--- include/viennaps/cellSet/csSegmentCells.hpp | 8 +-- include/viennaps/cellSet/csSurfaceCellSet.hpp | 10 ++-- 4 files changed, 54 insertions(+), 43 deletions(-) diff --git a/include/viennaps/cellSet/csDenseCellSet.hpp b/include/viennaps/cellSet/csDenseCellSet.hpp index 9ec5ddf..2f25b8d 100644 --- a/include/viennaps/cellSet/csDenseCellSet.hpp +++ b/include/viennaps/cellSet/csDenseCellSet.hpp @@ -6,14 +6,13 @@ #include "../psLogger.hpp" #include "../psMaterials.hpp" -#include "../psSmartPointer.hpp" -#include "../psVTKWriter.hpp" #include #include #include #include #include +#include #include #include @@ -24,15 +23,15 @@ */ template class csDenseCellSet { private: - using gridType = psSmartPointer>; + using gridType = lsSmartPointer>; using levelSetsType = - psSmartPointer>>>; - using materialMapType = psSmartPointer; + lsSmartPointer>>>; + using materialMapType = lsSmartPointer; levelSetsType levelSets = nullptr; gridType cellGrid = nullptr; - psSmartPointer> surface = nullptr; - psSmartPointer> BVH = nullptr; + lsSmartPointer> surface = nullptr; + lsSmartPointer> BVH = nullptr; materialMapType materialMap = nullptr; T gridDelta; @@ -67,25 +66,25 @@ template class csDenseCellSet { materialMap = passedMaterialMap; if (cellGrid == nullptr) - cellGrid = psSmartPointer>::New(); + cellGrid = lsSmartPointer>::New(); if (surface == nullptr) - surface = psSmartPointer>::New(levelSets->back()); + surface = lsSmartPointer>::New(levelSets->back()); else surface->deepCopy(levelSets->back()); gridDelta = surface->getGrid().getGridDelta(); depth = passedDepth; - std::vector>> levelSetsInOrder; - auto plane = psSmartPointer>::New(surface->getGrid()); + std::vector>> levelSetsInOrder; + auto plane = lsSmartPointer>::New(surface->getGrid()); { T origin[D] = {0.}; T normal[D] = {0.}; origin[D - 1] = depth; normal[D - 1] = 1.; lsMakeGeometry(plane, - psSmartPointer>::New(origin, normal)) + lsSmartPointer>::New(origin, normal)) .apply(); } if (!cellSetAboveSurface) @@ -103,12 +102,12 @@ template class csDenseCellSet { #ifndef NDEBUG int db_ls = 0; for (auto &ls : levelSetsInOrder) { - auto mesh = psSmartPointer>::New(); + auto mesh = lsSmartPointer>::New(); lsToSurfaceMesh(ls, mesh).apply(); - psVTKWriter(mesh, "cellSet_debug_" + std::to_string(db_ls++) + ".vtp") + lsVTKWriter(mesh, "cellSet_debug_" + std::to_string(db_ls++) + ".vtp") .apply(); } - psVTKWriter(cellGrid, "cellSet_debug_init.vtu").apply(); + lsVTKWriter(cellGrid, "cellSet_debug_init.vtu").apply(); #endif adjustMaterialIds(); @@ -139,7 +138,7 @@ template class csDenseCellSet { minExtent /= 2; } - BVH = psSmartPointer>::New(getBoundingBox(), BVHlayers); + BVH = lsSmartPointer>::New(getBoundingBox(), BVHlayers); buildBVH(); } @@ -192,9 +191,9 @@ template class csDenseCellSet { return cellGrid->template getElements<(1 << D)>()[idx]; } - psSmartPointer> getSurface() { return surface; } + lsSmartPointer> getSurface() { return surface; } - psSmartPointer> getCellGrid() { return cellGrid; } + lsSmartPointer> getCellGrid() { return cellGrid; } levelSetsType getLevelSets() const { return levelSets; } @@ -309,7 +308,7 @@ template class csDenseCellSet { // Write the cell set as .vtu file void writeVTU(std::string fileName) { - psVTKWriter(cellGrid, fileName).apply(); + lsVTKWriter(cellGrid, fileName).apply(); } // Save cell set data in simple text format @@ -403,15 +402,15 @@ template class csDenseCellSet { auto materialIds = getScalarData("Material"); // create overlay material - std::vector>> levelSetsInOrder; - auto plane = psSmartPointer>::New(surface->getGrid()); + std::vector>> levelSetsInOrder; + auto plane = lsSmartPointer>::New(surface->getGrid()); { T origin[D] = {0.}; T normal[D] = {0.}; origin[D - 1] = depth; normal[D - 1] = 1.; lsMakeGeometry(plane, - psSmartPointer>::New(origin, normal)) + lsSmartPointer>::New(origin, normal)) .apply(); } if (!cellSetAboveSurface) @@ -487,18 +486,18 @@ template class csDenseCellSet { // Updates the surface of the cell set. The new surface should be below the // old surface as this function can only remove cells from the cell set. void updateSurface() { - auto updateCellGrid = psSmartPointer>::New(); + auto updateCellGrid = lsSmartPointer>::New(); lsToVoxelMesh voxelConverter(updateCellGrid); { - auto plane = psSmartPointer>::New(surface->getGrid()); + auto plane = lsSmartPointer>::New(surface->getGrid()); T origin[D] = {0.}; T normal[D] = {0.}; origin[D - 1] = depth; normal[D - 1] = 1.; lsMakeGeometry(plane, - psSmartPointer>::New(origin, normal)) + lsSmartPointer>::New(origin, normal)) .apply(); voxelConverter.insertNextLevelSet(plane); } @@ -737,7 +736,7 @@ template class csDenseCellSet { } void calculateMinMaxIndex( - const std::vector>> &levelSetsInOrder) { + const std::vector>> &levelSetsInOrder) { // set to zero for (unsigned i = 0; i < D; ++i) { minIndex[i] = std::numeric_limits::max(); diff --git a/include/viennaps/cellSet/csMeanFreePath.hpp b/include/viennaps/cellSet/csMeanFreePath.hpp index 486fe97..954b36b 100644 --- a/include/viennaps/cellSet/csMeanFreePath.hpp +++ b/include/viennaps/cellSet/csMeanFreePath.hpp @@ -1,12 +1,20 @@ #pragma once -#include "psDomain.hpp" +#include "csDenseCellSet.hpp" + #include "psLogger.hpp" +#include +#include #include template class csMeanFreePath { +private: + using levelSetsType = + lsSmartPointer>>>; + using cellSetType = lsSmartPointer>; + public: csMeanFreePath() : traceDevice(rtcNewDevice("hugepages=1")) { static_assert(D == 2 && @@ -15,8 +23,12 @@ template class csMeanFreePath { ~csMeanFreePath() { rtcReleaseDevice(traceDevice); } - void setDomain(const psSmartPointer> passedDomain) { - domain = passedDomain; + void setLevelSets(levelSetsType passedLevelSets) { + levelSets = passedLevelSets; + } + + void setCellSet(cellSetType passedCellSet) { + cellSet = passedCellSet; } void setBulkLambda(const NumericType passedBulkLambda) { @@ -35,7 +47,7 @@ template class csMeanFreePath { void apply() { psLogger::getInstance().addInfo("Calculating mean free path ...").print(); - domain->getCellSet()->addScalarData("MeanFreePath", 0.); + cellSet->addScalarData("MeanFreePath", 0.); runKernel(); } @@ -47,7 +59,6 @@ template class csMeanFreePath { _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON); _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); #endif - auto &cellSet = domain->getCellSet(); unsigned numCells = cellSet->getElements().size(); auto data = cellSet->getScalarData("MeanFreePath"); auto materials = cellSet->getScalarData("Material"); @@ -55,11 +66,11 @@ template class csMeanFreePath { auto traceGeometry = rayGeometry(); { auto mesh = lsSmartPointer>::New(); - lsToDiskMesh(domain->getLevelSets()->back(), mesh) + lsToDiskMesh(levelSets->back(), mesh) .apply(); auto &points = mesh->getNodes(); auto normals = mesh->getCellData().getVectorData("Normals"); - gridDelta = domain->getGrid().getGridDelta(); + gridDelta = levelSets->back()->getGrid().getGridDelta(); traceGeometry.initGeometry(traceDevice, points, *normals, gridDelta * rayInternal::DiskFactor); } @@ -171,7 +182,8 @@ template class csMeanFreePath { } private: - psSmartPointer> domain = nullptr; + levelSetsType levelSets = nullptr; + cellSetType cellSet = nullptr; RTCDevice traceDevice; NumericType gridDelta = 0; diff --git a/include/viennaps/cellSet/csSegmentCells.hpp b/include/viennaps/cellSet/csSegmentCells.hpp index 2f4c48b..af626b0 100644 --- a/include/viennaps/cellSet/csSegmentCells.hpp +++ b/include/viennaps/cellSet/csSegmentCells.hpp @@ -3,23 +3,23 @@ #include "csDenseCellSet.hpp" template class csSegmentCells { - psSmartPointer> cellSet = nullptr; + lsSmartPointer> cellSet = nullptr; std::string cellTypeString = "CellType"; psMaterial bulkMaterial = psMaterial::GAS; public: csSegmentCells( - const psSmartPointer> &passedCellSet) + const lsSmartPointer> &passedCellSet) : cellSet(passedCellSet) {} csSegmentCells( - const psSmartPointer> &passedCellSet, + const lsSmartPointer> &passedCellSet, const std::string cellTypeString, const psMaterial passedBulkMaterial) : cellSet(passedCellSet), cellTypeString(cellTypeString), bulkMaterial(passedBulkMaterial) {} void setCellSet( - const psSmartPointer> &passedCellSet) { + const lsSmartPointer> &passedCellSet) { cellSet = passedCellSet; } diff --git a/include/viennaps/cellSet/csSurfaceCellSet.hpp b/include/viennaps/cellSet/csSurfaceCellSet.hpp index 599aa55..9c906a9 100644 --- a/include/viennaps/cellSet/csSurfaceCellSet.hpp +++ b/include/viennaps/cellSet/csSurfaceCellSet.hpp @@ -17,8 +17,8 @@ template class cellBase { }; template > class csSurfaceCellSet { - psSmartPointer> domain = nullptr; - psSmartPointer> BVH = nullptr; + lsSmartPointer> domain = nullptr; + lsSmartPointer> BVH = nullptr; std::vector surfaceCells; std::vector innerCells; @@ -26,7 +26,7 @@ template > class csSurfaceCellSet { public: csSurfaceCellSet() = default; - csSurfaceCellSet(const psSmartPointer> &passedDomain, + csSurfaceCellSet(const lsSmartPointer> &passedDomain, const T extraHeight = 0.) { domain = passedDomain; auto highestPoint = domain->getBoundingBox()[1][D - 1]; @@ -70,7 +70,7 @@ template > class csSurfaceCellSet { return surfaceCells[0]; } - psSmartPointer> getDomain() { return domain; } + lsSmartPointer> getDomain() { return domain; } void makeSurfaceCells(const std::vector &surfaceCellIds) { for (auto id : surfaceCellIds) { @@ -182,7 +182,7 @@ template > class csSurfaceCellSet { } BVH = - psSmartPointer>::New(cellSet->getBoundingBox(), BVHlayers); + lsSmartPointer>::New(cellSet->getBoundingBox(), BVHlayers); psUtils::Timer timer; timer.start(); From cdfd002744920d7203a662428dc1369b852c4303 Mon Sep 17 00:00:00 2001 From: Felix S Date: Sat, 4 May 2024 22:31:58 +0200 Subject: [PATCH 2/6] Decouple cell set from psMaterial (BROKEN Example oxideRegrowth) --- .../atomicLayerDeposition.cpp | 4 +++- include/viennaps/cellSet/csDenseCellSet.hpp | 21 +++++++++---------- include/viennaps/cellSet/csMeanFreePath.hpp | 9 ++++---- include/viennaps/cellSet/csSegmentCells.hpp | 16 +++++++------- include/viennaps/psDomain.hpp | 6 ++++-- 5 files changed, 31 insertions(+), 25 deletions(-) diff --git a/examples/atomicLayerDeposition/atomicLayerDeposition.cpp b/examples/atomicLayerDeposition/atomicLayerDeposition.cpp index 6a84599..681a615 100644 --- a/examples/atomicLayerDeposition/atomicLayerDeposition.cpp +++ b/examples/atomicLayerDeposition/atomicLayerDeposition.cpp @@ -42,7 +42,9 @@ int main(int argc, char *argv[]) { psMaterial::GAS, true); auto &cellSet = domain->getCellSet(); // Segment the cells into surface, material, and gas cells - csSegmentCells(cellSet).apply(); + csSegmentCells segmentation(cellSet); + segmentation.setBulkMaterial(psMaterial::GAS); + segmentation.apply(); cellSet->writeVTU("initial.vtu"); diff --git a/include/viennaps/cellSet/csDenseCellSet.hpp b/include/viennaps/cellSet/csDenseCellSet.hpp index 2f25b8d..1bec5bb 100644 --- a/include/viennaps/cellSet/csDenseCellSet.hpp +++ b/include/viennaps/cellSet/csDenseCellSet.hpp @@ -5,7 +5,6 @@ #include "csUtil.hpp" #include "../psLogger.hpp" -#include "../psMaterials.hpp" #include #include @@ -26,7 +25,7 @@ template class csDenseCellSet { using gridType = lsSmartPointer>; using levelSetsType = lsSmartPointer>>>; - using materialMapType = lsSmartPointer; + using materialMapType = lsSmartPointer; levelSetsType levelSets = nullptr; gridType cellGrid = nullptr; @@ -43,7 +42,7 @@ template class csDenseCellSet { hrleVectorType minIndex, maxIndex; bool cellSetAboveSurface = false; - psMaterial coverMaterial = psMaterial::None; + int coverMaterial = -1; std::bitset periodicBoundary; std::vector *fillingFractions; @@ -259,8 +258,9 @@ template class csDenseCellSet { cellSetAboveSurface = passedCellSetPosition; } - void setCoverMaterial(const psMaterial passedCoverMaterial) { - coverMaterial = passedCoverMaterial; + template + void setCoverMaterial(const Material passedCoverMaterial) { + coverMaterial = static_cast(passedCoverMaterial); } bool getCellSetPosition() const { return cellSetAboveSurface; } @@ -467,8 +467,8 @@ template class csDenseCellSet { if (isVoxel) { if (matMapPtr) { - auto material = matMapPtr->getMaterialAtIdx(materialId); - materialIds->at(cellIdx++) = static_cast(material); + int material = matMapPtr->getMaterialId(materialId); + materialIds->at(cellIdx++) = material; } else { materialIds->at(cellIdx++) = materialId; } @@ -665,7 +665,7 @@ template class csDenseCellSet { if (!materialMap) return; - auto numMaterials = materialMap->size(); + auto numMaterials = materialMap->getNumberOfLayers(); #pragma omp parallel for for (int i = 0; i < matIds->size(); i++) { @@ -673,10 +673,9 @@ template class csDenseCellSet { if (!cellSetAboveSurface) materialId--; if (materialId >= 0 && materialId < numMaterials) { - matIds->at(i) = - static_cast(materialMap->getMaterialAtIdx(materialId)); + matIds->at(i) = materialMap->getMaterialId(materialId); } else { - matIds->at(i) = static_cast(coverMaterial); + matIds->at(i) = coverMaterial; } } } diff --git a/include/viennaps/cellSet/csMeanFreePath.hpp b/include/viennaps/cellSet/csMeanFreePath.hpp index 954b36b..96b5f7a 100644 --- a/include/viennaps/cellSet/csMeanFreePath.hpp +++ b/include/viennaps/cellSet/csMeanFreePath.hpp @@ -35,8 +35,9 @@ template class csMeanFreePath { bulkLambda = passedBulkLambda; } - void setMaterial(const psMaterial passedMaterial) { - material = passedMaterial; + template + void setMaterial(const Material passedMaterial) { + material = static_cast(passedMaterial); } void setNumRaysPerCell(const int passedNumRaysPerCell) { @@ -99,7 +100,7 @@ template class csMeanFreePath { #pragma omp for for (int idx = 0; idx < numCells; ++idx) { - if (!psMaterialMap::isMaterial(materials->at(idx), material)) + if (static_cast(materials->at(idx)) != material) continue; auto cellCenter = cellSet->getCellCenter(idx); @@ -190,5 +191,5 @@ template class csMeanFreePath { NumericType bulkLambda = 0; NumericType maxLambda = 0.; long numRaysPerCell = 100; - psMaterial material = psMaterial::GAS; + int material = -1; }; diff --git a/include/viennaps/cellSet/csSegmentCells.hpp b/include/viennaps/cellSet/csSegmentCells.hpp index af626b0..3ac7ae7 100644 --- a/include/viennaps/cellSet/csSegmentCells.hpp +++ b/include/viennaps/cellSet/csSegmentCells.hpp @@ -5,18 +5,19 @@ template class csSegmentCells { lsSmartPointer> cellSet = nullptr; std::string cellTypeString = "CellType"; - psMaterial bulkMaterial = psMaterial::GAS; + int bulkMaterial = -1; public: csSegmentCells( const lsSmartPointer> &passedCellSet) : cellSet(passedCellSet) {} + template csSegmentCells( const lsSmartPointer> &passedCellSet, - const std::string cellTypeString, const psMaterial passedBulkMaterial) + const std::string cellTypeString, const Material passedBulkMaterial) : cellSet(passedCellSet), cellTypeString(cellTypeString), - bulkMaterial(passedBulkMaterial) {} + bulkMaterial(static_cast(passedBulkMaterial)) {} void setCellSet( const lsSmartPointer> &passedCellSet) { @@ -27,8 +28,9 @@ template class csSegmentCells { cellTypeString = passedCellTypeString; } - void setBulkMaterial(const psMaterial passedBulkMaterial) { - bulkMaterial = passedBulkMaterial; + template + void setBulkMaterial(const Material passedBulkMaterial) { + bulkMaterial = static_cast(passedBulkMaterial); } void apply() { @@ -38,11 +40,11 @@ template class csSegmentCells { #pragma omp parallel for for (int i = 0; i < materials->size(); ++i) { - if (!psMaterialMap::isMaterial(materials->at(i), bulkMaterial)) { + if (static_cast(materials->at(i)) != bulkMaterial) { auto neighbors = cellSet->getNeighbors(i); for (auto n : neighbors) { if (n >= 0 && - psMaterialMap::isMaterial(materials->at(n), bulkMaterial)) { + static_cast(materials->at(n)) == bulkMaterial) { cellType->at(i) = 0.; break; } diff --git a/include/viennaps/psDomain.hpp b/include/viennaps/psDomain.hpp index f5f1c2e..afb01b4 100644 --- a/include/viennaps/psDomain.hpp +++ b/include/viennaps/psDomain.hpp @@ -76,7 +76,8 @@ template class psDomain { } if (passedDomain->cellSet) { auto cellSetDepth = passedDomain->getCellSet()->getDepth(); - cellSet = csDomainType::New(levelSets, materialMap, cellSetDepth); + auto matMap = materialMap ? materialMap->getMaterialMap() : nullptr; + cellSet = csDomainType::New(levelSets, matMap, cellSetDepth); } else { cellSet = nullptr; } @@ -168,7 +169,8 @@ template class psDomain { cellSet = csDomainType::New(); cellSet->setCellSetPosition(isAboveSurface); cellSet->setCoverMaterial(coverMaterial); - cellSet->fromLevelSets(levelSets, materialMap, position); + auto matMap = materialMap ? materialMap->getMaterialMap() : nullptr; + cellSet->fromLevelSets(levelSets, matMap, position); } void setMaterialMap(materialMapType passedMaterialMap) { From 6a2473940ad71fd866cbef9f0aa6e347834201ca Mon Sep 17 00:00:00 2001 From: Felix S Date: Mon, 6 May 2024 10:54:12 +0200 Subject: [PATCH 3/6] Add missing coverMaterial default in cell set (Example oxideRegrowth works again) --- include/viennaps/cellSet/csDenseCellSet.hpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/include/viennaps/cellSet/csDenseCellSet.hpp b/include/viennaps/cellSet/csDenseCellSet.hpp index 1bec5bb..71c9437 100644 --- a/include/viennaps/cellSet/csDenseCellSet.hpp +++ b/include/viennaps/cellSet/csDenseCellSet.hpp @@ -467,7 +467,12 @@ template class csDenseCellSet { if (isVoxel) { if (matMapPtr) { - int material = matMapPtr->getMaterialId(materialId); + int index = materialId; + if (!cellSetAboveSurface) + --index; + int material = coverMaterial; + if (index >= 0 && index < matMapPtr->getNumberOfLayers()) + material = matMapPtr->getMaterialId(index); materialIds->at(cellIdx++) = material; } else { materialIds->at(cellIdx++) = materialId; From 174b2408e395f6cdb409293291dfd611f1394d2c Mon Sep 17 00:00:00 2001 From: Felix S Date: Tue, 7 May 2024 12:43:35 +0200 Subject: [PATCH 4/6] Replace psLogger with csLogger in cell set --- include/viennaps/cellSet/csDenseCellSet.hpp | 17 ++- include/viennaps/cellSet/csLogger.hpp | 148 ++++++++++++++++++++ include/viennaps/cellSet/csMeanFreePath.hpp | 5 +- include/viennaps/cellSet/csTracing.hpp | 2 +- include/viennaps/cellSet/csUtil.hpp | 22 +++ 5 files changed, 181 insertions(+), 13 deletions(-) create mode 100644 include/viennaps/cellSet/csLogger.hpp diff --git a/include/viennaps/cellSet/csDenseCellSet.hpp b/include/viennaps/cellSet/csDenseCellSet.hpp index 71c9437..e84f757 100644 --- a/include/viennaps/cellSet/csDenseCellSet.hpp +++ b/include/viennaps/cellSet/csDenseCellSet.hpp @@ -3,8 +3,7 @@ #include "csBVH.hpp" #include "csTracePath.hpp" #include "csUtil.hpp" - -#include "../psLogger.hpp" +#include "csLogger.hpp" #include #include @@ -339,7 +338,7 @@ template class csDenseCellSet { std::string line; if (!file.is_open()) { - psLogger::getInstance() + csLogger::getInstance() .addWarning("Could not open file " + fileName) .print(); return; @@ -347,7 +346,7 @@ template class csDenseCellSet { std::getline(file, line); if (std::stoi(line) != numberOfCells) { - psLogger::getInstance().addWarning("Incompatible cell set data.").print(); + csLogger::getInstance().addWarning("Incompatible cell set data.").print(); return; } @@ -554,7 +553,7 @@ template class csDenseCellSet { if (!cellNeighbors.empty() && !forceRebuild) return; - psUtils::Timer timer; + csUtil::Timer timer; timer.start(); const auto &cells = cellGrid->template getElements<(1 << D)>(); const auto &nodes = cellGrid->getNodes(); @@ -630,7 +629,7 @@ template class csDenseCellSet { } } timer.finish(); - psLogger::getInstance() + csLogger::getInstance() .addTiming("Building cell set neighborhood structure took", timer.currentDuration * 1e-9) .print(); @@ -717,7 +716,7 @@ template class csDenseCellSet { } void buildBVH() { - psUtils::Timer timer; + csUtil::Timer timer; timer.start(); auto &elems = cellGrid->template getElements<(1 << D)>(); auto &nodes = cellGrid->getNodes(); @@ -728,13 +727,13 @@ template class csDenseCellSet { auto &node = nodes[elems[elemIdx][n]]; auto cell = BVH->getCellIds(node); if (cell == nullptr) { - psLogger::getInstance().addError("BVH building error.").print(); + csLogger::getInstance().addError("BVH building error.").print(); } cell->insert(elemIdx); } } timer.finish(); - psLogger::getInstance() + csLogger::getInstance() .addTiming("Building cell set BVH took", timer.currentDuration * 1e-9) .print(); } diff --git a/include/viennaps/cellSet/csLogger.hpp b/include/viennaps/cellSet/csLogger.hpp new file mode 100644 index 0000000..36ef22c --- /dev/null +++ b/include/viennaps/cellSet/csLogger.hpp @@ -0,0 +1,148 @@ +#pragma once + +#include "csUtil.hpp" + +#include + +// verbosity levels: +// 0 errors +// 1 + warnings +// 2 + info +// 3 + timings +// 4 + intermediate output (meshes) +// 5 + debug +enum class csLogLevel : unsigned { + ERROR = 0, + WARNING = 1, + INFO = 2, + TIMING = 3, + INTERMEDIATE = 4, + DEBUG = 5 +}; + +/// Singleton class for thread-safe logging. The logger can be accessed via +/// csLogger::getInstance(). The logger can be configured to print messages of a +/// certain level or lower. The default level is INFO. The different logging +/// levels are: ERROR, WARNING, INFO, TIMING, INTERMEDIATE, DEBUG. The logger +/// can also be used to print timing information. +class csLogger { + std::string message; + + bool error = false; + const unsigned tabWidth = 4; + static csLogLevel logLevel; + + csLogger() {} + +public: + // delete constructors to result in better error messages by compilers + csLogger(const csLogger &) = delete; + void operator=(const csLogger &) = delete; + + // Set the log level for all instances of the logger. + static void setLogLevel(const csLogLevel passedLogLevel) { + logLevel = passedLogLevel; + } + + static unsigned getLogLevel() { return static_cast(logLevel); } + + static csLogger &getInstance() { + static csLogger instance; + return instance; + } + + // Add debug message if log level is high enough. + csLogger &addDebug(std::string s) { + if (getLogLevel() < 5) + return *this; +#pragma omp critical + { message += std::string(tabWidth, ' ') + "DEBUG: " + s + "\n"; } + return *this; + } + + // Add timing message if log level is high enough. + template + csLogger &addTiming(std::string s, csUtil::Timer &timer) { + if (getLogLevel() < 3) + return *this; +#pragma omp critical + { + message += std::string(tabWidth, ' ') + s + + " took: " + std::to_string(timer.currentDuration * 1e-9) + + " s \n"; + } + return *this; + } + + csLogger &addTiming(std::string s, double timeInSeconds) { + if (getLogLevel() < 3) + return *this; +#pragma omp critical + { + message += std::string(tabWidth, ' ') + s + ": " + + std::to_string(timeInSeconds) + " s \n"; + } + return *this; + } + + csLogger &addTiming(std::string s, double timeInSeconds, + double totalTimeInSeconds) { + if (getLogLevel() < 3) + return *this; +#pragma omp critical + { + message += std::string(tabWidth, ' ') + s + ": " + + std::to_string(timeInSeconds) + " s\n" + + std::string(tabWidth, ' ') + "Percent of total time: " + + std::to_string(timeInSeconds / totalTimeInSeconds * 100) + + "\n"; + } + return *this; + } + + // Add info message if log level is high enough. + csLogger &addInfo(std::string s) { + if (getLogLevel() < 2) + return *this; +#pragma omp critical + { message += std::string(tabWidth, ' ') + s + "\n"; } + return *this; + } + + // Add warning message if log level is high enough. + csLogger &addWarning(std::string s) { + if (getLogLevel() < 1) + return *this; +#pragma omp critical + { message += "\n" + std::string(tabWidth, ' ') + "WARNING: " + s + "\n"; } + return *this; + } + + // Add error message if log level is high enough. + csLogger &addError(std::string s, bool shouldAbort = true) { +#pragma omp critical + { + message += "\n" + std::string(tabWidth, ' ') + "ERROR: " + s + "\n"; + // always abort once error message should be printed + error = true; + } + // abort now if asked + if (shouldAbort) + print(); + return *this; + } + + // Print message to std::cout if log level is high enough. + void print(std::ostream &out = std::cout) { +#pragma omp critical + { + out << message; + message.clear(); + if (error) + abort(); + } + } +}; + +// initialize static member of logger +inline csLogLevel csLogger::logLevel = csLogLevel::INTERMEDIATE; diff --git a/include/viennaps/cellSet/csMeanFreePath.hpp b/include/viennaps/cellSet/csMeanFreePath.hpp index 96b5f7a..c9f0ab1 100644 --- a/include/viennaps/cellSet/csMeanFreePath.hpp +++ b/include/viennaps/cellSet/csMeanFreePath.hpp @@ -1,8 +1,7 @@ #pragma once #include "csDenseCellSet.hpp" - -#include "psLogger.hpp" +#include "csLogger.hpp" #include #include @@ -47,7 +46,7 @@ template class csMeanFreePath { NumericType getMaxLambda() const { return maxLambda; } void apply() { - psLogger::getInstance().addInfo("Calculating mean free path ...").print(); + csLogger::getInstance().addInfo("Calculating mean free path ...").print(); cellSet->addScalarData("MeanFreePath", 0.); runKernel(); } diff --git a/include/viennaps/cellSet/csTracing.hpp b/include/viennaps/cellSet/csTracing.hpp index 43e398d..ef30893 100644 --- a/include/viennaps/cellSet/csTracing.hpp +++ b/include/viennaps/cellSet/csTracing.hpp @@ -61,7 +61,7 @@ template class csTracing { std::array, 3> orthoBasis; if (usePrimaryDirection) { orthoBasis = rayInternal::getOrthonormalBasis(primaryDirection); - psLogger::getInstance() + csLogger::getInstance() .addInfo("Using primary direction: " + std::to_string(primaryDirection[0]) + " " + std::to_string(primaryDirection[1]) + " " + diff --git a/include/viennaps/cellSet/csUtil.hpp b/include/viennaps/cellSet/csUtil.hpp index da078a5..fa3899c 100644 --- a/include/viennaps/cellSet/csUtil.hpp +++ b/include/viennaps/cellSet/csUtil.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include #include @@ -21,6 +22,27 @@ template struct csVolumeParticle { }; namespace csUtil { + +template struct Timer { + using TimePoint = typename Clock::time_point; + + TimePoint mStart; + typename Clock::duration::rep totalDuration = 0.; // in ns + typename Clock::duration::rep currentDuration; // in ns + + void start() { mStart = Clock::now(); } + void finish() { + TimePoint end = Clock::now(); + typename Clock::duration dur(end - mStart); + currentDuration = dur.count(); + totalDuration += currentDuration; + } + void reset() { + currentDuration = 0.; + totalDuration = 0.; + } +}; + template void printTriple(const csTriple &p) { std::cout << "[" << p[0] << ", " << p[1] << ", " << p[2] << "]\n"; } From aa4cac2ca2842058e8c949c7bcbc2e607fdc44b0 Mon Sep 17 00:00:00 2001 From: Felix S Date: Mon, 27 May 2024 14:37:08 +0200 Subject: [PATCH 5/6] Move cell set to ViennaCS library --- .../atomicLayerDeposition.cpp | 4 +- include/viennaps/cellSet/csBVH.hpp | 36 - include/viennaps/cellSet/csBoundingVolume.hpp | 206 ----- include/viennaps/cellSet/csDenseCellSet.hpp | 779 ------------------ include/viennaps/cellSet/csLogger.hpp | 148 ---- include/viennaps/cellSet/csMeanFreePath.hpp | 194 ----- include/viennaps/cellSet/csPointSource.hpp | 41 - include/viennaps/cellSet/csSegmentCells.hpp | 57 -- include/viennaps/cellSet/csTracePath.hpp | 41 - include/viennaps/cellSet/csTracing.hpp | 260 ------ include/viennaps/cellSet/csTracingKernel.hpp | 297 ------- .../viennaps/cellSet/csTracingParticle.hpp | 52 -- include/viennaps/cellSet/csUtil.hpp | 183 ---- include/viennaps/models/psOxideRegrowth.hpp | 4 +- include/viennaps/models/psPlasmaDamage.hpp | 5 +- include/viennaps/psDomain.hpp | 4 +- 16 files changed, 9 insertions(+), 2302 deletions(-) delete mode 100644 include/viennaps/cellSet/csBVH.hpp delete mode 100644 include/viennaps/cellSet/csBoundingVolume.hpp delete mode 100644 include/viennaps/cellSet/csDenseCellSet.hpp delete mode 100644 include/viennaps/cellSet/csLogger.hpp delete mode 100644 include/viennaps/cellSet/csMeanFreePath.hpp delete mode 100644 include/viennaps/cellSet/csPointSource.hpp delete mode 100644 include/viennaps/cellSet/csSegmentCells.hpp delete mode 100644 include/viennaps/cellSet/csTracePath.hpp delete mode 100644 include/viennaps/cellSet/csTracing.hpp delete mode 100644 include/viennaps/cellSet/csTracingKernel.hpp delete mode 100644 include/viennaps/cellSet/csTracingParticle.hpp delete mode 100644 include/viennaps/cellSet/csUtil.hpp diff --git a/examples/atomicLayerDeposition/atomicLayerDeposition.cpp b/examples/atomicLayerDeposition/atomicLayerDeposition.cpp index 681a615..fc9d4de 100644 --- a/examples/atomicLayerDeposition/atomicLayerDeposition.cpp +++ b/examples/atomicLayerDeposition/atomicLayerDeposition.cpp @@ -1,5 +1,5 @@ -#include -#include +#include +#include #include #include diff --git a/include/viennaps/cellSet/csBVH.hpp b/include/viennaps/cellSet/csBVH.hpp deleted file mode 100644 index 835763b..0000000 --- a/include/viennaps/cellSet/csBVH.hpp +++ /dev/null @@ -1,36 +0,0 @@ -#pragma once - -#include "csBoundingVolume.hpp" -#include - -/// Helper class to quickly determine the cell in which a given point resides -/// in. To do so, an octree is built around the cell set structure. -template class csBVH { -private: - using BVPtrType = lsSmartPointer>; - using BoundsType = csPair>; - using CellIdsPtr = std::set *; - - unsigned numLayers = 1; - BVPtrType BV = nullptr; - -public: - csBVH(const BoundsType &domainBounds, unsigned layers = 1) - : numLayers(layers) { - BV = BVPtrType::New(domainBounds, numLayers - 1); - } - - BVPtrType getTopBV() { return BV; } - - void getLowestBVBounds(const std::array &point) { - BV->getBoundingVolumeBounds(point); - } - - CellIdsPtr getCellIds(const std::array &point) { - return BV->getCellIds(point); - } - - void clearCellIds() { BV->clear(); } - - size_t getTotalCellCount() { return BV->getTotalCellCounts(); } -}; diff --git a/include/viennaps/cellSet/csBoundingVolume.hpp b/include/viennaps/cellSet/csBoundingVolume.hpp deleted file mode 100644 index f3a9d6e..0000000 --- a/include/viennaps/cellSet/csBoundingVolume.hpp +++ /dev/null @@ -1,206 +0,0 @@ -#pragma once - -#include "csUtil.hpp" -#include - -#include -#include - -template class csBoundingVolume { -private: - using BVPtrType = lsSmartPointer>; - using BoundsType = csPair>; - using CellIdsPtr = std::set *; - - static constexpr int numCells = 1 << D; - std::array, numCells> cellIds; - std::array bounds; - std::array links; - int layer = -1; - -public: - csBoundingVolume() {} - csBoundingVolume(const BoundsType &outerBound, int thisLayer) - : layer(thisLayer) { - - if constexpr (D == 3) - buildBounds3D(outerBound); - else - buildBounds2D(outerBound); - - if (layer > 0) { - for (size_t i = 0; i < numCells; i++) { - links[i] = BVPtrType::New(bounds[i], layer - 1); - } - } - } - - CellIdsPtr getCellIds(const std::array &point) { - auto vid = getVolumeIndex(point); - // assert(vid < numCells && "Point in invalid BV"); - - if (vid == numCells) - return nullptr; - - if (layer == 0) { - return &cellIds[vid]; - } - - return links[vid]->getCellIds(point); - } - - void getBoundingVolumeBounds(const std::array &point) { - auto vid = getVolumeIndex(point); - assert(vid < numCells && "Point in invalid BV"); - - if (layer == 0) { - printBound(vid); - return; - } - - links[vid]->getBoundingVolumeBounds(point); - } - - size_t getTotalCellCounts() { - size_t count = 0; - if (layer == 0) { - for (size_t i = 0; i < numCells; i++) { - count += cellIds[i].size(); - } - return count; - } else { - for (size_t i = 0; i < numCells; i++) { - count += links[i]->getTotalCellCounts(); - } - } - return count; - } - - void clear() { - if (layer == 0) { - for (size_t i = 0; i < numCells; i++) { - cellIds[i].clear(); - } - } else { - for (size_t i = 0; i < numCells; i++) { - links[i]->clear(); - } - } - } - - BVPtrType getLink(const std::array &point) { - auto vid = getVolumeIndex(point); - return getLink(vid); - } - - BVPtrType getLink(size_t vid) { return links[vid]; } - - size_t getVolumeIndex(const std::array &point) { - size_t vid = numCells; - for (size_t idx = 0; idx < numCells; idx++) { - if (insideVolume(point, idx)) { - vid = idx; - break; - } - } - return vid; - } - - bool insideVolume(const std::array &p, const size_t idx) const { - if constexpr (D == 3) { - return p[0] > bounds[idx][0][0] && p[0] <= bounds[idx][1][0] && - p[1] > bounds[idx][0][1] && p[1] <= bounds[idx][1][1] && - p[2] > bounds[idx][0][2] && p[2] <= bounds[idx][1][2]; - } else { - return p[0] > bounds[idx][0][0] && p[0] <= bounds[idx][1][0] && - p[1] > bounds[idx][0][1] && p[1] <= bounds[idx][1][1]; - } - } - - void printBound(size_t vid) { - std::cout << "Bounding volume span: ["; - for (int i = 0; i < D; i++) - std::cout << bounds[vid][0][i] << ", "; - std::cout << "] - ["; - for (int i = 0; i < D; i++) - std::cout << bounds[vid][1][i] << ", "; - std::cout << "]\n"; - } - -private: - void buildBounds2D(const BoundsType &outerBound) { - auto xExt = (outerBound[1][0] - outerBound[0][0]) / 2.; - auto yExt = (outerBound[1][1] - outerBound[0][1]) / 2.; - - const auto BVH1 = - csPair>{outerBound[0][0], outerBound[0][1], - outerBound[0][0] + xExt, outerBound[0][1] + yExt}; - const auto BVH2 = - csPair>{outerBound[0][0] + xExt, outerBound[0][1], - outerBound[0][0] + 2 * xExt, outerBound[0][1] + yExt}; - const auto BVH3 = csPair>{ - outerBound[0][0] + xExt, - outerBound[0][1] + yExt, - outerBound[0][0] + 2 * xExt, - outerBound[0][1] + 2 * yExt, - }; - const auto BVH4 = csPair>{ - outerBound[0][0], - outerBound[0][1] + yExt, - outerBound[0][0] + xExt, - outerBound[0][1] + 2 * yExt, - }; - - bounds = - std::array>, numCells>{BVH1, BVH2, BVH3, BVH4}; - } - - void buildBounds3D(const BoundsType &outerBound) { - auto xExt = (outerBound[1][0] - outerBound[0][0]) / T(2); - auto yExt = (outerBound[1][1] - outerBound[0][1]) / T(2); - auto zExt = (outerBound[1][2] - outerBound[0][2]) / T(2); - - const auto BVH1 = - csPair>{outerBound[0][0], outerBound[0][1], - outerBound[0][2], outerBound[0][0] + xExt, - outerBound[0][1] + yExt, outerBound[0][2] + zExt}; - const auto BVH2 = csPair>{ - outerBound[0][0] + xExt, outerBound[0][1], - outerBound[0][2], outerBound[0][0] + 2 * xExt, - outerBound[0][1] + yExt, outerBound[0][2] + zExt}; - const auto BVH3 = csPair>{outerBound[0][0] + xExt, - outerBound[0][1] + yExt, - outerBound[0][2], - outerBound[0][0] + 2 * xExt, - outerBound[0][1] + 2 * yExt, - outerBound[0][2] + zExt}; - const auto BVH4 = csPair>{outerBound[0][0], - outerBound[0][1] + yExt, - outerBound[0][2], - outerBound[0][0] + xExt, - outerBound[0][1] + 2 * yExt, - outerBound[0][2] + zExt}; - - // top - const auto BVH5 = csPair>{ - outerBound[0][0], outerBound[0][1], - outerBound[0][2] + zExt, outerBound[0][0] + xExt, - outerBound[0][1] + yExt, outerBound[0][2] + 2 * zExt}; - const auto BVH6 = csPair>{ - outerBound[0][0] + xExt, outerBound[0][1], - outerBound[0][2] + zExt, outerBound[0][0] + 2 * xExt, - outerBound[0][1] + yExt, outerBound[0][2] + 2 * zExt}; - const auto BVH7 = csPair>{ - outerBound[0][0] + xExt, outerBound[0][1] + yExt, - outerBound[0][2] + zExt, outerBound[0][0] + 2 * xExt, - outerBound[0][1] + 2 * yExt, outerBound[0][2] + 2 * zExt}; - const auto BVH8 = csPair>{outerBound[0][0], - outerBound[0][1] + yExt, - outerBound[0][2] + zExt, - outerBound[0][0] + xExt, - outerBound[0][1] + 2 * yExt, - outerBound[0][2] + 2 * zExt}; - bounds = std::array>, numCells>{ - BVH1, BVH2, BVH3, BVH4, BVH5, BVH6, BVH7, BVH8}; - } -}; diff --git a/include/viennaps/cellSet/csDenseCellSet.hpp b/include/viennaps/cellSet/csDenseCellSet.hpp deleted file mode 100644 index e84f757..0000000 --- a/include/viennaps/cellSet/csDenseCellSet.hpp +++ /dev/null @@ -1,779 +0,0 @@ -#pragma once - -#include "csBVH.hpp" -#include "csTracePath.hpp" -#include "csUtil.hpp" -#include "csLogger.hpp" - -#include -#include -#include -#include -#include -#include -#include - -#include - -/** - This class represents a cell-based voxel implementation of a volume. The - depth of the cell set in z-direction can be specified. -*/ -template class csDenseCellSet { -private: - using gridType = lsSmartPointer>; - using levelSetsType = - lsSmartPointer>>>; - using materialMapType = lsSmartPointer; - - levelSetsType levelSets = nullptr; - gridType cellGrid = nullptr; - lsSmartPointer> surface = nullptr; - lsSmartPointer> BVH = nullptr; - materialMapType materialMap = nullptr; - - T gridDelta; - T depth = 0.; - std::size_t numberOfCells; - int BVHlayers = 0; - - std::vector> cellNeighbors; // -x, x, -y, y, -z, z - hrleVectorType minIndex, maxIndex; - - bool cellSetAboveSurface = false; - int coverMaterial = -1; - std::bitset periodicBoundary; - - std::vector *fillingFractions; - const T eps = 1e-4; - -public: - csDenseCellSet() {} - - csDenseCellSet(levelSetsType passedLevelSets, - materialMapType passedMaterialMap = nullptr, - T passedDepth = 0., bool passedCellSetPosition = false) - : levelSets(passedLevelSets), cellSetAboveSurface(passedCellSetPosition) { - fromLevelSets(passedLevelSets, passedMaterialMap, passedDepth); - } - - void fromLevelSets(levelSetsType passedLevelSets, - materialMapType passedMaterialMap = nullptr, - T passedDepth = 0.) { - levelSets = passedLevelSets; - materialMap = passedMaterialMap; - - if (cellGrid == nullptr) - cellGrid = lsSmartPointer>::New(); - - if (surface == nullptr) - surface = lsSmartPointer>::New(levelSets->back()); - else - surface->deepCopy(levelSets->back()); - - gridDelta = surface->getGrid().getGridDelta(); - - depth = passedDepth; - std::vector>> levelSetsInOrder; - auto plane = lsSmartPointer>::New(surface->getGrid()); - { - T origin[D] = {0.}; - T normal[D] = {0.}; - origin[D - 1] = depth; - normal[D - 1] = 1.; - lsMakeGeometry(plane, - lsSmartPointer>::New(origin, normal)) - .apply(); - } - if (!cellSetAboveSurface) - levelSetsInOrder.push_back(plane); - for (auto ls : *levelSets) - levelSetsInOrder.push_back(ls); - if (cellSetAboveSurface) { - levelSetsInOrder.push_back(plane); - } - - calculateMinMaxIndex(levelSetsInOrder); - lsToVoxelMesh(levelSetsInOrder, cellGrid).apply(); - // lsToVoxelMesh also saves the extent in the cell grid - -#ifndef NDEBUG - int db_ls = 0; - for (auto &ls : levelSetsInOrder) { - auto mesh = lsSmartPointer>::New(); - lsToSurfaceMesh(ls, mesh).apply(); - lsVTKWriter(mesh, "cellSet_debug_" + std::to_string(db_ls++) + ".vtp") - .apply(); - } - lsVTKWriter(cellGrid, "cellSet_debug_init.vtu").apply(); -#endif - - adjustMaterialIds(); - - // create filling fractions as default scalar cell data - numberOfCells = cellGrid->template getElements<(1 << D)>().size(); - std::vector fillingFractionsTemp(numberOfCells, 0.); - - cellGrid->getCellData().insertNextScalarData( - std::move(fillingFractionsTemp), "fillingFraction"); - fillingFractions = cellGrid->getCellData().getScalarData("fillingFraction"); - - // calculate number of BVH layers - for (unsigned i = 0; i < D; ++i) { - cellGrid->minimumExtent[i] -= eps; - cellGrid->maximumExtent[i] += eps; - } - auto minExtent = cellGrid->maximumExtent[0] - cellGrid->minimumExtent[0]; - minExtent = std::min(minExtent, cellGrid->maximumExtent[1] - - cellGrid->minimumExtent[1]); - if constexpr (D == 3) - minExtent = std::min(minExtent, cellGrid->maximumExtent[2] - - cellGrid->minimumExtent[2]); - - BVHlayers = 0; - while (minExtent / 2 > gridDelta) { - BVHlayers++; - minExtent /= 2; - } - - BVH = lsSmartPointer>::New(getBoundingBox(), BVHlayers); - buildBVH(); - } - - csPair> getBoundingBox() const { - if constexpr (D == 3) - return csPair>{cellGrid->minimumExtent, - cellGrid->maximumExtent}; - else - return csPair>{ - cellGrid->minimumExtent[0], cellGrid->minimumExtent[1], - cellGrid->maximumExtent[0], cellGrid->maximumExtent[1]}; - } - - void setPeriodicBoundary(std::array isPeriodic) { - for (int i = 0; i < D; i++) { - periodicBoundary[i] = isPeriodic[i]; - } - } - - std::vector *addScalarData(std::string name, T initValue = 0.) { - if (cellGrid->getCellData().getScalarData(name) != nullptr) { - auto data = cellGrid->getCellData().getScalarData(name); - data->resize(numberOfCells, initValue); - std::fill(data->begin(), data->end(), initValue); - return data; - } - std::vector newData(numberOfCells, initValue); - cellGrid->getCellData().insertNextScalarData(std::move(newData), name); - fillingFractions = cellGrid->getCellData().getScalarData("fillingFraction"); - return cellGrid->getCellData().getScalarData(name); - } - - T getDepth() const { return depth; } - - T getGridDelta() const { return gridDelta; } - - std::vector> &getNodes() const { - return cellGrid->getNodes(); - } - - const std::array &getNode(unsigned int idx) const { - return cellGrid->getNodes()[idx]; - } - - std::vector> &getElements() const { - return cellGrid->template getElements<(1 << D)>(); - } - - const std::array &getElement(unsigned int idx) const { - return cellGrid->template getElements<(1 << D)>()[idx]; - } - - lsSmartPointer> getSurface() { return surface; } - - lsSmartPointer> getCellGrid() { return cellGrid; } - - levelSetsType getLevelSets() const { return levelSets; } - - size_t getNumberOfCells() const { return numberOfCells; } - - std::vector *getFillingFractions() const { return fillingFractions; } - - T getFillingFraction(const std::array &point) { - csTriple point3 = {0., 0., 0.}; - for (int i = 0; i < D; i++) - point3[i] = point[i]; - auto idx = findIndex(point3); - if (idx < 0) - return -1.; - - return getFillingFractions()->at(idx); - } - - T getAverageFillingFraction(const std::array &point, - const T radius) const { - T sum = 0.; - int count = 0; - for (int i = 0; i < numberOfCells; i++) { - auto &cell = cellGrid->template getElements<(1 << D)>()[i]; - auto node = cellGrid->getNodes()[cell[0]]; - for (int j = 0; j < D; j++) - node[j] += gridDelta / 2.; - if (csUtil::distance(node, point) < radius) { - sum += fillingFractions->at(i); - count++; - } - } - return sum / count; - } - - std::array getCellCenter(unsigned long idx) const { - auto center = - cellGrid - ->getNodes()[cellGrid->template getElements<(1 << D)>()[idx][0]]; - for (int i = 0; i < D; i++) - center[i] += gridDelta / 2.; - return center; - } - - int getIndex(const std::array &point) const { return findIndex(point); } - - std::vector *getScalarData(std::string name) { - return cellGrid->getCellData().getScalarData(name); - } - - std::vector getScalarDataLabels() const { - std::vector labels; - auto numScalarData = cellGrid->getCellData().getScalarDataSize(); - for (int i = 0; i < numScalarData; i++) { - labels.push_back(cellGrid->getCellData().getScalarDataLabel(i)); - } - return labels; - } - - // Set whether the cell set should be created below (false) or above (true) - // the surface. - void setCellSetPosition(const bool passedCellSetPosition) { - cellSetAboveSurface = passedCellSetPosition; - } - - template - void setCoverMaterial(const Material passedCoverMaterial) { - coverMaterial = static_cast(passedCoverMaterial); - } - - bool getCellSetPosition() const { return cellSetAboveSurface; } - - // Sets the filling fraction at given cell index. - bool setFillingFraction(const int idx, const T fill) { - if (idx < 0) - return false; - - getFillingFractions()->at(idx) = fill; - return true; - } - - // Sets the filling fraction for cell which contains given point. - bool setFillingFraction(const std::array &point, const T fill) { - auto idx = findIndex(point); - return setFillingFraction(idx, fill); - } - - // Add to the filling fraction at given cell index. - bool addFillingFraction(const int idx, const T fill) { - if (idx < 0) - return false; - - fillingFractions->at(idx) += fill; - return true; - } - - // Add to the filling fraction for cell which contains given point. - bool addFillingFraction(const std::array &point, T fill) { - auto idx = findIndex(point); - return addFillingFraction(idx, fill); - } - - // Add to the filling fraction for cell which contains given point only if the - // cell has the specified material ID. - bool addFillingFractionInMaterial(const std::array &point, T fill, - int materialId) { - auto idx = findIndex(point); - if (getScalarData("Material")->at(idx) == materialId) - return addFillingFraction(idx, fill); - else - return false; - } - - // Write the cell set as .vtu file - void writeVTU(std::string fileName) { - lsVTKWriter(cellGrid, fileName).apply(); - } - - // Save cell set data in simple text format - void writeCellSetData(std::string fileName) const { - auto numScalarData = cellGrid->getCellData().getScalarDataSize(); - - std::ofstream file(fileName); - file << numberOfCells << "\n"; - for (int i = 0; i < numScalarData; i++) { - auto label = cellGrid->getCellData().getScalarDataLabel(i); - file << label << ","; - } - file << "\n"; - - for (size_t j = 0; j < numberOfCells; j++) { - for (int i = 0; i < numScalarData; i++) { - file << cellGrid->getCellData().getScalarData(i)->at(j) << ","; - } - file << "\n"; - } - - file.close(); - } - - // Read cell set data from text - void readCellSetData(std::string fileName) { - std::ifstream file(fileName); - std::string line; - - if (!file.is_open()) { - csLogger::getInstance() - .addWarning("Could not open file " + fileName) - .print(); - return; - } - - std::getline(file, line); - if (std::stoi(line) != numberOfCells) { - csLogger::getInstance().addWarning("Incompatible cell set data.").print(); - return; - } - - std::vector labels; - std::getline(file, line); - { - std::stringstream ss(line); - std::string label; - while (std::getline(ss, label, ',')) { - labels.push_back(label); - } - } - - std::vector *> cellDataP; - for (int i = 0; i < labels.size(); i++) { - auto dataP = getScalarData(labels[i]); - if (dataP == nullptr) { - dataP = addScalarData(labels[i], 0.); - } - } - - for (int i = 0; i < labels.size(); i++) { - cellDataP.push_back(getScalarData(labels[i])); - } - - std::size_t j = 0; - while (std::getline(file, line)) { - std::stringstream ss(line); - std::size_t i = 0; - std::string value; - while (std::getline(ss, value, ',')) - cellDataP[i++]->at(j) = std::stod(value); - - j++; - } - assert(j == numberOfCells && "Data incompatible"); - - file.close(); - } - - // Clear the filling fractions - void clear() { - auto ff = getFillingFractions(); - std::fill(ff->begin(), ff->end(), 0.); - } - - // Update the material IDs of the cell set. This function should be called if - // the level sets, the cell set is made out of, have changed. This does not - // work if the surface of the volume has changed. In this case, call the - // function "updateSurface" first. - void updateMaterials() { - auto materialIds = getScalarData("Material"); - - // create overlay material - std::vector>> levelSetsInOrder; - auto plane = lsSmartPointer>::New(surface->getGrid()); - { - T origin[D] = {0.}; - T normal[D] = {0.}; - origin[D - 1] = depth; - normal[D - 1] = 1.; - lsMakeGeometry(plane, - lsSmartPointer>::New(origin, normal)) - .apply(); - } - if (!cellSetAboveSurface) - levelSetsInOrder.push_back(plane); - for (auto ls : *levelSets) - levelSetsInOrder.push_back(ls); - if (cellSetAboveSurface) - levelSetsInOrder.push_back(plane); - - // set up iterators for all materials - std::vector::DomainType>> - iterators; - for (auto it = levelSetsInOrder.begin(); it != levelSetsInOrder.end(); - ++it) { - iterators.push_back( - hrleConstDenseCellIterator::DomainType>( - (*it)->getDomain(), minIndex)); - } - - // move iterator for lowest material id and then adjust others if they are - // needed - const materialMapType matMapPtr = materialMap; - unsigned cellIdx = 0; - for (; iterators.front().getIndices() < maxIndex; - iterators.front().next()) { - // go over all materials - for (unsigned materialId = 0; materialId < levelSetsInOrder.size(); - ++materialId) { - - auto &cellIt = iterators[materialId]; - cellIt.goToIndicesSequential(iterators.front().getIndices()); - - // find out whether the centre of the box is inside - T centerValue = 0.; - for (int i = 0; i < (1 << D); ++i) { - centerValue += cellIt.getCorner(i).getValue(); - } - - if (centerValue <= 0.) { - bool isVoxel; - // check if voxel is in bounds - for (unsigned i = 0; i < (1 << D); ++i) { - hrleVectorType index; - isVoxel = true; - for (unsigned j = 0; j < D; ++j) { - index[j] = - cellIt.getIndices(j) + cellIt.getCorner(i).getOffset()[j]; - if (index[j] > maxIndex[j]) { - isVoxel = false; - break; - } - } - } - - if (isVoxel) { - if (matMapPtr) { - int index = materialId; - if (!cellSetAboveSurface) - --index; - int material = coverMaterial; - if (index >= 0 && index < matMapPtr->getNumberOfLayers()) - material = matMapPtr->getMaterialId(index); - materialIds->at(cellIdx++) = material; - } else { - materialIds->at(cellIdx++) = materialId; - } - } - - // jump out of material for loop - break; - } - } - } - assert(cellIdx == numberOfCells && - "Cell set changed in `updateMaterials()'"); - } - - // Updates the surface of the cell set. The new surface should be below the - // old surface as this function can only remove cells from the cell set. - void updateSurface() { - auto updateCellGrid = lsSmartPointer>::New(); - - lsToVoxelMesh voxelConverter(updateCellGrid); - { - auto plane = lsSmartPointer>::New(surface->getGrid()); - T origin[D] = {0.}; - T normal[D] = {0.}; - origin[D - 1] = depth; - normal[D - 1] = 1.; - - lsMakeGeometry(plane, - lsSmartPointer>::New(origin, normal)) - .apply(); - voxelConverter.insertNextLevelSet(plane); - } - voxelConverter.insertNextLevelSet(levelSets->back()); - voxelConverter.insertNextLevelSet(surface); - voxelConverter.apply(); - - auto cutMatIds = updateCellGrid->getCellData().getScalarData("Material"); - auto &elements = cellGrid->template getElements<(1 << D)>(); - - const auto nCutCells = - updateCellGrid->template getElements<(1 << D)>().size(); - - auto numScalarData = cellGrid->getCellData().getScalarDataSize(); - - for (int elIdx = nCutCells - 1; elIdx >= 0; elIdx--) { - if (cutMatIds->at(elIdx) == 2) { - for (int i = 0; i < numScalarData; i++) { - auto data = cellGrid->getCellData().getScalarData(i); - data->erase(data->begin() + elIdx); - } - elements.erase(elements.begin() + elIdx); - } - } - numberOfCells = elements.size(); - surface->deepCopy(levelSets->back()); - - buildBVH(); - } - - // Merge a trace path to the cell set. - void mergePath(csTracePath &path, T factor = 1.) { - auto ff = getFillingFractions(); - if (!path.getData().empty()) { - for (const auto it : path.getData()) { - ff->at(it.first) += it.second / factor; - } - } - - if (!path.getGridData().empty()) { - const auto &data = path.getGridData(); - for (size_t idx = 0; idx < numberOfCells; idx++) { - ff->at(idx) += data[idx] / factor; - } - } - } - - void buildNeighborhood(bool forceRebuild = false) { - if (!cellNeighbors.empty() && !forceRebuild) - return; - - csUtil::Timer timer; - timer.start(); - const auto &cells = cellGrid->template getElements<(1 << D)>(); - const auto &nodes = cellGrid->getNodes(); - unsigned const numNodes = nodes.size(); - unsigned const numCells = cells.size(); - cellNeighbors.resize(numCells); - const bool usePeriodicBoundary = periodicBoundary.any(); - - std::vector> nodeCellConnections(numNodes); - - // for each node, store which cells are connected with the node - for (unsigned cellIdx = 0; cellIdx < numCells; cellIdx++) { - for (unsigned cellNodeIdx = 0; cellNodeIdx < (1 << D); cellNodeIdx++) { - nodeCellConnections[cells[cellIdx][cellNodeIdx]].push_back(cellIdx); - } - } - -#pragma omp parallel for - for (int cellIdx = 0; cellIdx < numCells; cellIdx++) { - auto coord = nodes[cells[cellIdx][0]]; - for (int i = 0; i < D; i++) { - coord[i] += gridDelta / 2.; - cellNeighbors[cellIdx][i] = -1; - cellNeighbors[cellIdx][i + D] = -1; - } - - if (usePeriodicBoundary) { - auto onBoundary = isBoundaryCell(coord); - if (onBoundary.any()) { - /*look for neighbor cells using BVH*/ - for (std::size_t i = 0; i < 2 * D; i++) { - auto neighborCoord = coord; - bool minBoundary = i % 2 == 0; - - if (onBoundary.test(i)) { - // wrap around boundary - if (!minBoundary) { - neighborCoord[i / 2] = - cellGrid->minimumExtent[i / 2] + gridDelta / 2.; - } else { - neighborCoord[i / 2] = - cellGrid->maximumExtent[i / 2] - gridDelta / 2.; - } - } else { - neighborCoord[i / 2] += minBoundary ? -gridDelta : gridDelta; - } - cellNeighbors[cellIdx][i] = findIndex(neighborCoord); - } - continue; - } - } - - for (unsigned cellNodeIdx = 0; cellNodeIdx < (1 << D); cellNodeIdx++) { - auto &cellsAtNode = nodeCellConnections[cells[cellIdx][cellNodeIdx]]; - - for (const auto &neighborCell : cellsAtNode) { - if (neighborCell != cellIdx) { - - auto neighborCoord = getCellCenter(neighborCell); - - if (csUtil::distance(coord, neighborCoord) < gridDelta + eps) { - - for (int i = 0; i < D; i++) { - if (coord[i] - neighborCoord[i] > gridDelta / 2.) { - cellNeighbors[cellIdx][i * 2] = neighborCell; - } else if (coord[i] - neighborCoord[i] < -gridDelta / 2.) { - cellNeighbors[cellIdx][i * 2 + 1] = neighborCell; - } - } - } - } - } - } - } - timer.finish(); - csLogger::getInstance() - .addTiming("Building cell set neighborhood structure took", - timer.currentDuration * 1e-9) - .print(); - } - - const std::array &getNeighbors(unsigned long cellIdx) const { - assert(cellIdx < numberOfCells && "Cell idx out of bounds"); - return cellNeighbors[cellIdx]; - } - - bool isPointInCell(const csTriple &point, unsigned int cellIdx) const { - const auto &elem = getElement(cellIdx); - const auto &cellMin = getNode(elem[0]); - return isPointInCell(point, cellMin); - } - -private: - int findIndex(const csTriple &point) const { - const auto &elems = cellGrid->template getElements<(1 << D)>(); - const auto &nodes = cellGrid->getNodes(); - int idx = -1; - - auto cellIds = BVH->getCellIds(point); - if (!cellIds) - return idx; - for (const auto cellId : *cellIds) { - if (isPointInCell(point, nodes[elems[cellId][0]])) { - idx = cellId; - break; - } - } - return idx; - } - - void adjustMaterialIds() { - auto matIds = getScalarData("Material"); - if (!materialMap) - return; - - auto numMaterials = materialMap->getNumberOfLayers(); - -#pragma omp parallel for - for (int i = 0; i < matIds->size(); i++) { - int materialId = static_cast(matIds->at(i)); - if (!cellSetAboveSurface) - materialId--; - if (materialId >= 0 && materialId < numMaterials) { - matIds->at(i) = materialMap->getMaterialId(materialId); - } else { - matIds->at(i) = coverMaterial; - } - } - } - - int findSurfaceHitPoint(csTriple &hitPoint, const csTriple &direction) { - // find surface hit point - auto idx = findIndex(hitPoint); - - if (idx > 0) - return idx; - - auto moveDirection = multNew(direction, gridDelta / 2.); - size_t sanityCounter = 0; - while (idx < 0) { - add(hitPoint, moveDirection); - if (++sanityCounter > 100 || !checkBoundsPeriodic(hitPoint)) { - return -1; - } - idx = findIndex(hitPoint); - } - - return idx; - } - - bool isPointInCell(const csTriple &point, - const csTriple &cellMin) const { - if constexpr (D == 3) - return point[0] >= cellMin[0] && point[0] <= (cellMin[0] + gridDelta) && - point[1] >= cellMin[1] && point[1] <= (cellMin[1] + gridDelta) && - point[2] >= cellMin[2] && point[2] <= (cellMin[2] + gridDelta); - else - return point[0] >= cellMin[0] && point[0] <= (cellMin[0] + gridDelta) && - point[1] >= cellMin[1] && point[1] <= (cellMin[1] + gridDelta); - } - - void buildBVH() { - csUtil::Timer timer; - timer.start(); - auto &elems = cellGrid->template getElements<(1 << D)>(); - auto &nodes = cellGrid->getNodes(); - BVH->clearCellIds(); - - for (size_t elemIdx = 0; elemIdx < elems.size(); elemIdx++) { - for (size_t n = 0; n < (1 << D); n++) { - auto &node = nodes[elems[elemIdx][n]]; - auto cell = BVH->getCellIds(node); - if (cell == nullptr) { - csLogger::getInstance().addError("BVH building error.").print(); - } - cell->insert(elemIdx); - } - } - timer.finish(); - csLogger::getInstance() - .addTiming("Building cell set BVH took", timer.currentDuration * 1e-9) - .print(); - } - - void calculateMinMaxIndex( - const std::vector>> &levelSetsInOrder) { - // set to zero - for (unsigned i = 0; i < D; ++i) { - minIndex[i] = std::numeric_limits::max(); - maxIndex[i] = std::numeric_limits::lowest(); - } - for (unsigned l = 0; l < levelSetsInOrder.size(); ++l) { - auto &grid = levelSetsInOrder[l]->getGrid(); - auto &domain = levelSetsInOrder[l]->getDomain(); - for (unsigned i = 0; i < D; ++i) { - minIndex[i] = std::min(minIndex[i], (grid.isNegBoundaryInfinite(i)) - ? domain.getMinRunBreak(i) - : grid.getMinBounds(i)); - - maxIndex[i] = std::max(maxIndex[i], (grid.isPosBoundaryInfinite(i)) - ? domain.getMaxRunBreak(i) - : grid.getMaxBounds(i)); - } - } - } - - std::bitset<2 * D> isBoundaryCell(const std::array &cellCoord) { - std::bitset<2 * D> onBoundary; - for (int i = 0; i < 2 * D; i += 2) { - if (!periodicBoundary[i / 2]) - continue; - if (cellCoord[i / 2] - cellGrid->minimumExtent[i / 2] < gridDelta) { - /* cell is at min boundary */ - onBoundary.set(i); - } - if (cellGrid->maximumExtent[i / 2] - cellCoord[i / 2] < gridDelta) { - /* cell is at max boundary */ - onBoundary.set(i + 1); - } - } - return onBoundary; - } -}; diff --git a/include/viennaps/cellSet/csLogger.hpp b/include/viennaps/cellSet/csLogger.hpp deleted file mode 100644 index 36ef22c..0000000 --- a/include/viennaps/cellSet/csLogger.hpp +++ /dev/null @@ -1,148 +0,0 @@ -#pragma once - -#include "csUtil.hpp" - -#include - -// verbosity levels: -// 0 errors -// 1 + warnings -// 2 + info -// 3 + timings -// 4 + intermediate output (meshes) -// 5 + debug -enum class csLogLevel : unsigned { - ERROR = 0, - WARNING = 1, - INFO = 2, - TIMING = 3, - INTERMEDIATE = 4, - DEBUG = 5 -}; - -/// Singleton class for thread-safe logging. The logger can be accessed via -/// csLogger::getInstance(). The logger can be configured to print messages of a -/// certain level or lower. The default level is INFO. The different logging -/// levels are: ERROR, WARNING, INFO, TIMING, INTERMEDIATE, DEBUG. The logger -/// can also be used to print timing information. -class csLogger { - std::string message; - - bool error = false; - const unsigned tabWidth = 4; - static csLogLevel logLevel; - - csLogger() {} - -public: - // delete constructors to result in better error messages by compilers - csLogger(const csLogger &) = delete; - void operator=(const csLogger &) = delete; - - // Set the log level for all instances of the logger. - static void setLogLevel(const csLogLevel passedLogLevel) { - logLevel = passedLogLevel; - } - - static unsigned getLogLevel() { return static_cast(logLevel); } - - static csLogger &getInstance() { - static csLogger instance; - return instance; - } - - // Add debug message if log level is high enough. - csLogger &addDebug(std::string s) { - if (getLogLevel() < 5) - return *this; -#pragma omp critical - { message += std::string(tabWidth, ' ') + "DEBUG: " + s + "\n"; } - return *this; - } - - // Add timing message if log level is high enough. - template - csLogger &addTiming(std::string s, csUtil::Timer &timer) { - if (getLogLevel() < 3) - return *this; -#pragma omp critical - { - message += std::string(tabWidth, ' ') + s + - " took: " + std::to_string(timer.currentDuration * 1e-9) + - " s \n"; - } - return *this; - } - - csLogger &addTiming(std::string s, double timeInSeconds) { - if (getLogLevel() < 3) - return *this; -#pragma omp critical - { - message += std::string(tabWidth, ' ') + s + ": " + - std::to_string(timeInSeconds) + " s \n"; - } - return *this; - } - - csLogger &addTiming(std::string s, double timeInSeconds, - double totalTimeInSeconds) { - if (getLogLevel() < 3) - return *this; -#pragma omp critical - { - message += std::string(tabWidth, ' ') + s + ": " + - std::to_string(timeInSeconds) + " s\n" + - std::string(tabWidth, ' ') + "Percent of total time: " + - std::to_string(timeInSeconds / totalTimeInSeconds * 100) + - "\n"; - } - return *this; - } - - // Add info message if log level is high enough. - csLogger &addInfo(std::string s) { - if (getLogLevel() < 2) - return *this; -#pragma omp critical - { message += std::string(tabWidth, ' ') + s + "\n"; } - return *this; - } - - // Add warning message if log level is high enough. - csLogger &addWarning(std::string s) { - if (getLogLevel() < 1) - return *this; -#pragma omp critical - { message += "\n" + std::string(tabWidth, ' ') + "WARNING: " + s + "\n"; } - return *this; - } - - // Add error message if log level is high enough. - csLogger &addError(std::string s, bool shouldAbort = true) { -#pragma omp critical - { - message += "\n" + std::string(tabWidth, ' ') + "ERROR: " + s + "\n"; - // always abort once error message should be printed - error = true; - } - // abort now if asked - if (shouldAbort) - print(); - return *this; - } - - // Print message to std::cout if log level is high enough. - void print(std::ostream &out = std::cout) { -#pragma omp critical - { - out << message; - message.clear(); - if (error) - abort(); - } - } -}; - -// initialize static member of logger -inline csLogLevel csLogger::logLevel = csLogLevel::INTERMEDIATE; diff --git a/include/viennaps/cellSet/csMeanFreePath.hpp b/include/viennaps/cellSet/csMeanFreePath.hpp deleted file mode 100644 index c9f0ab1..0000000 --- a/include/viennaps/cellSet/csMeanFreePath.hpp +++ /dev/null @@ -1,194 +0,0 @@ -#pragma once - -#include "csDenseCellSet.hpp" -#include "csLogger.hpp" - -#include -#include -#include - -template class csMeanFreePath { - -private: - using levelSetsType = - lsSmartPointer>>>; - using cellSetType = lsSmartPointer>; - -public: - csMeanFreePath() : traceDevice(rtcNewDevice("hugepages=1")) { - static_assert(D == 2 && - "Mean free path calculation only implemented for 2D"); - } - - ~csMeanFreePath() { rtcReleaseDevice(traceDevice); } - - void setLevelSets(levelSetsType passedLevelSets) { - levelSets = passedLevelSets; - } - - void setCellSet(cellSetType passedCellSet) { - cellSet = passedCellSet; - } - - void setBulkLambda(const NumericType passedBulkLambda) { - bulkLambda = passedBulkLambda; - } - - template - void setMaterial(const Material passedMaterial) { - material = static_cast(passedMaterial); - } - - void setNumRaysPerCell(const int passedNumRaysPerCell) { - numRaysPerCell = passedNumRaysPerCell; - } - - NumericType getMaxLambda() const { return maxLambda; } - - void apply() { - csLogger::getInstance().addInfo("Calculating mean free path ...").print(); - cellSet->addScalarData("MeanFreePath", 0.); - runKernel(); - } - -private: - void runKernel() { -#ifdef ARCH_X86 - // for best performance set FTZ and DAZ flags in MXCSR control and status - // register - _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON); - _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); -#endif - unsigned numCells = cellSet->getElements().size(); - auto data = cellSet->getScalarData("MeanFreePath"); - auto materials = cellSet->getScalarData("Material"); - - auto traceGeometry = rayGeometry(); - { - auto mesh = lsSmartPointer>::New(); - lsToDiskMesh(levelSets->back(), mesh) - .apply(); - auto &points = mesh->getNodes(); - auto normals = mesh->getCellData().getVectorData("Normals"); - gridDelta = levelSets->back()->getGrid().getGridDelta(); - traceGeometry.initGeometry(traceDevice, points, *normals, - gridDelta * rayInternal::DiskFactor); - } - - auto rtcScene = rtcNewScene(traceDevice); - rtcSetSceneFlags(rtcScene, RTC_SCENE_FLAG_NONE); - rtcSetSceneBuildQuality(rtcScene, RTC_BUILD_QUALITY_HIGH); - auto rtcGeometry = traceGeometry.getRTCGeometry(); - auto geometryID = rtcAttachGeometry(rtcScene, rtcGeometry); - assert(rtcGetDeviceError(traceDevice) == RTC_ERROR_NONE && - "Embree device error"); - - maxLambda = 0.; - -#pragma omp parallel reduction(max : maxLambda) - { - rtcJoinCommitScene(rtcScene); - - alignas(128) auto rayHit = - RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; - -#if VIENNARAY_EMBREE_VERSION < 4 - auto rtcContext = RTCIntersectContext{}; - rtcInitIntersectContext(&rtcContext); -#endif - -#pragma omp for - for (int idx = 0; idx < numCells; ++idx) { - if (static_cast(materials->at(idx)) != material) - continue; - - auto cellCenter = cellSet->getCellCenter(idx); - auto &ray = rayHit.ray; -#ifdef ARCH_X86 - reinterpret_cast<__m128 &>(ray) = - _mm_set_ps(1e-4f, (float)cellCenter[2], (float)cellCenter[1], - (float)cellCenter[0]); -#else - ray.org_x = (float)cellCenter[0]; - ray.org_y = (float)cellCenter[1]; - ray.org_z = (float)cellCenter[2]; - ray.tnear = 1e-4f; -#endif - -#ifdef VIENNARAY_USE_RAY_MASKING - ray.mask = -1; -#endif - for (unsigned cIdx = 0; cIdx < numRaysPerCell; ++cIdx) { - - auto direction = getDirection(cIdx); - -#ifdef ARCH_X86 - reinterpret_cast<__m128 &>(ray.dir_x) = - _mm_set_ps(0.0f, (float)direction[2], (float)direction[1], - (float)direction[0]); -#else - ray.dir_x = (float)direction[0]; - ray.dir_y = (float)direction[1]; - ray.dir_z = (float)direction[2]; - ray.time = 0.0f; -#endif - - ray.tfar = std::numeric_limits::max(); - rayHit.hit.instID[0] = RTC_INVALID_GEOMETRY_ID; - rayHit.hit.geomID = RTC_INVALID_GEOMETRY_ID; - -#if VIENNARAY_EMBREE_VERSION < 4 - // Run the intersection - rtcIntersect1(rtcScene, &rtcContext, &rayHit); -#else - rtcIntersect1(rtcScene, &rayHit); -#endif - - /* -------- No hit -------- */ - if (rayHit.hit.geomID == RTC_INVALID_GEOMETRY_ID) { - data->at(idx) += bulkLambda; - continue; - } - - /* -------- Geometry hit -------- */ - const auto rayDir = - rayTriple{ray.dir_x, ray.dir_y, ray.dir_z}; - auto geomNormal = traceGeometry.getPrimNormal(rayHit.hit.primID); - - /* -------- Hit at disk backside -------- */ - if (rayInternal::DotProduct(rayDir, geomNormal) > 0) { - continue; - } - - data->at(idx) += ray.tfar; - } - - data->at(idx) /= numRaysPerCell; - maxLambda = std::max(maxLambda, data->at(idx)); - } - } // end of parallel section - - traceGeometry.releaseGeometry(); - rtcReleaseScene(rtcScene); - } - - std::array getDirection(const unsigned int idx) { - std::array direction; - NumericType theta = idx * 2. * M_PI / numRaysPerCell; - direction[0] = std::cos(theta); - direction[1] = std::sin(theta); - direction[2] = 0.; - return direction; - } - -private: - levelSetsType levelSets = nullptr; - cellSetType cellSet = nullptr; - RTCDevice traceDevice; - - NumericType gridDelta = 0; - NumericType bulkLambda = 0; - NumericType maxLambda = 0.; - long numRaysPerCell = 100; - int material = -1; -}; diff --git a/include/viennaps/cellSet/csPointSource.hpp b/include/viennaps/cellSet/csPointSource.hpp deleted file mode 100644 index b89bdbc..0000000 --- a/include/viennaps/cellSet/csPointSource.hpp +++ /dev/null @@ -1,41 +0,0 @@ -#pragma once - -#include "csUtil.hpp" - -#include - -template -class csPointSource : public raySource> { - const unsigned mNumPoints; - const csTriple origin; - const csTriple direction; - -public: - csPointSource(csTriple passedOrigin, - csTriple passedDirection, - std::array &pTraceSettings, const size_t pNumPoints) - : origin(passedOrigin), direction(passedDirection), - mNumPoints(pNumPoints) {} - - void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) const { -#ifdef ARCH_X86 - reinterpret_cast<__m128 &>(ray) = - _mm_set_ps(1e-4f, (float)origin[2], (float)origin[1], (float)origin[0]); - - reinterpret_cast<__m128 &>(ray.dir_x) = _mm_set_ps( - 0.0f, (float)direction[2], (float)direction[1], (float)direction[0]); -#else - ray.org_x = (float)origin[0]; - ray.org_y = (float)origin[1]; - ray.org_z = (float)origin[2]; - ray.tnear = 1e-4f; - - ray.dir_x = (float)direction[0]; - ray.dir_y = (float)direction[1]; - ray.dir_z = (float)direction[2]; - ray.time = 0.0f; -#endif - } - - size_t getNumPoints() const { return mNumPoints; } -}; diff --git a/include/viennaps/cellSet/csSegmentCells.hpp b/include/viennaps/cellSet/csSegmentCells.hpp deleted file mode 100644 index 3ac7ae7..0000000 --- a/include/viennaps/cellSet/csSegmentCells.hpp +++ /dev/null @@ -1,57 +0,0 @@ -#pragma once - -#include "csDenseCellSet.hpp" - -template class csSegmentCells { - lsSmartPointer> cellSet = nullptr; - std::string cellTypeString = "CellType"; - int bulkMaterial = -1; - -public: - csSegmentCells( - const lsSmartPointer> &passedCellSet) - : cellSet(passedCellSet) {} - - template - csSegmentCells( - const lsSmartPointer> &passedCellSet, - const std::string cellTypeString, const Material passedBulkMaterial) - : cellSet(passedCellSet), cellTypeString(cellTypeString), - bulkMaterial(static_cast(passedBulkMaterial)) {} - - void setCellSet( - const lsSmartPointer> &passedCellSet) { - cellSet = passedCellSet; - } - - void setCellTypeString(const std::string passedCellTypeString) { - cellTypeString = passedCellTypeString; - } - - template - void setBulkMaterial(const Material passedBulkMaterial) { - bulkMaterial = static_cast(passedBulkMaterial); - } - - void apply() { - auto cellType = cellSet->addScalarData(cellTypeString, -1.); - cellSet->buildNeighborhood(); - auto materials = cellSet->getScalarData("Material"); - -#pragma omp parallel for - for (int i = 0; i < materials->size(); ++i) { - if (static_cast(materials->at(i)) != bulkMaterial) { - auto neighbors = cellSet->getNeighbors(i); - for (auto n : neighbors) { - if (n >= 0 && - static_cast(materials->at(n)) == bulkMaterial) { - cellType->at(i) = 0.; - break; - } - } - } else { - cellType->at(i) = 1.; - } - } - } -}; diff --git a/include/viennaps/cellSet/csTracePath.hpp b/include/viennaps/cellSet/csTracePath.hpp deleted file mode 100644 index 52bb55d..0000000 --- a/include/viennaps/cellSet/csTracePath.hpp +++ /dev/null @@ -1,41 +0,0 @@ -#pragma once - -#include -#include - -template class csTracePath { -private: - std::unordered_map data; - std::vector gridData; - -public: - std::unordered_map &getData() { return data; } - - std::vector &getGridData() { return gridData; } - - T getGridValue(int idx) const { return gridData[idx]; } - - void addPoint(int idx, T value) { - auto search = data.find(idx); - if (search != data.end()) { - data[idx] += value; - } else { - data.insert({idx, value}); - } - } - - void addPath(csTracePath &path) { - for (const auto it : path.getData()) { - addPoint(it.first, it.second); - } - } - - void useGridData(size_t numCells) { gridData.resize(numCells, 0.); } - - void addGridData(int idx, T value) { gridData[idx] += value; } - - void clear() { - data.clear(); - gridData.clear(); - } -}; diff --git a/include/viennaps/cellSet/csTracing.hpp b/include/viennaps/cellSet/csTracing.hpp deleted file mode 100644 index ef30893..0000000 --- a/include/viennaps/cellSet/csTracing.hpp +++ /dev/null @@ -1,260 +0,0 @@ -#pragma once - -#include "csDenseCellSet.hpp" -#include "csPointSource.hpp" -#include "csTracingKernel.hpp" -#include "csTracingParticle.hpp" - -#include -#include -#include -#include -#include -#include - -template class csTracing { -private: - lsSmartPointer> cellSet = nullptr; - std::unique_ptr> mParticle = nullptr; - - RTCDevice mDevice; - rayGeometry mGeometry; - size_t mNumberOfRaysPerPoint = 0; - size_t mNumberOfRaysFixed = 1000; - T mGridDelta = 0; - rayBoundaryCondition mBoundaryConditions[D] = {}; - const rayTraceDirection mSourceDirection = - D == 2 ? rayTraceDirection::POS_Y : rayTraceDirection::POS_Z; - bool mUseRandomSeeds = true; - bool usePrimaryDirection = false; - rayTriple primaryDirection = {0.}; - size_t mRunNumber = 0; - int excludeMaterialId = -1; - bool usePointSource = false; - csTriple pointSourceOrigin = {0.}; - csTriple pointSourceDirection = {0.}; - -public: - csTracing() : mDevice(rtcNewDevice("hugepages=1")) { - // TODO: currently only periodic boundary conditions are implemented in - // csTracingKernel - for (int i = 0; i < D; i++) - mBoundaryConditions[i] = rayBoundaryCondition::PERIODIC; - } - - ~csTracing() { - mGeometry.releaseGeometry(); - rtcReleaseDevice(mDevice); - } - - void apply() { - createGeometry(); - initMemoryFlags(); - auto boundingBox = mGeometry.getBoundingBox(); - rayInternal::adjustBoundingBox( - boundingBox, mSourceDirection, mGridDelta * rayInternal::DiskFactor); - auto traceSettings = rayInternal::getTraceSettings(mSourceDirection); - - auto boundary = rayBoundary(mDevice, boundingBox, mBoundaryConditions, - traceSettings); - - std::array, 3> orthoBasis; - if (usePrimaryDirection) { - orthoBasis = rayInternal::getOrthonormalBasis(primaryDirection); - csLogger::getInstance() - .addInfo("Using primary direction: " + - std::to_string(primaryDirection[0]) + " " + - std::to_string(primaryDirection[1]) + " " + - std::to_string(primaryDirection[2])) - .print(); - } - - if (usePointSource) { - auto raySource = - csPointSource(pointSourceOrigin, pointSourceDirection, - traceSettings, mGeometry.getNumPoints()); - - csTracingKernel(mDevice, mGeometry, boundary, raySource, mParticle, - mNumberOfRaysPerPoint, mNumberOfRaysFixed, - mUseRandomSeeds, mRunNumber++, cellSet, - excludeMaterialId - 1) - .apply(); - } else { - auto raySource = raySourceRandom( - boundingBox, mParticle->getSourceDistributionPower(), traceSettings, - mGeometry.getNumPoints(), usePrimaryDirection, orthoBasis); - - csTracingKernel(mDevice, mGeometry, boundary, raySource, mParticle, - mNumberOfRaysPerPoint, mNumberOfRaysFixed, - mUseRandomSeeds, mRunNumber++, cellSet, - excludeMaterialId - 1) - .apply(); - } - - averageNeighborhood(); - boundary.releaseGeometry(); - } - - void setCellSet(lsSmartPointer> passedCellSet) { - cellSet = passedCellSet; - } - - void setPointSource(const csTriple &passedOrigin, - const csTriple &passedDirection) { - usePointSource = true; - pointSourceOrigin = passedOrigin; - pointSourceDirection = passedDirection; - } - - template - void setParticle(std::unique_ptr &p) { - static_assert(std::is_base_of, ParticleType>::value && - "Particle object does not interface correct class"); - mParticle = p->clone(); - } - - void setTotalNumberOfRays(const size_t passedNumber) { - mNumberOfRaysFixed = passedNumber; - mNumberOfRaysPerPoint = 0; - } - - void setNumberOfRaysPerPoint(const size_t passedNumber) { - mNumberOfRaysPerPoint = passedNumber; - mNumberOfRaysFixed = 0; - } - - void setPrimaryDirection(const rayTriple pPrimaryDirection) { - primaryDirection = pPrimaryDirection; - usePrimaryDirection = true; - } - - void setExcludeMaterialId(int passedId) { excludeMaterialId = passedId; } - - lsSmartPointer> getCellSet() const { return cellSet; } - - void averageNeighborhood() { - auto data = cellSet->getFillingFractions(); - auto materialIds = cellSet->getScalarData("Material"); - const auto &elems = cellSet->getElements(); - const auto &nodes = cellSet->getNodes(); - std::vector average(data->size(), 0.); - -#pragma omp parallel for - for (int i = 0; i < data->size(); i++) { - if (data->at(i) < 0) { - average[i] = -1.; - continue; - } - if (materialIds->at(i) == excludeMaterialId) - continue; - - int numNeighbors = 1; - average[i] += data->at(i); - - for (int d = 0; d < D; d++) { - auto mid = calcMidPoint(nodes[elems[i][0]]); - mid[d] -= mGridDelta; - auto elemId = cellSet->getIndex(mid); - if (elemId >= 0) { - if (data->at(elemId) >= 0 && - materialIds->at(elemId) != excludeMaterialId) { - average[i] += data->at(elemId); - numNeighbors++; - } - } - - mid[d] += 2 * mGridDelta; - elemId = cellSet->getIndex(mid); - if (elemId >= 0) { - if (data->at(elemId) >= 0 && - materialIds->at(elemId) != excludeMaterialId) { - average[i] += data->at(elemId); - numNeighbors++; - } - } - } - average[i] /= static_cast(numNeighbors); - } - -#pragma omp parallel for - for (int i = 0; i < data->size(); i++) { - data->at(i) = average[i]; - } - } - - void averageNeighborhoodSingleMaterial(int materialId) { - auto data = cellSet->getFillingFractions(); - auto materialIds = cellSet->getScalarData("Material"); - const auto &elems = cellSet->getElements(); - const auto &nodes = cellSet->getNodes(); - std::vector average(data->size(), 0.); - -#pragma omp parallel for - for (int i = 0; i < data->size(); i++) { - if (materialIds->at(i) != materialId) - continue; - - int numNeighbors = 1; - average[i] += data->at(i); - for (int d = 0; d < D; d++) { - auto mid = calcMidPoint(nodes[elems[i][0]]); - mid[d] -= mGridDelta; - auto elemId = cellSet->getIndex(mid); - if (elemId >= 0) { - if (data->at(elemId) >= 0 && materialIds->at(elemId) == materialId) { - average[i] += data->at(elemId); - numNeighbors++; - } - } - - mid[d] += 2 * mGridDelta; - elemId = cellSet->getIndex(mid); - if (elemId >= 0) { - if (data->at(elemId) >= 0 && materialIds->at(elemId) == materialId) { - average[i] += data->at(elemId); - numNeighbors++; - } - } - } - average[i] /= static_cast(numNeighbors); - } - -#pragma omp parallel for - for (int i = 0; i < data->size(); i++) { - data->at(i) = average[i]; - } - } - -private: - void createGeometry() { - auto levelSets = cellSet->getLevelSets(); - auto diskMesh = lsSmartPointer>::New(); - lsToDiskMesh converter(diskMesh); - for (auto ls : *levelSets) { - converter.insertNextLevelSet(ls); - } - converter.apply(); - auto points = diskMesh->getNodes(); - auto normals = *diskMesh->getCellData().getVectorData("Normals"); - auto materialIds = *diskMesh->getCellData().getScalarData("MaterialIds"); - mGridDelta = levelSets->back()->getGrid().getGridDelta(); - mGeometry.initGeometry(mDevice, points, normals, - mGridDelta * rayInternal::DiskFactor); - mGeometry.setMaterialIds(materialIds); - } - - inline csTriple calcMidPoint(const csTriple &minNode) { - return csTriple{minNode[0] + mGridDelta / T(2), - minNode[1] + mGridDelta / T(2), - minNode[2] + mGridDelta / T(2)}; - } - - void initMemoryFlags() { -#ifdef ARCH_X86 - // for best performance set FTZ and DAZ flags in MXCSR control and status - // register - _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON); - _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); -#endif - } -}; diff --git a/include/viennaps/cellSet/csTracingKernel.hpp b/include/viennaps/cellSet/csTracingKernel.hpp deleted file mode 100644 index f80481d..0000000 --- a/include/viennaps/cellSet/csTracingKernel.hpp +++ /dev/null @@ -1,297 +0,0 @@ -#pragma once - -#include "csDenseCellSet.hpp" -#include "csTracePath.hpp" -#include "csTracingParticle.hpp" - -#include -#include -#include -#include -#include -#include - -template class csTracingKernel { -public: - csTracingKernel(RTCDevice &pDevice, rayGeometry &pRTCGeometry, - rayBoundary &pRTCBoundary, - raySource &pSource, - std::unique_ptr> &pParticle, - const size_t pNumOfRayPerPoint, const size_t pNumOfRayFixed, - const bool pUseRandomSeed, const size_t pRunNumber, - lsSmartPointer> passedCellSet, - int passedExclude) - : mDevice(pDevice), mGeometry(pRTCGeometry), mBoundary(pRTCBoundary), - mSource(pSource), mParticle(pParticle->clone()), - mNumRays(pNumOfRayFixed == 0 - ? pSource.getNumPoints() * pNumOfRayPerPoint - : pNumOfRayFixed), - mUseRandomSeeds(pUseRandomSeed), mRunNumber(pRunNumber), - cellSet(passedCellSet), excludeMaterial(passedExclude), - mGridDelta(cellSet->getGridDelta()) { - assert(rtcGetDeviceProperty(mDevice, RTC_DEVICE_PROPERTY_VERSION) >= - 30601 && - "Error: The minimum version of Embree is 3.6.1"); - } - - void apply() { - auto rtcScene = rtcNewScene(mDevice); - rtcSetSceneFlags(rtcScene, RTC_SCENE_FLAG_NONE); - auto bbquality = RTC_BUILD_QUALITY_HIGH; - rtcSetSceneBuildQuality(rtcScene, bbquality); - auto rtcGeometry = mGeometry.getRTCGeometry(); - auto rtcBoundary = mBoundary.getRTCGeometry(); - - auto boundaryID = rtcAttachGeometry(rtcScene, rtcBoundary); - auto geometryID = rtcAttachGeometry(rtcScene, rtcGeometry); - assert(rtcGetDeviceError(mDevice) == RTC_ERROR_NONE && - "Embree device error"); - - const csPair meanFreePath = mParticle->getMeanFreePath(); - - auto myCellSet = cellSet; - -#pragma omp parallel shared(myCellSet) - { - rtcJoinCommitScene(rtcScene); - - alignas(128) auto rayHit = - RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; - - const int threadID = omp_get_thread_num(); - unsigned int seed = mRunNumber; - if (mUseRandomSeeds) { - std::random_device rd; - seed = static_cast(rd()); - } - - // thread-local particle object - auto particle = mParticle->clone(); - - // thread local path - csTracePath path; - // if (!traceOnPath) - path.useGridData(myCellSet->getNumberOfCells()); - -#if VIENNARAY_EMBREE_VERSION < 4 - auto rtcContext = RTCIntersectContext{}; - rtcInitIntersectContext(&rtcContext); -#endif - -#pragma omp for schedule(dynamic) - for (long long idx = 0; idx < mNumRays; ++idx) { - // particle specific RNG seed - auto particleSeed = rayInternal::tea<3>(idx, seed); - rayRNG RngState(particleSeed); - - particle->initNew(RngState); - - mSource.fillRay(rayHit.ray, idx, RngState); // fills also tnear - -#ifdef VIENNARAY_USE_RAY_MASKING - rayHit.ray.mask = -1; -#endif - - bool reflect = false; - bool hitFromBack = false; - do { - rayHit.ray.tfar = - std::numeric_limits::max(); - rayHit.hit.instID[0] = RTC_INVALID_GEOMETRY_ID; - rayHit.hit.geomID = RTC_INVALID_GEOMETRY_ID; - -#if VIENNARAY_EMBREE_VERSION < 4 - // Run the intersection - rtcIntersect1(rtcScene, &rtcContext, &rayHit); -#else - rtcIntersect1(rtcScene, &rayHit); -#endif - - /* -------- No hit -------- */ - if (rayHit.hit.geomID == RTC_INVALID_GEOMETRY_ID) { - reflect = false; - break; - } - - /* -------- Boundary hit -------- */ - if (rayHit.hit.geomID == boundaryID) { - mBoundary.processHit(rayHit, reflect); - continue; - } - - // Calculate point of impact - const auto &ray = rayHit.ray; - const rayInternal::rtcNumericType xx = - ray.org_x + ray.dir_x * ray.tfar; - const rayInternal::rtcNumericType yy = - ray.org_y + ray.dir_y * ray.tfar; - const rayInternal::rtcNumericType zz = - ray.org_z + ray.dir_z * ray.tfar; - - /* -------- Hit from back -------- */ - const auto rayDir = rayTriple{ray.dir_x, ray.dir_y, ray.dir_z}; - const auto geomNormal = mGeometry.getPrimNormal(rayHit.hit.primID); - if (rayInternal::DotProduct(rayDir, geomNormal) > 0) { - // If the dot product of the ray direction and the surface normal is - // greater than zero, then we hit the back face of the disk. - if (hitFromBack) { - // if hitFromback == true, then the ray hits the back of a disk - // the second time. In this case we ignore the ray. - break; - } - hitFromBack = true; - // Let ray through, i.e., continue. - reflect = true; -#ifdef ARCH_X86 - reinterpret_cast<__m128 &>(rayHit.ray) = - _mm_set_ps(1e-4f, zz, yy, xx); -#else - rayHit.ray.org_x = xx; - rayHit.ray.org_y = yy; - rayHit.ray.org_z = zz; - rayHit.ray.tnear = 1e-4f; -#endif - // keep ray direction as it is - continue; - } - - /* -------- Surface hit -------- */ - assert(rayHit.hit.geomID == geometryID && "Geometry hit ID invalid"); - - // get fill and reflection - const auto fillnDirection = - particle->surfaceHit(rayDir, geomNormal, reflect, RngState); - - if (mGeometry.getMaterialId(rayHit.hit.primID) != excludeMaterial) { - // trace in cell set - auto hitPoint = std::array{xx, yy, zz}; - std::vector> particleStack; - std::normal_distribution normalDist{meanFreePath[0], - meanFreePath[1]}; - - particleStack.emplace_back(csVolumeParticle{ - hitPoint, rayDir, fillnDirection.first, 0., -1, 0}); - - while (!particleStack.empty()) { - auto volumeParticle = std::move(particleStack.back()); - particleStack.pop_back(); - - // trace particle - while (volumeParticle.energy >= 0) { - volumeParticle.distance = -1; - while (volumeParticle.distance < 0) - volumeParticle.distance = normalDist(RngState); - auto travelDist = csUtil::multNew(volumeParticle.direction, - volumeParticle.distance); - csUtil::add(volumeParticle.position, travelDist); - - if (!checkBoundsPeriodic(volumeParticle.position)) - break; - - auto newIdx = myCellSet->getIndex(volumeParticle.position); - if (newIdx < 0) - break; - - if (newIdx != volumeParticle.cellId) { - volumeParticle.cellId = newIdx; - auto fill = particle->collision(volumeParticle, RngState, - particleStack); - path.addGridData(newIdx, fill); - } - } - } - } - - if (!reflect) { - break; - } - - // Update ray direction and origin -#ifdef ARCH_X86 - reinterpret_cast<__m128 &>(rayHit.ray) = - _mm_set_ps(1e-4f, zz, yy, xx); - reinterpret_cast<__m128 &>(rayHit.ray.dir_x) = _mm_set_ps( - 0.0f, (rayInternal::rtcNumericType)fillnDirection.second[2], - (rayInternal::rtcNumericType)fillnDirection.second[1], - (rayInternal::rtcNumericType)fillnDirection.second[0]); -#else - rayHit.ray.org_x = xx; - rayHit.ray.org_y = yy; - rayHit.ray.org_z = zz; - rayHit.ray.tnear = 1e-4f; - - rayHit.ray.dir_x = - (rayInternal::rtcNumericType)fillnDirection.second[0]; - rayHit.ray.dir_y = - (rayInternal::rtcNumericType)fillnDirection.second[1]; - rayHit.ray.dir_z = - (rayInternal::rtcNumericType)fillnDirection.second[2]; - rayHit.ray.time = 0.0f; -#endif - } while (reflect); - } // end ray tracing for loop - -#pragma omp critical - myCellSet->mergePath(path, mNumRays); - } // end parallel section - - rtcReleaseGeometry(rtcGeometry); - rtcReleaseGeometry(rtcBoundary); - } - -private: - bool checkBounds(const csTriple &hitPoint) const { - const auto &min = cellSet->getCellGrid()->minimumExtent; - const auto &max = cellSet->getCellGrid()->maximumExtent; - - return hitPoint[0] >= min[0] && hitPoint[0] <= max[0] && - hitPoint[1] >= min[1] && hitPoint[1] <= max[1] && - hitPoint[2] >= min[2] && hitPoint[2] <= max[2]; - } - - bool checkBoundsPeriodic(csTriple &hitPoint) const { - const auto &min = cellSet->getCellGrid()->minimumExtent; - const auto &max = cellSet->getCellGrid()->maximumExtent; - - if constexpr (D == 3) { - if (hitPoint[2] < min[2] || hitPoint[2] > max[2]) - return false; - - if (hitPoint[0] < min[0]) { - hitPoint[0] = max[0] - mGridDelta / 2.; - } else if (hitPoint[0] > max[0]) { - hitPoint[0] = min[0] + mGridDelta / 2.; - } - - if (hitPoint[1] < min[1]) { - hitPoint[1] = max[1] - mGridDelta / 2.; - } else if (hitPoint[1] > max[1]) { - hitPoint[1] = min[1] + mGridDelta / 2.; - } - } else { - if (hitPoint[1] < min[1] || hitPoint[1] > max[1]) - return false; - - if (hitPoint[0] < min[0]) { - hitPoint[0] = max[0] - mGridDelta / 2.; - } else if (hitPoint[0] > max[0]) { - hitPoint[0] = min[0] + mGridDelta / 2.; - } - } - - return true; - } - -private: - RTCDevice &mDevice; - rayGeometry &mGeometry; - rayBoundary &mBoundary; - raySource &mSource; - std::unique_ptr> const mParticle = nullptr; - const long long mNumRays; - const bool mUseRandomSeeds; - const size_t mRunNumber; - lsSmartPointer> cellSet = nullptr; - const T mGridDelta = 0.; - const int excludeMaterial = -1; -}; diff --git a/include/viennaps/cellSet/csTracingParticle.hpp b/include/viennaps/cellSet/csTracingParticle.hpp deleted file mode 100644 index bc391ce..0000000 --- a/include/viennaps/cellSet/csTracingParticle.hpp +++ /dev/null @@ -1,52 +0,0 @@ -#pragma once - -#include "csUtil.hpp" - -#include -#include -#include - -template class csAbstractParticle { -public: - virtual ~csAbstractParticle() = default; - virtual std::unique_ptr clone() const = 0; - - virtual void initNew(rayRNG &Rng) = 0; - - virtual std::pair> surfaceHit(const rayTriple &rayDir, - const rayTriple &geomNormal, - bool &reflect, rayRNG &Rng) = 0; - virtual T getSourceDistributionPower() const = 0; - virtual csPair getMeanFreePath() const = 0; - virtual T collision(csVolumeParticle &particle, rayRNG &RNG, - std::vector> &particleStack) = 0; -}; - -template -class csParticle : public csAbstractParticle { -public: - std::unique_ptr> clone() const override final { - return std::make_unique(static_cast(*this)); - } - virtual void initNew(rayRNG &Rng) override {} - virtual std::pair> surfaceHit(const rayTriple &rayDir, - const rayTriple &geomNormal, - bool &reflect, - rayRNG &Rng) override { - reflect = false; - return std::pair>{1., rayTriple{0., 0., 0.}}; - } - virtual T getSourceDistributionPower() const override { return 1.; } - virtual csPair getMeanFreePath() const override { return {1., 1.}; } - virtual T - collision(csVolumeParticle &particle, rayRNG &RNG, - std::vector> &particleStack) override { - return 0.; - } - -protected: - // We make clear csParticle class needs to be inherited - csParticle() = default; - csParticle(const csParticle &) = default; - csParticle(csParticle &&) = default; -}; diff --git a/include/viennaps/cellSet/csUtil.hpp b/include/viennaps/cellSet/csUtil.hpp deleted file mode 100644 index fa3899c..0000000 --- a/include/viennaps/cellSet/csUtil.hpp +++ /dev/null @@ -1,183 +0,0 @@ -#pragma once - -#include -#include -#include -#include -#include -#include -#include - -template using csPair = std::array; - -template using csTriple = std::array; - -template struct csVolumeParticle { - csTriple position; - csTriple direction; - T energy; - T distance; - int cellId; - int scattered; -}; - -namespace csUtil { - -template struct Timer { - using TimePoint = typename Clock::time_point; - - TimePoint mStart; - typename Clock::duration::rep totalDuration = 0.; // in ns - typename Clock::duration::rep currentDuration; // in ns - - void start() { mStart = Clock::now(); } - void finish() { - TimePoint end = Clock::now(); - typename Clock::duration dur(end - mStart); - currentDuration = dur.count(); - totalDuration += currentDuration; - } - void reset() { - currentDuration = 0.; - totalDuration = 0.; - } -}; - -template void printTriple(const csTriple &p) { - std::cout << "[" << p[0] << ", " << p[1] << ", " << p[2] << "]\n"; -} - -template -inline T dot(const std::array &v1, const std::array &v2) { - T res = 0; - for (int i = 0; i < D; i++) - res += v1[i] * v2[i]; - return res; -} - -template -inline void mult(std::array &v1, const T fac) { - for (int i = 0; i < D; i++) - v1[i] *= fac; -} - -template -inline std::array multNew(const std::array &v1, const T fac) { - std::array res; - for (int i = 0; i < D; i++) - res[i] = v1[i] * fac; - return res; -} - -template -inline void add(std::array &vec, const std::array &vec2) { - for (size_t i = 0; i < D; i++) - vec[i] += vec2[i]; -} - -template -inline void sub(std::array &v1, const std::array &v2) { - for (int i = 0; i < D; i++) - v1[i] -= v2[i]; -} - -template -inline void multAdd(std::array &result, const std::array &mult, - const std::array &add, const T fac) { - for (int i = 0; i < D; i++) - result[i] = add[i] + mult[i] * fac; -} - -template -inline T distance(const std::array &p1, const std::array &p2) { - T res = 0; - for (int i = 0; i < D; i++) - res += (p1[i] - p2[i]) * (p1[i] - p2[i]); - return std::sqrt(res); -} - -template inline T norm(const std::array &vec) { - T res = 0; - for (int i = 0; i < D; i++) - res += vec[i] * vec[i]; - return std::sqrt(res); -} - -template void normalize(std::array &vec) { - T vecNorm = 1. / norm(vec); - if (vecNorm == 1.) - return; - std::for_each(vec.begin(), vec.end(), - [&vecNorm](T &entry) { entry *= vecNorm; }); -} - -template -void scaleToLength(std::array &vec, const T length) { - const auto vecLength = norm(vec); - for (size_t i = 0; i < D; i++) - vec[i] *= length / vecLength; -} - -template -csTriple crossProd(const csTriple &pVecA, const csTriple &pVecB) { - csTriple rr; - rr[0] = pVecA[1] * pVecB[2] - pVecA[2] * pVecB[1]; - rr[1] = pVecA[2] * pVecB[0] - pVecA[0] * pVecB[2]; - rr[2] = pVecA[0] * pVecB[1] - pVecA[1] * pVecB[0]; - return rr; -} - -#ifdef ARCH_X86 -[[nodiscard]] static inline float DotProductSse(__m128 const &x, - __m128 const &y) { - // __m128 mulRes; - // mulRes = _mm_mul_ps(x, y); - // return SumSse(mulRes); - return _mm_cvtss_f32(_mm_dp_ps(x, y, 0x77)); -} - -[[nodiscard]] static inline float SumSse(__m128 const &v) { - __m128 shufReg, sumsReg; - // Calculates the sum of SSE Register - - // https://stackoverflow.com/a/35270026/195787 - shufReg = _mm_movehdup_ps(v); // Broadcast elements 3,1 to 2,0 - sumsReg = _mm_add_ps(v, shufReg); - shufReg = _mm_movehl_ps(shufReg, sumsReg); // High Half -> Low Half - sumsReg = _mm_add_ss(sumsReg, shufReg); - return _mm_cvtss_f32(sumsReg); // Result in the lower part of the SSE Register -} - -[[nodiscard]] inline static __m128 CrossProductSse(__m128 const &vec0, - __m128 const &vec1) { - // from geometrian.com/programming/tutorials/cross-product/index.php - __m128 tmp0 = _mm_shuffle_ps(vec0, vec0, _MM_SHUFFLE(3, 0, 2, 1)); - __m128 tmp1 = _mm_shuffle_ps(vec1, vec1, _MM_SHUFFLE(3, 1, 0, 2)); - __m128 tmp2 = _mm_mul_ps(tmp0, vec1); - __m128 tmp3 = _mm_mul_ps(tmp0, tmp1); - __m128 tmp4 = _mm_shuffle_ps(tmp2, tmp2, _MM_SHUFFLE(3, 0, 2, 1)); - return _mm_sub_ps(tmp3, tmp4); -} - -// Norm of 3D Vector using SSE -// fastcpp.blogspot.com/2012/02/calculating-length-of-3d-vector-using.html -[[nodiscard]] static inline float NormSse(__m128 const &v) { - return _mm_cvtss_f32(_mm_sqrt_ss(_mm_dp_ps(v, v, 0x71))); -} - -[[nodiscard]] static inline __m128 NormalizeSse(__m128 const &v) { - __m128 inverse_norm = _mm_rsqrt_ps(_mm_dp_ps(v, v, 0x77)); - return _mm_mul_ps(v, inverse_norm); -} - -[[nodiscard]] static inline __m128 NormalizeAccurateSse(__m128 const &v) { - __m128 norm = _mm_sqrt_ps(_mm_dp_ps(v, v, 0x7F)); - return _mm_div_ps(v, norm); -} - -template [[nodiscard]] rayTriple ConvertSse(__m128 const &vec) { - alignas(16) float result[4]; - _mm_store_ps(&result[0], vec); - return csTriple{result[0], result[1], result[2]}; -} -#endif -} // namespace csUtil diff --git a/include/viennaps/models/psOxideRegrowth.hpp b/include/viennaps/models/psOxideRegrowth.hpp index 765aa26..2cbef6d 100644 --- a/include/viennaps/models/psOxideRegrowth.hpp +++ b/include/viennaps/models/psOxideRegrowth.hpp @@ -1,11 +1,11 @@ #pragma once -#include "../cellSet/csDenseCellSet.hpp" - #include "../psAdvectionCallback.hpp" #include "../psProcessModel.hpp" #include "../psToDiskMesh.hpp" +#include + #include namespace OxideRegrowthImplementation { diff --git a/include/viennaps/models/psPlasmaDamage.hpp b/include/viennaps/models/psPlasmaDamage.hpp index dc26ef7..d0b8767 100644 --- a/include/viennaps/models/psPlasmaDamage.hpp +++ b/include/viennaps/models/psPlasmaDamage.hpp @@ -1,9 +1,10 @@ #pragma once -#include "../cellSet/csTracing.hpp" -#include "../cellSet/csTracingParticle.hpp" #include "../psProcessModel.hpp" +#include +#include + #include namespace PlasmaDamageImplementation { diff --git a/include/viennaps/psDomain.hpp b/include/viennaps/psDomain.hpp index afb01b4..f24d459 100644 --- a/include/viennaps/psDomain.hpp +++ b/include/viennaps/psDomain.hpp @@ -1,12 +1,12 @@ #pragma once -#include "cellSet/csDenseCellSet.hpp" - #include "psMaterials.hpp" #include "psSmartPointer.hpp" #include "psSurfacePointValuesToLevelSet.hpp" #include "psVTKWriter.hpp" +#include + #include #include #include From 290952c1ace2947579ea990eb6a1559cb56bb709 Mon Sep 17 00:00:00 2001 From: Felix S Date: Mon, 27 May 2024 17:28:36 +0200 Subject: [PATCH 6/6] Move cell set to ViennaCS library - CMake changes --- CMakeLists.txt | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6f31e8a..674dd4b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -80,6 +80,7 @@ enable_language(C) find_package(ViennaRay QUIET) find_package(ViennaLS QUIET) +find_package(ViennaCS QUIET) set(VIENNAPS_SYSTEM_VIENNARAY ${ViennaRay_FOUND} @@ -89,6 +90,10 @@ set(VIENNAPS_SYSTEM_VIENNALS ${ViennaLS_FOUND} CACHE INTERNAL "") +set(VIENNAPS_SYSTEM_VIENNACS + ${ViennaCS_FOUND} + CACHE INTERNAL "") + # -------------------------------------------------------------------------------------------------------- # Setup Dependencies # -------------------------------------------------------------------------------------------------------- @@ -112,7 +117,13 @@ CPMFindPackage( GIT_REPOSITORY "https://github.com/ViennaTools/ViennaLS" EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON}) -target_link_libraries(${PROJECT_NAME} INTERFACE ViennaTools::ViennaRay ViennaTools::ViennaLS) +CPMFindPackage( + NAME ViennaCS + VERSION 1.0.0 + GIT_REPOSITORY "https://github.com/ViennaTools/ViennaCS" + EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON}) + +target_link_libraries(${PROJECT_NAME} INTERFACE ViennaTools::ViennaCS) # -------------------------------------------------------------------------------------------------------- # Setup Application @@ -164,4 +175,4 @@ packageProject( INCLUDE_DIR ${PROJECT_SOURCE_DIR}/include/viennaps INCLUDE_DESTINATION include/viennaps-${PROJECT_VERSION} COMPATIBILITY SameMajorVersion - DEPENDENCIES "ViennaLS;ViennaRay") + DEPENDENCIES "ViennaCS;ViennaLS;ViennaRay")