32 enum struct GeometryEnum :
unsigned {
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;
49 std::array<bool, 3> ignoreBoundaryConditions{
false,
false,
false};
55 : levelSet(passedLevelSet) {}
59 : levelSet(passedLevelSet), sphere(passedSphere) {
60 geometry = GeometryEnum::SPHERE;
65 : levelSet(passedLevelSet), plane(passedPlane) {
66 geometry = GeometryEnum::PLANE;
71 : levelSet(passedLevelSet), box(passedBox) {
72 geometry = GeometryEnum::BOX;
77 : levelSet(passedLevelSet), cylinder(passedCylinder) {
78 geometry = GeometryEnum::CYLINDER;
83 : levelSet(passedLevelSet), pointCloud(passedPointCloud) {
84 geometry = GeometryEnum::CUSTOM;
88 levelSet = passedlsDomain;
93 sphere = passedSphere;
94 geometry = GeometryEnum::SPHERE;
100 geometry = GeometryEnum::PLANE;
106 geometry = GeometryEnum::BOX;
111 cylinder = passedCylinder;
112 geometry = GeometryEnum::CYLINDER;
118 pointCloud = passedPointCloud;
119 geometry = GeometryEnum::CUSTOM;
125 for (
unsigned i = 0; i <
D; ++i) {
126 ignoreBoundaryConditions[i] = passedIgnoreBoundaryConditions;
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];
142 if (levelSet ==
nullptr) {
143 Logger::getInstance()
144 .addWarning(
"No level set was passed to MakeGeometry.")
150 case GeometryEnum::SPHERE:
151 makeSphere(sphere->origin, sphere->radius);
153 case GeometryEnum::PLANE:
154 makePlane(plane->origin, plane->normal);
156 case GeometryEnum::BOX:
157 makeBox(box->minCorner, box->maxCorner);
159 case GeometryEnum::CYLINDER:
160 makeCylinder(cylinder);
162 case GeometryEnum::CUSTOM:
163 makeCustom(pointCloud);
166 Logger::getInstance()
167 .addWarning(
"Invalid geometry type was specified for MakeGeometry. "
168 "Not creating geometry.")
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.")
183 auto &grid = levelSet->getGrid();
184 viennahrle::CoordType gridDelta = grid.getGridDelta();
187 viennahrle::Index<D> index;
188 viennahrle::Index<D> endIndex;
190 for (
unsigned i = 0; i <
D; ++i) {
191 index[i] = (origin[i] - radius) / gridDelta - 1;
192 endIndex[i] = (origin[i] + radius) / gridDelta + 1;
195 constexpr double initialWidth = 2.;
196 const T valueLimit = initialWidth * 0.5 *
gridDelta + 1e-5;
197 const T radius2 = radius * radius;
199 pointDataType pointData;
200 const viennahrle::Index<D> minIndex = index;
202 while (index < endIndex) {
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];
208 if constexpr (
D == 3)
209 z = (index[(i + 2) %
D] * gridDelta) -
origin[(i + 2) %
D];
210 T x = radius2 - y * y - z * z;
214 std::abs((index[i] * gridDelta) - origin[i]) - std::sqrt(x);
215 if (std::abs(dirRadius) < std::abs(distance))
216 distance = dirRadius;
219 if (std::abs(distance) <= valueLimit) {
220 pointData.push_back(std::make_pair(index, distance / gridDelta));
223 for (; dim <
D - 1; ++dim) {
224 if (index[dim] < endIndex[dim])
226 index[dim] = minIndex[dim];
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]);
242 levelSet->insertPoints(pointData);
243 levelSet->getDomain().segment();
244 levelSet->finalize(initialWidth);
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.")
258 auto &grid = levelSet->getGrid();
259 viennahrle::CoordType
gridDelta = grid.getGridDelta();
263 VectorType<T, D> normal = passedNormal;
264 for (
unsigned i = 0; i <
D; ++i) {
265 modulus += normal[i] * normal[i];
267 modulus = std::sqrt(modulus);
268 for (
unsigned i = 0; i <
D; ++i) {
269 normal[i] /= modulus;
274 bool infDimSet =
false;
275 for (
unsigned n = 0; n <
D; ++n) {
276 if (grid.getBoundaryConditions(n) ==
277 viennahrle::BoundaryType::INFINITE_BOUNDARY) {
282 Logger::getInstance().addError(
283 "Planes can only be created with one Infinite Boundary "
284 "Condition. More than one found!");
289 Logger::getInstance().addError(
"Planes require exactly one Infinite "
290 "Boundary Condition. None found!");
293 if (passedNormal[i] == 0.) {
294 Logger::getInstance().addError(
295 "MakeGeometry: Plane cannot be parallel to Infinite Boundary "
301 std::vector<Vec3D<T>> cornerPoints;
302 cornerPoints.resize(2 * (
D - 1));
305 unsigned j = (i + 1) %
D;
306 unsigned k = (i + 2) %
D;
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);
321 cornerPoints[0][j] = minCoord[0];
322 cornerPoints[1][j] = maxCoord[0];
324 if constexpr (
D == 3) {
325 cornerPoints[0][k] = minCoord[1];
326 cornerPoints[1][k] = maxCoord[1];
328 cornerPoints[2][j] = minCoord[0];
329 cornerPoints[2][k] = maxCoord[1];
330 cornerPoints[3][j] = maxCoord[0];
331 cornerPoints[3][k] = minCoord[1];
335 auto mesh = SmartPointer<Mesh<T>>::New();
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];
342 cornerPoints[n][2] = 0.;
343 cornerPoints[n][i] =
origin[i] - numerator / normal[i];
344 mesh->insertNextNode(cornerPoints[n]);
348 std::array<unsigned, 2> line = {0, 1};
350 std::swap(line[0], line[1]);
351 mesh->insertNextLine(line);
353 std::array<unsigned, 3> triangle = {0, 1, 2};
355 std::swap(triangle[0], triangle[1]);
356 mesh->insertNextTriangle(triangle);
357 triangle = {0, 3, 1};
359 std::swap(triangle[0], triangle[1]);
360 mesh->insertNextTriangle(triangle);
364 static unsigned planeCounter = 0;
365 VTKWriter<T>(mesh,
"plane" + std::to_string(planeCounter++) +
".vtk")
370 FromSurfaceMesh<T, D>(levelSet, mesh).apply();
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.")
383 std::vector<Vec3D<T>> corners;
384 corners.resize(std::pow(2,
D), Vec3D<T>{0, 0, 0});
387 for (
unsigned i = 0; i <
D; ++i)
388 corners[0][i] = minCorner[i];
391 for (
unsigned i = 0; i <
D; ++i)
395 corners[1] = corners[0];
396 corners[1][0] = corners.back()[0];
398 corners[2] = corners[0];
399 corners[2][1] = corners.back()[1];
401 if constexpr (
D == 3) {
402 corners[3] = corners.back();
403 corners[3][2] = corners[0][2];
405 corners[4] = corners[0];
406 corners[4][2] = corners.back()[2];
408 corners[5] = corners.back();
409 corners[5][1] = corners[0][1];
411 corners[6] = corners.back();
412 corners[6][0] = corners[0][0];
416 auto mesh = SmartPointer<Mesh<T>>::New();
417 for (
unsigned i = 0; i < corners.size(); ++i) {
418 mesh->insertNextNode(corners[i]);
422 std::array<unsigned, 2> lines[4] = {{0, 2}, {2, 3}, {3, 1}, {1, 0}};
423 for (
auto &line : lines)
424 mesh->insertNextLine(line);
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);
434 FromSurfaceMesh<T, D> mesher(levelSet, mesh);
435 mesher.setRemoveBoundaryTriangles(ignoreBoundaryConditions);
439 void makeCylinder(SmartPointer<Cylinder<T, D>> cylinder) {
441 Logger::getInstance()
442 .addWarning(
"MakeGeometry: Cylinder can only be created in 3D!")
449 auto gridDelta = levelSet->getGrid().getGridDelta();
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);
456 auto mesh = SmartPointer<Mesh<T>>::New();
458 mesh->insertNextNode(Vec3D<T>{0.0, 0.0, 0.0});
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;
466 for (
double angle = 0.; angle < limit; angle += smallAngle) {
468 point[0] = cylinder->radius * std::cos(angle);
469 point[1] = cylinder->radius * std::sin(angle);
471 points.push_back(point);
472 mesh->insertNextNode(point);
476 mesh->insertNextNode(Vec3D<T>{0.0, 0.0, cylinder->height});
479 for (
unsigned i = 0; i < numPoints; ++i) {
481 std::array<unsigned, 3> triangle{};
482 triangle[0] = (i + 1) % numPoints + 1;
485 mesh->insertNextTriangle(triangle);
490 if (cylinder->topRadius) {
491 points[i][0] = cylinder->topRadius * std::cos(angle);
492 points[i][1] = cylinder->topRadius * std::sin(angle);
495 points[i][2] = cylinder->height;
496 mesh->insertNextNode(points[i]);
499 triangle[0] = numPoints + 1;
500 triangle[1] = numPoints + i + 2;
501 triangle[2] = (i + 1) % numPoints + 2 + numPoints;
502 mesh->insertNextTriangle(triangle);
506 for (
unsigned i = 0; i < numPoints; ++i) {
507 std::array<unsigned, 3> triangle{};
509 triangle[1] = (i + 1) % numPoints + 1;
510 triangle[2] = i + numPoints + 2;
511 mesh->insertNextTriangle(triangle);
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);
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;
529 Vec3D<T> rotAxis = {-cylinderAxis[1], cylinderAxis[0], 0.0};
531 T rotationAngle = std::acos(cylinderAxis[2]);
538 Vec3D<T> translationVector;
539 for (
int i = 0; i < 3; ++i) {
540 translationVector[i] = cylinder->origin[i];
546 FromSurfaceMesh<T, D> mesher(levelSet, mesh);
547 mesher.setRemoveBoundaryTriangles(ignoreBoundaryConditions);
551 void makeCustom(SmartPointer<PointCloud<T, D>> pointCloud) {
553 auto mesh = SmartPointer<Mesh<T>>::New();
554 ConvexHull<T, D>(mesh, pointCloud).apply();
557 FromSurfaceMesh<T, D> mesher(levelSet, mesh);
558 mesher.setRemoveBoundaryTriangles(ignoreBoundaryConditions);