ViennaLS
Loading...
Searching...
No Matches
lsCompareArea.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <limits>
4#include <unordered_map>
5
6#include <hrleDenseCellIterator.hpp>
7#include <lsDomain.hpp>
8#include <lsMesh.hpp>
10
11namespace viennals {
12
13using namespace viennacore;
14
23template <class T, int D = 2> class CompareArea {
24 using hrleDomainType = typename Domain<T, D>::DomainType;
25 using hrleIndexType = viennahrle::IndexType;
26
27 SmartPointer<Domain<T, D>> levelSetTarget = nullptr;
28 SmartPointer<Domain<T, D>> levelSetSample = nullptr;
29 viennahrle::Index<D> minIndex, maxIndex;
30
31 unsigned long int differentCellsCount = 0;
32 unsigned long int customDifferentCellCount = 0;
33
34 hrleIndexType xRangeMin = std::numeric_limits<hrleIndexType>::lowest();
35 hrleIndexType xRangeMax = std::numeric_limits<hrleIndexType>::max();
36 hrleIndexType yRangeMin = std::numeric_limits<hrleIndexType>::lowest();
37 hrleIndexType yRangeMax = std::numeric_limits<hrleIndexType>::max();
38 bool useCustomXIncrement = false;
39 bool useCustomYIncrement = false;
40
41 unsigned short int customXIncrement = 0;
42 unsigned short int customYIncrement = 0;
43 unsigned short int defaultIncrement = 1;
44
45 double gridDelta = 0.0;
46
47 // Mesh output related members
48 SmartPointer<Mesh<T>> outputMesh = nullptr;
49
50 bool checkAndCalculateBounds() {
51 if (levelSetTarget == nullptr || levelSetSample == nullptr) {
52 Logger::getInstance()
53 .addError("Missing level set in CompareArea.")
54 .print();
55 return false;
56 }
57
58 // Check if the grids are compatible
59 const auto &gridTarget = levelSetTarget->getGrid();
60 const auto &gridSample = levelSetSample->getGrid();
61
62 if (gridTarget.getGridDelta() != gridSample.getGridDelta()) {
63 Logger::getInstance()
64 .addError("Grid delta mismatch in CompareArea. The grid deltas of "
65 "the two level sets must be equal.")
66 .print();
67 return false;
68 } else {
69 gridDelta = gridTarget.getGridDelta();
70 }
71
72 // Initialize min and max indices
73 for (unsigned i = 0; i < D; ++i) {
74 minIndex[i] = std::numeric_limits<hrleIndexType>::max();
75 maxIndex[i] = std::numeric_limits<hrleIndexType>::lowest();
76 }
77
78 // Get the boundaries of the two level sets
79 auto &domainTarget = levelSetTarget->getDomain();
80 auto &domainSample = levelSetSample->getDomain();
81
82 // Calculate actual bounds
83 for (unsigned i = 0; i < D; ++i) {
84 minIndex[i] = std::min({minIndex[i],
85 (gridTarget.isNegBoundaryInfinite(i))
86 ? domainTarget.getMinRunBreak(i)
87 : gridTarget.getMinBounds(i),
88 (gridSample.isNegBoundaryInfinite(i))
89 ? domainSample.getMinRunBreak(i)
90 : gridSample.getMinBounds(i)});
91
92 maxIndex[i] = std::max({maxIndex[i],
93 (gridTarget.isPosBoundaryInfinite(i))
94 ? domainTarget.getMaxRunBreak(i)
95 : gridTarget.getMaxBounds(i),
96 (gridSample.isPosBoundaryInfinite(i))
97 ? domainSample.getMaxRunBreak(i)
98 : gridSample.getMaxBounds(i)});
99 }
100
101 return true;
102 }
103
104public:
106 static_assert(
107 D == 2 &&
108 "CompareArea is currently only implemented for 2D level sets.");
109 }
110
111 CompareArea(SmartPointer<Domain<T, D>> passedLevelSetTarget,
112 SmartPointer<Domain<T, D>> passedLevelSetSample)
113 : levelSetTarget(passedLevelSetTarget),
114 levelSetSample(passedLevelSetSample) {
115 static_assert(
116 D == 2 &&
117 "CompareArea is currently only implemented for 2D level sets.");
118 }
119
121 void setLevelSetTarget(SmartPointer<Domain<T, D>> passedLevelSet) {
122 levelSetTarget = passedLevelSet;
123 }
124
126 void setLevelSetSample(SmartPointer<Domain<T, D>> passedLevelSet) {
127 levelSetSample = passedLevelSet;
128 }
129
131 void setDefaultIncrement(unsigned short int increment) {
132 defaultIncrement = increment;
133 }
134
136 void setXRangeAndIncrement(hrleIndexType minXRange, hrleIndexType maxXRange,
137 unsigned short int Xincrement) {
138 xRangeMin = minXRange;
139 xRangeMax = maxXRange;
140 customXIncrement = Xincrement;
141 useCustomXIncrement = true;
142 }
143
145 void setYRangeAndIncrement(hrleIndexType minYRange, hrleIndexType maxYRange,
146 unsigned short int Yincrement) {
147 yRangeMin = minYRange;
148 yRangeMax = maxYRange;
149 customYIncrement = Yincrement;
150 useCustomYIncrement = true;
151 }
152
157 void setOutputMesh(SmartPointer<Mesh<T>> passedMesh) {
158 outputMesh = passedMesh;
159 }
160
162 double getAreaMismatch() const {
163 return static_cast<double>(differentCellsCount) * gridDelta * gridDelta;
164 }
165
167 double getCustomAreaMismatch() const {
168 return static_cast<double>(customDifferentCellCount) * gridDelta *
169 gridDelta;
170 }
171
173 unsigned long int getCellCount() const { return differentCellsCount; }
174
177 unsigned long int getCustomCellCount() const {
178 return customDifferentCellCount;
179 }
180
182 void apply() {
183 // Calculate the bounds for iteration
184 if (!checkAndCalculateBounds()) {
185 differentCellsCount =
186 std::numeric_limits<unsigned long int>::max(); // Error indicator
187 customDifferentCellCount =
188 std::numeric_limits<unsigned long int>::max(); // Error indicator
189 return;
190 }
191
192 // Set up dense cell iterators for both level sets
193 viennahrle::ConstDenseCellIterator<hrleDomainType> itTarget(
194 levelSetTarget->getDomain(), minIndex);
195 viennahrle::ConstDenseCellIterator<hrleDomainType> itSample(
196 levelSetSample->getDomain(), minIndex);
197
198 differentCellsCount = 0;
199 customDifferentCellCount = 0;
200
201 // Initialize mesh-related variables if generating a mesh
202 const bool generateMesh = outputMesh != nullptr;
203 if (generateMesh) {
204 // Save the extent of the resulting mesh
205 outputMesh->clear();
206 for (unsigned i = 0; i < D; ++i) {
207 outputMesh->minimumExtent[i] = std::numeric_limits<T>::max();
208 outputMesh->maximumExtent[i] = std::numeric_limits<T>::lowest();
209 }
210 }
211
212 // Vector for cell differences and increment values to be inserted in the
213 // mesh
214 std::vector<T> cellDifference;
215 std::vector<T> incrementValues;
216 size_t currentPointId = 0;
217 std::unordered_map<viennahrle::Index<D>, size_t,
218 typename viennahrle::Index<D>::hash>
219 pointIdMapping;
220
221 // Iterate through the domain defined by the bounding box
222 for (; itTarget.getIndices() < maxIndex; itTarget.next()) {
223 itSample.goToIndicesSequential(itTarget.getIndices());
224
225 // Compare cell values at the current position
226 T centerValueTarget = 0.;
227 T centerValueSample = 0.;
228
229 // Calculate the average value of the corners for both level sets
230 for (int i = 0; i < (1 << D); ++i) {
231 centerValueTarget += itTarget.getCorner(i).getValue();
232 centerValueSample += itSample.getCorner(i).getValue();
233 }
234 centerValueTarget /= (1 << D);
235 centerValueSample /= (1 << D);
236
237 // Determine if cell is inside for each level set
238 bool insideTarget = (centerValueTarget <= 0.);
239 bool insideSample = (centerValueSample <= 0.);
240 bool isDifferent = insideTarget != insideSample;
241
242 // Calculate range checks once per cell
243 bool inXRange = useCustomXIncrement &&
244 (itTarget.getIndices()[0] * gridDelta >= xRangeMin &&
245 itTarget.getIndices()[0] * gridDelta <= xRangeMax);
246 bool inYRange = useCustomYIncrement &&
247 (itTarget.getIndices()[1] * gridDelta >= yRangeMin &&
248 itTarget.getIndices()[1] * gridDelta <= yRangeMax);
249
250 // Calculate increment to add based on ranges
251 unsigned short int incrementToAdd = defaultIncrement;
252 if (inXRange && inYRange) {
253 // Apply both increments
254 incrementToAdd = customXIncrement + customYIncrement;
255 } else if (inXRange) {
256 incrementToAdd = customXIncrement;
257 } else if (inYRange) {
258 incrementToAdd = customYIncrement;
259 }
260
261 // If cells differ, update the counters
262 if (isDifferent) {
263 // Always increment simple cell count by 1
264 differentCellsCount += 1;
265
266 // For custom cell count, apply the calculated increment
267 customDifferentCellCount += incrementToAdd;
268 }
269
270 // For mesh generation, process all cells where at least one level set is
271 // inside
272 if (generateMesh && (insideTarget || insideSample)) {
273 // Material ID: 0 for matching inside, 1 for mismatched
274 int difference = isDifferent ? 1 : 0;
275
276 std::array<unsigned, 1 << D> voxel;
277 bool addVoxel = true;
278 // TODO: I think this check whether Voxel is added or not can be removed
279
280 // Insert all points of voxel into pointList
281 for (unsigned i = 0; i < (1 << D); ++i) {
282 viennahrle::Index<D> index;
283 for (unsigned j = 0; j < D; ++j) {
284 index[j] =
285 itTarget.getIndices(j) + itTarget.getCorner(i).getOffset()[j];
286 if (index[j] > maxIndex[j]) {
287 addVoxel = false;
288 break;
289 }
290 }
291 if (addVoxel) {
292 auto pointIdValue = std::make_pair(index, currentPointId);
293 auto pointIdPair = pointIdMapping.insert(pointIdValue);
294 voxel[i] = pointIdPair.first->second;
295 if (pointIdPair.second) {
296 ++currentPointId;
297 }
298 } else {
299 break;
300 }
301 }
302
303 if (addVoxel) {
304 if constexpr (D == 3) {
305 std::array<unsigned, 8> hexa{voxel[0], voxel[1], voxel[3],
306 voxel[2], voxel[4], voxel[5],
307 voxel[7], voxel[6]};
308 outputMesh->hexas.push_back(hexa);
309 } else if constexpr (D == 2) {
310 std::array<unsigned, 4> quad{voxel[0], voxel[1], voxel[3],
311 voxel[2]};
312 outputMesh->tetras.push_back(quad);
313 }
314
315 cellDifference.push_back(difference);
316
317 // Store increment value (0 if not different)
318 incrementValues.push_back(isDifferent ? static_cast<T>(incrementToAdd)
319 : 0);
320 }
321 }
322 }
323
324 // Finalize mesh if generating one
325 if (generateMesh && !pointIdMapping.empty()) {
326 // Insert points into the mesh
327 outputMesh->nodes.resize(pointIdMapping.size());
328 for (auto it = pointIdMapping.begin(); it != pointIdMapping.end(); ++it) {
329 Vec3D<T> coords;
330 for (unsigned i = 0; i < D; ++i) {
331 coords[i] = gridDelta * it->first[i];
332
333 if (coords[i] < outputMesh->minimumExtent[i]) {
334 outputMesh->minimumExtent[i] = coords[i];
335 } else if (coords[i] > outputMesh->maximumExtent[i]) {
336 outputMesh->maximumExtent[i] = coords[i];
337 }
338 }
339 outputMesh->nodes[it->second] = coords;
340 }
341
342 // Add cell data to the mesh
343 assert(cellDifference.size() == incrementValues.size());
344 assert(incrementValues.size() ==
345 outputMesh->template getElements<1 << D>().size());
346 outputMesh->cellData.insertNextScalarData(cellDifference, "Difference");
347 outputMesh->cellData.insertNextScalarData(incrementValues,
348 "CustomIncrement");
349 }
350 }
351};
352
353// Precompile for common precision and dimensions
355
356} // namespace viennals
CompareArea(SmartPointer< Domain< T, D > > passedLevelSetTarget, SmartPointer< Domain< T, D > > passedLevelSetSample)
Definition lsCompareArea.hpp:111
unsigned long int getCustomCellCount() const
Returns the number of cells where the level sets differ, with custom increments applied.
Definition lsCompareArea.hpp:177
void setXRangeAndIncrement(hrleIndexType minXRange, hrleIndexType maxXRange, unsigned short int Xincrement)
Sets the x-range and custom increment value.
Definition lsCompareArea.hpp:136
void setLevelSetTarget(SmartPointer< Domain< T, D > > passedLevelSet)
Sets the target level set.
Definition lsCompareArea.hpp:121
unsigned long int getCellCount() const
Returns the number of cells where the level sets differ.
Definition lsCompareArea.hpp:173
double getCustomAreaMismatch() const
Returns the computed area mismatch, with custom increments applied.
Definition lsCompareArea.hpp:167
void setLevelSetSample(SmartPointer< Domain< T, D > > passedLevelSet)
Sets the sample level set.
Definition lsCompareArea.hpp:126
void apply()
Computes the area difference between the two level sets.
Definition lsCompareArea.hpp:182
void setOutputMesh(SmartPointer< Mesh< T > > passedMesh)
Set the output mesh where difference areas will be stored for visualization. Each cell in the mesh wi...
Definition lsCompareArea.hpp:157
void setDefaultIncrement(unsigned short int increment)
Set default increment value.
Definition lsCompareArea.hpp:131
void setYRangeAndIncrement(hrleIndexType minYRange, hrleIndexType maxYRange, unsigned short int Yincrement)
Sets the y-range and custom increment value.
Definition lsCompareArea.hpp:145
double getAreaMismatch() const
Returns the computed area mismatch.
Definition lsCompareArea.hpp:162
CompareArea()
Definition lsCompareArea.hpp:105
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
#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