ViennaLS
Loading...
Searching...
No Matches
lsMakeGeometry.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <cassert>
4
6
7#include <hrleTypes.hpp>
8
9#include <lsConvexHull.hpp>
10#include <lsDomain.hpp>
11#include <lsFromSurfaceMesh.hpp>
12#include <lsGeometries.hpp>
13#include <lsMesh.hpp>
14#include <lsTransformMesh.hpp>
15
16#include <vcVectorType.hpp>
17
18#ifndef NDEBUG
19#include <lsVTKWriter.hpp>
20#endif
21
22namespace viennals {
23
24using namespace viennacore;
25
27template <class T, int D> class MakeGeometry {
28 typedef typename Domain<T, D>::PointValueVectorType pointDataType;
29
32 enum struct GeometryEnum : unsigned {
33 SPHERE = 0,
34 PLANE = 1,
35 BOX = 2,
36 CUSTOM = 3,
37 CYLINDER = 4
38 };
39
40 SmartPointer<Domain<T, D>> levelSet;
41 GeometryEnum geometry = GeometryEnum::SPHERE;
42 SmartPointer<Sphere<T, D>> sphere;
43 SmartPointer<Plane<T, D>> plane;
44 SmartPointer<Box<T, D>> box;
45 SmartPointer<Cylinder<T, D>> cylinder;
46 SmartPointer<PointCloud<T, D>> pointCloud;
47 const double numericEps = 1e-9;
48 // bool ignoreBoundaryConditions = false;
49 std::array<bool, 3> ignoreBoundaryConditions{false, false, false};
50
51public:
52 MakeGeometry() = default;
53
54 MakeGeometry(SmartPointer<Domain<T, D>> passedLevelSet)
55 : levelSet(passedLevelSet) {}
56
57 MakeGeometry(SmartPointer<Domain<T, D>> passedLevelSet,
58 SmartPointer<Sphere<T, D>> passedSphere)
59 : levelSet(passedLevelSet), sphere(passedSphere) {
60 geometry = GeometryEnum::SPHERE;
61 }
62
63 MakeGeometry(SmartPointer<Domain<T, D>> passedLevelSet,
64 SmartPointer<Plane<T, D>> passedPlane)
65 : levelSet(passedLevelSet), plane(passedPlane) {
66 geometry = GeometryEnum::PLANE;
67 }
68
69 MakeGeometry(SmartPointer<Domain<T, D>> passedLevelSet,
70 SmartPointer<Box<T, D>> passedBox)
71 : levelSet(passedLevelSet), box(passedBox) {
72 geometry = GeometryEnum::BOX;
73 }
74
75 MakeGeometry(SmartPointer<Domain<T, D>> passedLevelSet,
76 SmartPointer<Cylinder<T, D>> passedCylinder)
77 : levelSet(passedLevelSet), cylinder(passedCylinder) {
78 geometry = GeometryEnum::CYLINDER;
79 }
80
81 MakeGeometry(SmartPointer<Domain<T, D>> passedLevelSet,
82 SmartPointer<PointCloud<T, D>> passedPointCloud)
83 : levelSet(passedLevelSet), pointCloud(passedPointCloud) {
84 geometry = GeometryEnum::CUSTOM;
85 }
86
87 void setLevelSet(SmartPointer<Domain<T, D>> passedlsDomain) {
88 levelSet = passedlsDomain;
89 }
90
92 void setGeometry(SmartPointer<Sphere<T, D>> passedSphere) {
93 sphere = passedSphere;
94 geometry = GeometryEnum::SPHERE;
95 }
96
98 void setGeometry(SmartPointer<Plane<T, D>> passedPlane) {
99 plane = passedPlane;
100 geometry = GeometryEnum::PLANE;
101 }
102
104 void setGeometry(SmartPointer<Box<T, D>> passedBox) {
105 box = passedBox;
106 geometry = GeometryEnum::BOX;
107 }
108
110 void setGeometry(SmartPointer<Cylinder<T, D>> passedCylinder) {
111 cylinder = passedCylinder;
112 geometry = GeometryEnum::CYLINDER;
113 }
114
117 void setGeometry(SmartPointer<PointCloud<T, D>> passedPointCloud) {
118 pointCloud = passedPointCloud;
119 geometry = GeometryEnum::CUSTOM;
120 }
121
124 void setIgnoreBoundaryConditions(bool passedIgnoreBoundaryConditions) {
125 for (unsigned i = 0; i < D; ++i) {
126 ignoreBoundaryConditions[i] = passedIgnoreBoundaryConditions;
127 }
128 }
129
133 template <std::size_t N>
135 std::array<bool, N> passedIgnoreBoundaryConditions) {
136 for (unsigned i = 0; i < D && i < N; ++i) {
137 ignoreBoundaryConditions[i] = passedIgnoreBoundaryConditions[i];
138 }
139 }
140
141 void apply() {
142 if (levelSet == nullptr) {
143 Logger::getInstance()
144 .addWarning("No level set was passed to MakeGeometry.")
145 .print();
146 return;
147 }
148
149 switch (geometry) {
150 case GeometryEnum::SPHERE:
151 makeSphere(sphere->origin, sphere->radius);
152 break;
153 case GeometryEnum::PLANE:
154 makePlane(plane->origin, plane->normal);
155 break;
156 case GeometryEnum::BOX:
157 makeBox(box->minCorner, box->maxCorner);
158 break;
159 case GeometryEnum::CYLINDER:
160 makeCylinder(cylinder);
161 break;
162 case GeometryEnum::CUSTOM:
163 makeCustom(pointCloud);
164 break;
165 default:
166 Logger::getInstance()
167 .addWarning("Invalid geometry type was specified for MakeGeometry. "
168 "Not creating geometry.")
169 .print();
170 }
171 }
172
173private:
174 void makeSphere(VectorType<T, D> origin, T radius) {
175 if (levelSet == nullptr) {
176 Logger::getInstance()
177 .addWarning("No level set was passed to MakeGeometry.")
178 .print();
179 return;
180 }
181
182 // TODO, this is a stupid algorithm and scales with volume, which is madness
183 auto &grid = levelSet->getGrid();
184 viennahrle::CoordType gridDelta = grid.getGridDelta();
185
186 // calculate indices from sphere size
187 viennahrle::Index<D> index;
188 viennahrle::Index<D> endIndex;
189
190 for (unsigned i = 0; i < D; ++i) {
191 index[i] = (origin[i] - radius) / gridDelta - 1;
192 endIndex[i] = (origin[i] + radius) / gridDelta + 1;
193 }
194
195 constexpr double initialWidth = 2.;
196 const T valueLimit = initialWidth * 0.5 * gridDelta + 1e-5;
197 const T radius2 = radius * radius;
198
199 pointDataType pointData;
200 const viennahrle::Index<D> minIndex = index;
201
202 while (index < endIndex) {
203 // take the shortest manhattan distance to gridline intersection
204 T distance = std::numeric_limits<T>::max();
205 for (unsigned i = 0; i < D; ++i) {
206 T y = (index[(i + 1) % D] * gridDelta) - origin[(i + 1) % D];
207 T z = 0;
208 if constexpr (D == 3)
209 z = (index[(i + 2) % D] * gridDelta) - origin[(i + 2) % D];
210 T x = radius2 - y * y - z * z;
211 if (x < 0.)
212 continue;
213 T dirRadius =
214 std::abs((index[i] * gridDelta) - origin[i]) - std::sqrt(x);
215 if (std::abs(dirRadius) < std::abs(distance))
216 distance = dirRadius;
217 }
218
219 if (std::abs(distance) <= valueLimit) {
220 pointData.push_back(std::make_pair(index, distance / gridDelta));
221 }
222 int dim = 0;
223 for (; dim < D - 1; ++dim) {
224 if (index[dim] < endIndex[dim])
225 break;
226 index[dim] = minIndex[dim];
227 }
228 ++index[dim];
229 }
230
231 // Mirror indices correctly into domain, unless boundary conditions
232 // are ignored
233 for (unsigned i = 0; i < pointData.size(); ++i) {
234 for (unsigned j = 0; j < D; ++j) {
235 if (!ignoreBoundaryConditions[j] && grid.isBoundaryPeriodic(j)) {
236 pointData[i].first[j] =
237 grid.globalIndex2LocalIndex(j, pointData[i].first[j]);
238 }
239 }
240 }
241
242 levelSet->insertPoints(pointData);
243 levelSet->getDomain().segment();
244 levelSet->finalize(initialWidth);
245 }
246
249 void makePlane(VectorType<T, D> origin,
250 VectorType<T, D> const &passedNormal) {
251 if (levelSet == nullptr) {
252 Logger::getInstance()
253 .addWarning("No level set was passed to MakeGeometry.")
254 .print();
255 return;
256 }
257
258 auto &grid = levelSet->getGrid();
259 viennahrle::CoordType gridDelta = grid.getGridDelta();
260
261 // normalise passedNormal
262 double modulus = 0.;
263 VectorType<T, D> normal = passedNormal;
264 for (unsigned i = 0; i < D; ++i) {
265 modulus += normal[i] * normal[i];
266 }
267 modulus = std::sqrt(modulus);
268 for (unsigned i = 0; i < D; ++i) {
269 normal[i] /= modulus;
270 }
271
272 // check that boundary conditions are correct
273 unsigned i = 0;
274 bool infDimSet = false;
275 for (unsigned n = 0; n < D; ++n) {
276 if (grid.getBoundaryConditions(n) ==
277 viennahrle::BoundaryType::INFINITE_BOUNDARY) {
278 if (!infDimSet) {
279 i = n;
280 infDimSet = true;
281 } else {
282 Logger::getInstance().addError(
283 "Planes can only be created with one Infinite Boundary "
284 "Condition. More than one found!");
285 }
286 }
287 }
288 if (!infDimSet) {
289 Logger::getInstance().addError("Planes require exactly one Infinite "
290 "Boundary Condition. None found!");
291 }
292
293 if (passedNormal[i] == 0.) {
294 Logger::getInstance().addError(
295 "MakeGeometry: Plane cannot be parallel to Infinite Boundary "
296 "direction!");
297 }
298
299 // find minimum and maximum points in infinite direction
300 // there are 2*(D-1) points in the corners of the simulation domain
301 std::vector<Vec3D<T>> cornerPoints;
302 cornerPoints.resize(2 * (D - 1));
303
304 // cyclic permutations
305 unsigned j = (i + 1) % D;
306 unsigned k = (i + 2) % D;
307
308 double minCoord[2];
309 double maxCoord[2];
310 // Find grid boundaries, there used to be a +-1 for the coords.
311 // If an error pops up here, probably has to do with that.
312 // But if +-1 is added here, the boundaries are exceeded and
313 // the correct boundary conditions will add stray points for
314 // tilted planes in lsFromSurfaceMesh later on.
315 for (unsigned n = 0; n < D - 1; ++n) {
316 minCoord[n] = gridDelta * (grid.getMinIndex((i + n + 1) % D) - 1);
317 maxCoord[n] = gridDelta * (grid.getMaxIndex((i + n + 1) % D) + 1);
318 }
319
320 // set corner points
321 cornerPoints[0][j] = minCoord[0];
322 cornerPoints[1][j] = maxCoord[0];
323
324 if constexpr (D == 3) {
325 cornerPoints[0][k] = minCoord[1];
326 cornerPoints[1][k] = maxCoord[1];
327
328 cornerPoints[2][j] = minCoord[0];
329 cornerPoints[2][k] = maxCoord[1];
330 cornerPoints[3][j] = maxCoord[0];
331 cornerPoints[3][k] = minCoord[1];
332 }
333
334 // now find i coordinate of points
335 auto mesh = SmartPointer<Mesh<T>>::New();
336
337 for (unsigned n = 0; n < cornerPoints.size(); ++n) {
338 double numerator = (cornerPoints[n][j] - origin[j]) * normal[j];
339 if constexpr (D == 3)
340 numerator += (cornerPoints[n][k] - origin[k]) * normal[k];
341 else
342 cornerPoints[n][2] = 0.;
343 cornerPoints[n][i] = origin[i] - numerator / normal[i];
344 mesh->insertNextNode(cornerPoints[n]);
345 }
346
347 if (D == 2) {
348 std::array<unsigned, 2> line = {0, 1};
349 if (normal[i] < 0.)
350 std::swap(line[0], line[1]);
351 mesh->insertNextLine(line);
352 } else {
353 std::array<unsigned, 3> triangle = {0, 1, 2};
354 if (normal[i] < 0.)
355 std::swap(triangle[0], triangle[1]);
356 mesh->insertNextTriangle(triangle);
357 triangle = {0, 3, 1};
358 if (normal[i] < 0.)
359 std::swap(triangle[0], triangle[1]);
360 mesh->insertNextTriangle(triangle);
361 }
362
363#ifndef NDEBUG
364 static unsigned planeCounter = 0;
365 VTKWriter<T>(mesh, "plane" + std::to_string(planeCounter++) + ".vtk")
366 .apply();
367#endif
368
369 // now convert mesh to levelset
370 FromSurfaceMesh<T, D>(levelSet, mesh).apply();
371 }
372
373 // This function creates a box starting in minCorner spanning to maxCorner
374 void makeBox(VectorType<T, D> minCorner, VectorType<T, D> maxCorner) {
375 if (levelSet == nullptr) {
376 Logger::getInstance()
377 .addWarning("No level set was passed to MakeGeometry.")
378 .print();
379 return;
380 }
381
382 // draw all triangles for the surface and then import from the mesh
383 std::vector<Vec3D<T>> corners;
384 corners.resize(std::pow(2, D), Vec3D<T>{0, 0, 0});
385
386 // first corner is the minCorner
387 for (unsigned i = 0; i < D; ++i)
388 corners[0][i] = minCorner[i];
389
390 // last corner is maxCorner
391 for (unsigned i = 0; i < D; ++i)
392 corners.back()[i] = maxCorner[i];
393
394 // calculate all missing corners
395 corners[1] = corners[0];
396 corners[1][0] = corners.back()[0];
397
398 corners[2] = corners[0];
399 corners[2][1] = corners.back()[1];
400
401 if constexpr (D == 3) {
402 corners[3] = corners.back();
403 corners[3][2] = corners[0][2];
404
405 corners[4] = corners[0];
406 corners[4][2] = corners.back()[2];
407
408 corners[5] = corners.back();
409 corners[5][1] = corners[0][1];
410
411 corners[6] = corners.back();
412 corners[6][0] = corners[0][0];
413 }
414
415 // now add all corners to mesh
416 auto mesh = SmartPointer<Mesh<T>>::New();
417 for (unsigned i = 0; i < corners.size(); ++i) {
418 mesh->insertNextNode(corners[i]);
419 }
420
421 if (D == 2) {
422 std::array<unsigned, 2> lines[4] = {{0, 2}, {2, 3}, {3, 1}, {1, 0}};
423 for (auto &line : lines)
424 mesh->insertNextLine(line);
425 } else {
426 std::array<unsigned, 3> triangles[12] = {
427 {0, 3, 1}, {0, 2, 3}, {0, 1, 5}, {0, 5, 4}, {0, 4, 2}, {4, 6, 2},
428 {7, 6, 4}, {7, 4, 5}, {7, 2, 6}, {7, 3, 2}, {1, 3, 5}, {3, 7, 5}};
429 for (auto &triangle : triangles)
430 mesh->insertNextTriangle(triangle);
431 }
432
433 // now convert mesh to levelset
434 FromSurfaceMesh<T, D> mesher(levelSet, mesh);
435 mesher.setRemoveBoundaryTriangles(ignoreBoundaryConditions);
436 mesher.apply();
437 }
438
439 void makeCylinder(SmartPointer<Cylinder<T, D>> cylinder) {
440 if (D != 3) {
441 Logger::getInstance()
442 .addWarning("MakeGeometry: Cylinder can only be created in 3D!")
443 .print();
444 return;
445 }
446 // generate the points on the edges of the cylinders and mesh
447 // them manually
448 // cylinder axis will be (0,0,1)
449 auto gridDelta = levelSet->getGrid().getGridDelta();
450
451 auto points = SmartPointer<PointCloud<T, D>>::New();
452 const unsigned numPoints =
453 std::ceil(2 * M_PI * cylinder->radius / gridDelta);
454 const double smallAngle = 2.0 * M_PI / static_cast<double>(numPoints);
455
456 auto mesh = SmartPointer<Mesh<T>>::New();
457 // insert midpoint at base
458 mesh->insertNextNode(Vec3D<T>{0.0, 0.0, 0.0});
459 {
460 constexpr double limit = 2 * M_PI - 1e-6;
461 std::vector<Vec3D<T>> points;
462 if (cylinder->topRadius)
463 std::vector<Vec3D<T>> pointsTop;
464
465 // create and insert points at base
466 for (double angle = 0.; angle < limit; angle += smallAngle) {
467 Vec3D<T> point;
468 point[0] = cylinder->radius * std::cos(angle);
469 point[1] = cylinder->radius * std::sin(angle);
470 point[2] = 0.0;
471 points.push_back(point);
472 mesh->insertNextNode(point);
473 }
474
475 // insert midpoint at top
476 mesh->insertNextNode(Vec3D<T>{0.0, 0.0, cylinder->height});
477
478 double angle = 0;
479 for (unsigned i = 0; i < numPoints; ++i) {
480 // create triangles at base
481 std::array<unsigned, 3> triangle{};
482 triangle[0] = (i + 1) % numPoints + 1;
483 triangle[1] = i + 1;
484 triangle[2] = 0;
485 mesh->insertNextTriangle(triangle);
486
487 // insert points at top
488 // If topRadius is specified, update the first two coordinates of the
489 // points
490 if (cylinder->topRadius) {
491 points[i][0] = cylinder->topRadius * std::cos(angle);
492 points[i][1] = cylinder->topRadius * std::sin(angle);
493 angle += smallAngle;
494 }
495 points[i][2] = cylinder->height;
496 mesh->insertNextNode(points[i]);
497
498 // insert triangles at top
499 triangle[0] = numPoints + 1;
500 triangle[1] = numPoints + i + 2;
501 triangle[2] = (i + 1) % numPoints + 2 + numPoints;
502 mesh->insertNextTriangle(triangle);
503 }
504
505 // insert sidewall triangles
506 for (unsigned i = 0; i < numPoints; ++i) {
507 std::array<unsigned, 3> triangle{};
508 triangle[0] = i + 1;
509 triangle[1] = (i + 1) % numPoints + 1;
510 triangle[2] = i + numPoints + 2;
511 mesh->insertNextTriangle(triangle);
512
513 triangle[0] = (i + 1) % numPoints + 1;
514 triangle[1] = (i + 1) % numPoints + 2 + numPoints;
515 triangle[2] = i + numPoints + 2;
516 mesh->insertNextTriangle(triangle);
517 }
518 }
519
520 // rotate mesh
521 // normalise axis vector
522 T unit =
523 std::sqrt(DotProduct(cylinder->axisDirection, cylinder->axisDirection));
524 Vec3D<T> cylinderAxis;
525 for (int i = 0; i < 3; ++i) {
526 cylinderAxis[i] = cylinder->axisDirection[i] / unit;
527 }
528 // get rotation axis via cross product of (0,0,1) and axis of cylinder
529 Vec3D<T> rotAxis = {-cylinderAxis[1], cylinderAxis[0], 0.0};
530 // angle is acos of dot product
531 T rotationAngle = std::acos(cylinderAxis[2]);
532
533 // rotate mesh
534 TransformMesh<T>(mesh, TransformEnum::ROTATION, rotAxis, rotationAngle)
535 .apply();
536
537 // translate mesh
538 Vec3D<T> translationVector;
539 for (int i = 0; i < 3; ++i) {
540 translationVector[i] = cylinder->origin[i];
541 }
542 TransformMesh<T>(mesh, TransformEnum::TRANSLATION, translationVector)
543 .apply();
544
545 // read mesh from surface
546 FromSurfaceMesh<T, D> mesher(levelSet, mesh);
547 mesher.setRemoveBoundaryTriangles(ignoreBoundaryConditions);
548 mesher.apply();
549 }
550
551 void makeCustom(SmartPointer<PointCloud<T, D>> pointCloud) {
552 // create mesh from point cloud
553 auto mesh = SmartPointer<Mesh<T>>::New();
554 ConvexHull<T, D>(mesh, pointCloud).apply();
555
556 // read mesh from surface
557 FromSurfaceMesh<T, D> mesher(levelSet, mesh);
558 mesher.setRemoveBoundaryTriangles(ignoreBoundaryConditions);
559 mesher.apply();
560 }
561};
562
563// add all template specialisations for this class
565
566} // namespace viennals
Class describing a square box from one coordinate to another.
Definition lsGeometries.hpp:62
Class describing a square box from one coordinate to another.
Definition lsGeometries.hpp:87
Class containing all information about the level set, including the dimensions of the domain,...
Definition lsDomain.hpp:27
std::vector< std::pair< viennahrle::Index< D >, T > > PointValueVectorType
Definition lsDomain.hpp:34
Create level sets describing basic geometric forms.
Definition lsMakeGeometry.hpp:27
MakeGeometry(SmartPointer< Domain< T, D > > passedLevelSet, SmartPointer< Cylinder< T, D > > passedCylinder)
Definition lsMakeGeometry.hpp:75
void setIgnoreBoundaryConditions(bool passedIgnoreBoundaryConditions)
Ignore boundary conditions, meaning the parts of the generated geometry which are outside of the doma...
Definition lsMakeGeometry.hpp:124
void apply()
Definition lsMakeGeometry.hpp:141
void setLevelSet(SmartPointer< Domain< T, D > > passedlsDomain)
Definition lsMakeGeometry.hpp:87
void setGeometry(SmartPointer< Cylinder< T, D > > passedCylinder)
Set a cylinder to be created in the level set.
Definition lsMakeGeometry.hpp:110
void setGeometry(SmartPointer< PointCloud< T, D > > passedPointCloud)
Set a point cloud, which is used to create a geometry from its convex hull.
Definition lsMakeGeometry.hpp:117
void setIgnoreBoundaryConditions(std::array< bool, N > passedIgnoreBoundaryConditions)
Ignore boundary conditions, meaning the parts of the generated geometry which are outside of the doma...
Definition lsMakeGeometry.hpp:134
MakeGeometry(SmartPointer< Domain< T, D > > passedLevelSet, SmartPointer< Sphere< T, D > > passedSphere)
Definition lsMakeGeometry.hpp:57
MakeGeometry(SmartPointer< Domain< T, D > > passedLevelSet)
Definition lsMakeGeometry.hpp:54
MakeGeometry(SmartPointer< Domain< T, D > > passedLevelSet, SmartPointer< Plane< T, D > > passedPlane)
Definition lsMakeGeometry.hpp:63
MakeGeometry(SmartPointer< Domain< T, D > > passedLevelSet, SmartPointer< Box< T, D > > passedBox)
Definition lsMakeGeometry.hpp:69
void setGeometry(SmartPointer< Box< T, D > > passedBox)
Set a box to be created in the level set.
Definition lsMakeGeometry.hpp:104
void setGeometry(SmartPointer< Plane< T, D > > passedPlane)
Set a plane to be created in the level set.
Definition lsMakeGeometry.hpp:98
void setGeometry(SmartPointer< Sphere< T, D > > passedSphere)
Set sphere as geometry to be created in the level set.
Definition lsMakeGeometry.hpp:92
MakeGeometry(SmartPointer< Domain< T, D > > passedLevelSet, SmartPointer< PointCloud< T, D > > passedPointCloud)
Definition lsMakeGeometry.hpp:81
Class describing a plane via a point in it and the plane normal.
Definition lsGeometries.hpp:37
Class describing a point cloud, which can be used to create geometries from its convex hull mesh.
Definition lsGeometries.hpp:128
Class describing a sphere via origin and radius.
Definition lsGeometries.hpp:14
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
float gridDelta
Definition AirGapDeposition.py:21
tuple maxCorner
Definition AirGapDeposition.py:44
mesh
Definition AirGapDeposition.py:36
tuple origin
Definition AirGapDeposition.py:30
Definition lsAdvect.hpp:36
@ TRANSLATION
Definition lsTransformMesh.hpp:16
@ ROTATION
Definition lsTransformMesh.hpp:17
constexpr int D
Definition pyWrap.cpp:71
double T
Definition pyWrap.cpp:69