ViennaLS
Loading...
Searching...
No Matches
lsFromMesh.hpp
Go to the documentation of this file.
1#pragma once
2
4
5#include <lsDomain.hpp>
6#include <lsMesh.hpp>
7
8namespace viennals {
9
10using namespace viennacore;
11
16template <class T, int D> class FromMesh {
17 typedef typename Domain<T, D>::DomainType hrleDomainType;
18
19 SmartPointer<Domain<T, D>> levelSet = nullptr;
20 SmartPointer<Mesh<T>> mesh = nullptr;
21 bool sortPointList = true;
22
23public:
24 FromMesh() = default;
25
26 FromMesh(SmartPointer<Domain<T, D>> passedLevelSet,
27 SmartPointer<Mesh<T>> passedMesh)
28 : levelSet(passedLevelSet), mesh(passedMesh) {}
29
30 void setLevelSet(SmartPointer<Domain<T, D>> passedlsDomain) {
31 levelSet = passedlsDomain;
32 }
33
34 void setMesh(const SmartPointer<Mesh<T>> passedMesh) { mesh = passedMesh; }
35
36 void setSortPointList(bool passedSortPointList) {
37 sortPointList = passedSortPointList;
38 }
39
40 void apply() {
41 if (levelSet == nullptr) {
42 Logger::getInstance()
43 .addWarning("No level set was passed to FromMesh.")
44 .print();
45 return;
46 }
47 if (mesh == nullptr) {
48 Logger::getInstance()
49 .addWarning("No Mesh<> was supplied to FromMesh.")
50 .print();
51 return;
52 }
53
54 auto &domain = levelSet->getDomain();
55 auto &nodes = mesh->getNodes();
56 auto values = mesh->cellData.getScalarData("LSValues", true);
57
58 if (values == nullptr) {
59 Logger::getInstance()
60 .addWarning("Mesh does not contain level set values (\"LSValues\").")
61 .print();
62 return;
63 }
64
65 domain.initialize();
66
67 // if there are no points, just initialize an empty hrleDomain
68 if (nodes.empty()) {
69 return;
70 }
71
72 const viennahrle::Grid<D> &grid = domain.getGrid();
73 const T gridDelta = grid.getGridDelta();
74
75 viennahrle::Index<D> lastIndex;
76 for (unsigned i = 0; i < D; ++i) {
77 lastIndex[i] = nodes.front()[i] / gridDelta;
78 }
79
80 if (lastIndex != grid.getMinGridPoint()) {
81 domain.insertNextUndefinedPoint(0, grid.getMinGridPoint(),
82 (values->front() < 0)
85 }
86
87 auto pointDataIt = nodes.begin();
88 auto pointDataEnd = nodes.end();
89
90 auto valueIt = values->begin();
91 auto valueEnd = values->end();
92
93 VectorType<bool, D> signs;
94 std::fill(signs.begin(), signs.end(),
95 values->front() <= -std::numeric_limits<T>::epsilon());
96
97 while (pointDataIt != pointDataEnd && valueIt != valueEnd) {
98 // only read in points within the first 5 layers, to ignore
99 // undefined points
100 if (std::abs(*valueIt) > 2.5) {
101 ++pointDataIt;
102 ++valueIt;
103 continue;
104 }
105
106 bool setPoint = true;
107
108 viennahrle::Index<D> currentIndex;
109 for (unsigned i = 0; i < D; ++i) {
110 currentIndex[i] = std::round(pointDataIt->at(i) / gridDelta);
111 }
112
113 // if boundary conditions are infinite always set the point
114 // if not, check, whether it is inside of domain
115 for (unsigned i = 0; i < D; ++i) {
116 if (grid.getBoundaryConditions(i) !=
117 viennahrle::BoundaryType::INFINITE_BOUNDARY) {
118 if (currentIndex[i] > grid.getMaxBounds(i) ||
119 currentIndex[i] < grid.getMinBounds(i)) {
120 setPoint = false;
121 }
122 }
123 }
124
125 if (setPoint) {
126 // Add defined point as it appears in the list
127 domain.insertNextDefinedPoint(0, currentIndex, *valueIt);
128
129 // determine signs for next undefined runs
130 {
131 bool changeSign = false;
132 for (int i = D - 1; i >= 0; --i) {
133 changeSign = changeSign || (currentIndex[i] > lastIndex[i]);
134 if (changeSign) {
135 signs[i] = *valueIt <= -std::numeric_limits<T>::epsilon();
136 lastIndex[i] = currentIndex[i];
137 }
138 }
139 }
140 }
141
142 viennahrle::Index<D> nextIndex;
143
144 ++pointDataIt;
145 ++valueIt;
146
147 // choose correct next index
148 if (pointDataIt == pointDataEnd) {
149 nextIndex = grid.getMaxGridPoint();
150 ++nextIndex[D - 1];
151 } else {
152 for (unsigned i = 0; i < D; ++i) {
153 nextIndex[i] = std::round(pointDataIt->at(i) / gridDelta);
154 }
155 }
156
157 // move current index by one grid spacing and see if the next
158 // point has the same index, if not, there must be an undefined
159 // run inbetween
160 for (int q = 0; q < D; q++) {
161 viennahrle::Index<D> tmp = currentIndex;
162 ++tmp[q];
163 if (tmp[q] > grid.getMaxGridPoint(q))
164 continue;
165 for (int r = 0; r < q; ++r)
166 tmp[r] = grid.getMinGridPoint(r);
167
168 if (tmp >= nextIndex)
169 break;
170
171 domain.insertNextUndefinedPoint(0, tmp,
172 signs[q] ? Domain<T, D>::NEG_VALUE
174 }
175 }
176
177 domain.finalize();
178 }
179};
180
181// add all template specialisations for this class
183
184} // 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
static constexpr T NEG_VALUE
Definition lsDomain.hpp:52
static constexpr T POS_VALUE
Definition lsDomain.hpp:51
void setMesh(const SmartPointer< Mesh< T > > passedMesh)
Definition lsFromMesh.hpp:34
void setSortPointList(bool passedSortPointList)
Definition lsFromMesh.hpp:36
void setLevelSet(SmartPointer< Domain< T, D > > passedlsDomain)
Definition lsFromMesh.hpp:30
void apply()
Definition lsFromMesh.hpp:40
FromMesh(SmartPointer< Domain< T, D > > passedLevelSet, SmartPointer< Mesh< T > > passedMesh)
Definition lsFromMesh.hpp:26
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