ViennaLS
Loading...
Searching...
No Matches
lsToVoxelMesh.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <unordered_map>
4
6
7#include <hrleDenseCellIterator.hpp>
8#include <lsDomain.hpp>
9#include <lsMaterialMap.hpp>
10#include <lsMesh.hpp>
11
12namespace viennals {
13
14using namespace viennacore;
15
20template <class T, int D> class ToVoxelMesh {
21 typedef typename Domain<T, D>::DomainType hrleDomainType;
22
23 std::vector<SmartPointer<Domain<T, D>>> levelSets;
24 SmartPointer<Mesh<T>> mesh = nullptr;
25 SmartPointer<MaterialMap> materialMap = nullptr;
26 viennahrle::Index<D> minIndex, maxIndex;
27
28 void calculateBounds() {
29 // set to zero
30 for (unsigned i = 0; i < D; ++i) {
31 minIndex[i] = std::numeric_limits<viennahrle::IndexType>::max();
32 maxIndex[i] = std::numeric_limits<viennahrle::IndexType>::lowest();
33 }
34 for (unsigned l = 0; l < levelSets.size(); ++l) {
35 auto &grid = levelSets[l]->getGrid();
36 auto &domain = levelSets[l]->getDomain();
37 for (unsigned i = 0; i < D; ++i) {
38 minIndex[i] = std::min(minIndex[i], (grid.isNegBoundaryInfinite(i))
39 ? domain.getMinRunBreak(i)
40 : grid.getMinBounds(i));
41
42 maxIndex[i] = std::max(maxIndex[i], (grid.isPosBoundaryInfinite(i))
43 ? domain.getMaxRunBreak(i)
44 : grid.getMaxBounds(i));
45 }
46 }
47 }
48
49public:
50 ToVoxelMesh() = default;
51
52 ToVoxelMesh(SmartPointer<Mesh<T>> passedMesh) : mesh(passedMesh) {}
53
54 ToVoxelMesh(SmartPointer<Domain<T, D>> passedLevelSet,
55 SmartPointer<Mesh<T>> passedMesh)
56 : mesh(passedMesh) {
57 levelSets.push_back(passedLevelSet);
58 }
59
60 ToVoxelMesh(const std::vector<SmartPointer<Domain<T, D>>> passedLevelSets,
61 SmartPointer<Mesh<T>> passedMesh)
62 : mesh(passedMesh) {
63 levelSets = passedLevelSets;
64 }
65
70 void insertNextLevelSet(SmartPointer<Domain<T, D>> passedLevelSet) {
71 levelSets.push_back(passedLevelSet);
72 }
73
74 void setMesh(SmartPointer<Mesh<T>> passedMesh) { mesh = passedMesh; }
75
76 void setMaterialMap(SmartPointer<MaterialMap> passedMaterialMap) {
77 materialMap = passedMaterialMap;
78 }
79
80 void apply() {
81 if (levelSets.size() < 1) {
82 Logger::getInstance()
83 .addWarning("No Level Sets supplied to ToVoxelMesh! Not converting.")
84 .print();
85 }
86 if (mesh == nullptr) {
87 Logger::getInstance()
88 .addWarning("No mesh was passed to ToVoxelMesh.")
89 .print();
90 return;
91 }
92
93 mesh->clear();
94 auto &levelSet = *(levelSets.back());
95 auto &grid = levelSet.getGrid();
96
97 calculateBounds();
98
99 // save the extent of the resulting mesh
100 for (unsigned i = 0; i < D; ++i) {
101 mesh->minimumExtent[i] = std::numeric_limits<T>::max();
102 mesh->maximumExtent[i] = std::numeric_limits<T>::lowest();
103 }
104
105 std::unordered_map<viennahrle::Index<D>, size_t,
106 typename viennahrle::Index<D>::hash>
107 pointIdMapping;
108 size_t currentPointId = 0;
109
110 // prepare mesh for material ids
111 mesh->cellData.insertNextScalarData(typename PointData<T>::ScalarDataType(),
112 "Material");
113 auto &materialIds = *(mesh->cellData.getScalarData(0));
114 const bool useMaterialMap = materialMap != nullptr;
115
116 // set up iterators for all materials
117 std::vector<
118 viennahrle::ConstDenseCellIterator<typename Domain<T, D>::DomainType>>
119 iterators;
120 for (auto it = levelSets.begin(); it != levelSets.end(); ++it) {
121 iterators.emplace_back((*it)->getDomain(), minIndex);
122 }
123
124 // move iterator for lowest material id and then adjust others if they are
125 // needed
126 for (; iterators.front().getIndices() < maxIndex;
127 iterators.front().next()) {
128 // go over all materials
129 for (unsigned materialId = 0; materialId < levelSets.size();
130 ++materialId) {
131
132 auto &cellIt = iterators[materialId];
133
134 cellIt.goToIndicesSequential(iterators.front().getIndices());
135
136 // find out whether the centre of the box is inside
137 T centerValue = 0.;
138 for (int i = 0; i < (1 << D); ++i) {
139 centerValue += cellIt.getCorner(i).getValue();
140 }
141
142 if (centerValue <= 0.) {
143 std::array<unsigned, 1 << D> voxel;
144 bool addVoxel = false;
145 // now insert all points of voxel into pointList
146 for (unsigned i = 0; i < (1 << D); ++i) {
147 viennahrle::Index<D> index;
148 addVoxel = true;
149 for (unsigned j = 0; j < D; ++j) {
150 index[j] =
151 cellIt.getIndices(j) + cellIt.getCorner(i).getOffset()[j];
152 if (index[j] > maxIndex[j]) {
153 addVoxel = false;
154 break;
155 }
156 }
157 if (addVoxel) {
158 auto pointIdValue = std::make_pair(index, currentPointId);
159 auto pointIdPair = pointIdMapping.insert(pointIdValue);
160 voxel[i] = pointIdPair.first->second;
161 if (pointIdPair.second) {
162 ++currentPointId;
163 }
164 } else {
165 break;
166 }
167 }
168
169 // create element if inside domain bounds
170 if (addVoxel) {
171 int material = materialId;
172 if (useMaterialMap)
173 material = materialMap->getMaterialId(materialId);
174
175 if constexpr (D == 3) {
176 // reorder elements for hexas to be ordered correctly
177 std::array<unsigned, 8> hexa{voxel[0], voxel[1], voxel[3],
178 voxel[2], voxel[4], voxel[5],
179 voxel[7], voxel[6]};
180 mesh->hexas.push_back(hexa);
181 materialIds.push_back(material);
182 } else {
183 std::array<unsigned, 4> tetra{voxel[0], voxel[2], voxel[3],
184 voxel[1]};
185 mesh->tetras.push_back(tetra);
186 materialIds.push_back(material);
187 }
188 }
189 // jump out of material for loop
190 break;
191 }
192 }
193 }
194
195 // now insert points
196 double gridDelta = grid.getGridDelta();
197 mesh->nodes.resize(pointIdMapping.size());
198 for (auto it = pointIdMapping.begin(); it != pointIdMapping.end(); ++it) {
199 Vec3D<T> coords;
200 for (unsigned i = 0; i < D; ++i) {
201 coords[i] = gridDelta * it->first[i];
202
203 // save extent
204 if (coords[i] < mesh->minimumExtent[i]) {
205 mesh->minimumExtent[i] = coords[i];
206 } else if (coords[i] > mesh->maximumExtent[i]) {
207 mesh->maximumExtent[i] = coords[i];
208 }
209 }
210 mesh->nodes[it->second] = coords;
211 }
212 }
213};
214
215// add all template specializations for this class
217
218} // namespace viennals
Class containing all information about the level set, including the dimensions of the domain,...
Definition lsDomain.hpp:27
viennahrle::Domain< T, D > DomainType
Definition lsDomain.hpp:32
This class holds an explicit mesh, which is always given in 3 dimensions. If it describes a 2D mesh,...
Definition lsMesh.hpp:21
std::vector< T > ScalarDataType
Definition lsPointData.hpp:22
void setMesh(SmartPointer< Mesh< T > > passedMesh)
Definition lsToVoxelMesh.hpp:74
ToVoxelMesh(SmartPointer< Mesh< T > > passedMesh)
Definition lsToVoxelMesh.hpp:52
void setMaterialMap(SmartPointer< MaterialMap > passedMaterialMap)
Definition lsToVoxelMesh.hpp:76
ToVoxelMesh(SmartPointer< Domain< T, D > > passedLevelSet, SmartPointer< Mesh< T > > passedMesh)
Definition lsToVoxelMesh.hpp:54
void insertNextLevelSet(SmartPointer< Domain< T, D > > passedLevelSet)
Push level set to the list of level sets used for output. If more than one are specified,...
Definition lsToVoxelMesh.hpp:70
ToVoxelMesh(const std::vector< SmartPointer< Domain< T, D > > > passedLevelSets, SmartPointer< Mesh< T > > passedMesh)
Definition lsToVoxelMesh.hpp:60
void apply()
Definition lsToVoxelMesh.hpp:80
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
Definition lsAdvect.hpp:36
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68