Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2840-g502be23d4d 2025-03-14 06:30:00+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
grid_generator_pipe_junction.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2022 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
18
21
22
24
25namespace
26{
33 namespace PipeSegment
34 {
39 struct AdditionalData
40 {
41 double skeleton_length;
42
43 double cosecant_polar;
44 double cotangent_polar;
45 double cotangent_azimuth_half_right;
46 double cotangent_azimuth_half_left;
47 };
48
49
50
55 inline double
56 compute_z_expansion(const double x,
57 const double y,
58 const AdditionalData &data)
59 {
60 return
61 // Scale the unit cylinder to the correct length.
62 data.skeleton_length
63 // Next, adjust for the polar angle. This part will be zero if all
64 // openings and the bifurcation are located on a plane.
65 + x * data.cotangent_polar
66 // Last, adjust for the azimuth angle.
67 - std::abs(y) * data.cosecant_polar *
68 ((y > 0) ? data.cotangent_azimuth_half_right :
69 data.cotangent_azimuth_half_left);
70 }
71
72
73
80 template <int dim, int spacedim = dim>
81 class Manifold : public ChartManifold<dim, spacedim, 3>
82 {
83 public:
103 Manifold(const Tensor<1, spacedim> &normal_direction,
104 const Tensor<1, spacedim> &direction,
105 const Point<spacedim> &point_on_axis,
106 const AdditionalData &data,
107 const double tolerance = 1e-10);
108
112 virtual std::unique_ptr<::Manifold<dim, spacedim>>
113 clone() const override;
114
121 virtual Point<3>
122 pull_back(const Point<spacedim> &space_point) const override;
123
131 virtual Point<spacedim>
132 push_forward(const Point<3> &chart_point) const override;
133
134 private:
138 const Tensor<1, spacedim> normal_direction;
139
143 const Tensor<1, spacedim> direction;
144
148 const Point<spacedim> point_on_axis;
149
153 const AdditionalData data;
154
158 const double tolerance;
159
164 const Tensor<1, spacedim> dxn;
165 };
166
167
168
169 template <int dim, int spacedim>
170 Manifold<dim, spacedim>::Manifold(
171 const Tensor<1, spacedim> &normal_direction,
172 const Tensor<1, spacedim> &direction,
173 const Point<spacedim> &point_on_axis,
174 const AdditionalData &data,
175 const double tolerance)
177 , normal_direction(normal_direction)
178 , direction(direction)
179 , point_on_axis(point_on_axis)
180 , data(data)
181 , tolerance(tolerance)
182 , dxn(cross_product_3d(direction, normal_direction))
183 {
184 Assert(spacedim == 3,
186 "PipeSegment::Manifold can only be used for spacedim==3!"));
187
188 Assert(std::abs(normal_direction.norm() - 1) < tolerance,
189 ExcMessage("Normal direction must be unit vector."));
190 Assert(std::abs(direction.norm() - 1) < tolerance,
191 ExcMessage("Direction must be unit vector."));
192 Assert(normal_direction * direction < tolerance,
194 "Direction and normal direction must be perpendicular."));
195 }
196
197
198
199 template <int dim, int spacedim>
200 std::unique_ptr<::Manifold<dim, spacedim>>
202 {
203 return std::make_unique<Manifold<dim, spacedim>>(*this);
204 }
205
206
207
208 template <int dim, int spacedim>
210 Manifold<dim, spacedim>::pull_back(const Point<spacedim> &space_point) const
211 {
212 // First find the projection of the given point to the axis.
213 const Tensor<1, spacedim> normalized_point = space_point - point_on_axis;
214 double lambda = normalized_point * direction;
215 const Point<spacedim> projection = point_on_axis + direction * lambda;
216 const Tensor<1, spacedim> p_diff = space_point - projection;
217 const double r = p_diff.norm();
218
219 Assert(r > tolerance * data.skeleton_length,
221 "This class won't handle points on the direction axis."));
222
223 // Then compute the angle between the projection direction and
224 // another vector orthogonal to the direction vector.
225 const double phi =
227 p_diff,
228 /*axis=*/direction);
229
230 // Map the axial coordinate to a cylinder of height one.
231 lambda /= compute_z_expansion(r * std::cos(phi), r * std::sin(phi), data);
232
233 // Return distance from the axis, angle and signed distance on the axis.
234 return {r, phi, lambda};
235 }
236
237
238
239 template <int dim, int spacedim>
241 Manifold<dim, spacedim>::push_forward(const Point<3> &chart_point) const
242 {
243 // Rotate the orthogonal direction by the given angle.
244 const double sine_r = chart_point[0] * std::sin(chart_point[1]);
245 const double cosine_r = chart_point[0] * std::cos(chart_point[1]);
246
247 const Tensor<1, spacedim> intermediate =
248 normal_direction * cosine_r + dxn * sine_r;
249
250 // Map the axial coordinate back to the pipe segment.
251 const double lambda =
252 chart_point[2] * compute_z_expansion(cosine_r, sine_r, data);
253
254 // Finally, put everything together.
255 return point_on_axis + direction * lambda + intermediate;
256 }
257 } // namespace PipeSegment
258} // namespace
259
260
261
262namespace GridGenerator
263{
264 template <int dim, int spacedim>
265 void
267 const std::vector<std::pair<Point<spacedim>, double>> &,
268 const std::pair<Point<spacedim>, double> &,
269 const double)
270 {
272 }
273
274
275
276 // hide the template specialization from doxygen
277#ifndef DOXYGEN
278
279 template <>
280 void
282 const std::vector<std::pair<Point<3>, double>> &openings,
283 const std::pair<Point<3>, double> &bifurcation,
284 const double aspect_ratio)
285 {
286 constexpr unsigned int dim = 3;
287 constexpr unsigned int spacedim = 3;
288 using vector3d = Tensor<1, spacedim, double>;
289
290 constexpr unsigned int n_pipes = 3;
291 constexpr double tolerance = 1.e-12;
292
293 if constexpr (running_in_debug_mode())
294 {
295 // Verify user input.
296 Assert(bifurcation.second > 0,
297 ExcMessage("Invalid input: negative radius."));
298 Assert(openings.size() == n_pipes,
299 ExcMessage("Invalid input: only 3 openings allowed."));
300 for (const auto &opening : openings)
301 Assert(opening.second > 0,
302 ExcMessage("Invalid input: negative radius."));
303 }
304
305 // Each pipe segment will be identified by the index of its opening in the
306 // parameter array. To determine the next and previous entry in the array
307 // for a given index, we create auxiliary functions.
308 const auto cyclic = [](const unsigned int i) -> unsigned int {
309 constexpr unsigned int n_pipes = 3;
310 return (i < (n_pipes - 1)) ? i + 1 : 0;
311 };
312 const auto anticyclic = [](const unsigned int i) -> unsigned int {
313 constexpr unsigned int n_pipes = 3;
314 return (i > 0) ? i - 1 : n_pipes - 1;
315 };
316
317 // Cartesian base represented by unit vectors.
318 const std::array<vector3d, spacedim> directions = {
319 {vector3d({1., 0., 0.}), vector3d({0., 1., 0.}), vector3d({0., 0., 1.})}};
320
321 // The skeleton corresponds to the axis of symmetry in the center of each
322 // pipe segment. Each skeleton vector points from the associated opening to
323 // the common bifurcation point. For convenience, we also compute length and
324 // unit vector of every skeleton vector here.
325 std::array<vector3d, n_pipes> skeleton;
326 for (unsigned int p = 0; p < n_pipes; ++p)
327 skeleton[p] = bifurcation.first - openings[p].first;
328
329 std::array<double, n_pipes> skeleton_length;
330 for (unsigned int p = 0; p < n_pipes; ++p)
331 skeleton_length[p] = skeleton[p].norm();
332
333 // In many assertions that come up below, we will verify the integrity of
334 // the geometry. For this, we introduce a tolerance length which vectors
335 // must exceed to avoid being considered "too short". We relate this length
336 // to the longest pipe segment.
337 const double tolerance_length =
338 tolerance *
339 *std::max_element(skeleton_length.begin(), skeleton_length.end());
340
341 std::array<vector3d, n_pipes> skeleton_unit;
342 for (unsigned int p = 0; p < n_pipes; ++p)
343 {
344 Assert(skeleton_length[p] > tolerance_length,
345 ExcMessage("Invalid input: bifurcation matches opening."));
346 skeleton_unit[p] = skeleton[p] / skeleton_length[p];
347 }
348
349 // To determine the orientation of the pipe segments to each other, we will
350 // construct a plane: starting from the bifurcation point, we will move by
351 // the magnitude one in each of the skeleton directions and span a plane
352 // with the three points we reached.
353 //
354 // The normal vector of this particular plane then describes the edge at
355 // which all pipe segments meet. If we would interpret the bifurcation as a
356 // ball joint, the normal vector would correspond to the polar axis of the
357 // ball.
358 vector3d normal = cross_product_3d(skeleton_unit[1] - skeleton_unit[0],
359 skeleton_unit[2] - skeleton_unit[0]);
360 Assert(normal.norm() > tolerance_length,
361 ExcMessage("Invalid input: all three openings "
362 "are located on one line."));
363 normal /= normal.norm();
364
365 // Projections of all skeleton vectors perpendicular to the normal vector,
366 // or in other words, onto the plane described above.
367 std::array<vector3d, n_pipes> skeleton_plane;
368 for (unsigned int p = 0; p < n_pipes; ++p)
369 {
370 skeleton_plane[p] = skeleton[p] - (skeleton[p] * normal) * normal;
371 Assert(std::abs(skeleton_plane[p] * normal) <
372 tolerance * skeleton_plane[p].norm(),
374 Assert(skeleton_plane[p].norm() > tolerance_length,
375 ExcMessage("Invalid input."));
376 }
377
378 // Create a hyperball domain in 2d that will act as the reference cross
379 // section for each pipe segment.
380 Triangulation<dim - 1, spacedim - 1> tria_base;
382 /*center=*/Point<spacedim - 1>(),
383 /*radius=*/1.);
384
385 // Now move on to actually build the pipe junction geometry!
386 //
387 // For each pipe segment, we create a separate triangulation object which
388 // will be merged with the parameter triangulation in the end.
389 Assert(tria.n_cells() == 0,
390 ExcMessage("The output triangulation object needs to be empty."));
391
392 std::vector<PipeSegment::Manifold<dim, spacedim>> manifolds;
393 for (unsigned int p = 0; p < n_pipes; ++p)
394 {
396
397 //
398 // Step 1: create unit cylinder
399 //
400 // We create a unit cylinder by extrusion from the base cross section.
401 // The number of layers depends on the ratio of the length of the
402 // skeleton and half the minimal radius in the pipe segment. The latter
403 // corresponds to the length in radial direction of the smallest cell in
404 // the base cross section. Further, the aspect ratio of the extruded
405 // cells can be set individually with a function parameter.
406 const unsigned int n_slices =
407 1 + static_cast<unsigned int>(std::ceil(
408 aspect_ratio * skeleton_length[p] /
409 (0.5 * std::min(openings[p].second, bifurcation.second))));
411 n_slices,
412 /*height*/ 1.,
413 pipe);
414
415 // Set all material and manifold indicators on the unit cylinder, simply
416 // because they are easier to handle in this geometry. We will set
417 // boundary indicators at the end of the function. See general
418 // documentation of this function.
419 for (const auto &cell : pipe.active_cell_iterators())
420 {
421 cell->set_material_id(p);
422
423 for (const auto &face : cell->face_iterators())
424 if (face->at_boundary())
425 {
426 const auto center_z = face->center()[2];
427
428 if (std::abs(center_z) < tolerance)
429 {
430 // opening cross section
431 }
432 else if (std::abs(center_z - 1.) < tolerance)
433 {
434 // bifurcation cross section
435 }
436 else
437 {
438 // cone mantle
439 cell->set_all_manifold_ids(n_pipes);
440 face->set_all_manifold_ids(p);
441 }
442 }
443 }
444
445 //
446 // Step 2: transform unit cylinder to pipe segment
447 //
448 // For the given cylinder, we will interpret the base in the xy-plane as
449 // the cross section of the opening, and the top at z=1 as the surface
450 // where all pipe segments meet. On the latter surface, we assign the
451 // section in positive y-direction to face the next (right/cyclic) pipe
452 // segment, and allocate the domain in negative y-direction to border
453 // the previous (left/anticyclic) pipe segment.
454 //
455 // In the end, the transformed pipe segment will look like this:
456 // z z
457 // ^ ^
458 // left | right | /|
459 // anticyclic | cyclic |/ |
460 // /|\ /| |
461 // / | \ / | |
462 // | | | | | |
463 // | | | | | |
464 // ------+----->y ------+----->x
465
466 // Before transforming the unit cylinder however, we compute angle
467 // relations between the skeleton vectors viewed from the bifurcation
468 // point. For this purpose, we interpret the bifurcation as a ball joint
469 // as described above.
470 //
471 // In spherical coordinates, the polar angle describes the kink of the
472 // skeleton vector with respect to the polar axis. If all openings and
473 // the bifurcation are located on a plane, then this angle is pi/2 for
474 // every pipe segment.
475 const double polar_angle =
476 Physics::VectorRelations::angle(skeleton[p], normal);
477 Assert(std::abs(polar_angle) > tolerance &&
478 std::abs(polar_angle - numbers::PI) > tolerance,
479 ExcMessage("Invalid input."));
480
481 // Further, we compute the angles between this pipe segment to the other
482 // two. The angle corresponds to the azimuthal direction if we stick to
483 // the picture of the ball joint.
484 const double azimuth_angle_right =
486 skeleton_plane[cyclic(p)],
487 /*axis=*/normal);
488 Assert(std::abs(azimuth_angle_right) > tolerance,
489 ExcMessage("Invalid input: at least two openings located "
490 "in same direction from bifurcation"));
491
492 const double azimuth_angle_left =
494 skeleton_plane[anticyclic(p)],
495 /*axis=*/-normal);
496 Assert(std::abs(azimuth_angle_left) > tolerance,
497 ExcMessage("Invalid input: at least two openings located "
498 "in same direction from bifurcation"));
499
500 // We compute some trigonometric relations with these angles, and store
501 // them conveniently in a struct to be reused later.
502 PipeSegment::AdditionalData data;
503 data.skeleton_length = skeleton_length[p];
504 data.cosecant_polar = 1. / std::sin(polar_angle);
505 data.cotangent_polar = std::cos(polar_angle) * data.cosecant_polar;
506 data.cotangent_azimuth_half_right = std::cos(.5 * azimuth_angle_right) /
507 std::sin(.5 * azimuth_angle_right);
508 data.cotangent_azimuth_half_left =
509 std::cos(.5 * azimuth_angle_left) / std::sin(.5 * azimuth_angle_left);
510
511 // Now transform the cylinder as described above.
512 const auto pipe_segment_transform =
513 [&](const Point<spacedim> &pt) -> Point<spacedim> {
514 // We transform the cylinder in x- and y-direction to become a
515 // truncated cone, similarly to GridGenerator::truncated_cone().
516 const double r_factor =
517 (bifurcation.second - openings[p].second) * pt[2] +
518 openings[p].second;
519 const double x_new = r_factor * pt[0];
520 const double y_new = r_factor * pt[1];
521
522 // Further, to be able to smoothly merge all pipe segments at the
523 // bifurcation, we also need to transform in z-direction.
524 const double z_factor =
525 PipeSegment::compute_z_expansion(x_new, y_new, data);
526 Assert(z_factor > 0,
527 ExcMessage("Invalid input: at least one pipe segment "
528 "is not long enough in this configuration"));
529 const double z_new = z_factor * pt[2];
530
531 return {x_new, y_new, z_new};
532 };
533 GridTools::transform(pipe_segment_transform, pipe);
534
535 //
536 // Step 3: rotate pipe segment to match skeleton direction
537 //
538 // The symmetry axis of the pipe segment in its current state points in
539 // positive z-direction. We rotate the pipe segment that its symmetry
540 // axis matches the direction of the skeleton vector. For this purpose,
541 // we rotate the pipe segment around the axis that is described by the
542 // cross product of both vectors.
543 const double rotation_angle =
544 Physics::VectorRelations::angle(directions[2], skeleton_unit[p]);
545 const vector3d rotation_axis = [&]() {
546 const vector3d rotation_axis =
547 cross_product_3d(directions[2], skeleton_unit[p]);
548 const double norm = rotation_axis.norm();
549 if (norm < tolerance)
550 return directions[1];
551 else
552 return rotation_axis / norm;
553 }();
554 const Tensor<2, spacedim, double> rotation_matrix =
556 rotation_axis, rotation_angle);
558 [&](const Point<spacedim> &pt) { return rotation_matrix * pt; },
559 pipe);
560
561 //
562 // Step 4: rotate laterally to align pipe segments
563 //
564 // On the unit cylinder, we find that the edge on which all pipe
565 // segments meet is parallel to the x-axis. After the transformation to
566 // the pipe segment, we notice that this statement still holds for the
567 // projection of this edge onto the xy-plane, which corresponds to the
568 // cross section of the opening.
569 //
570 // With the latest rotation however, this is no longer the case. We
571 // rotate the unit vector in x-direction in the same fashion, which
572 // gives us the current direction of the projected edge.
573 const vector3d Rx = rotation_matrix * directions[0];
574
575 // To determine how far we need to rotate, we also need to project the
576 // polar axis of the bifurcation ball joint into the same plane.
577 const vector3d normal_projected_on_opening =
578 normal - (normal * skeleton_unit[p]) * skeleton_unit[p];
579
580 // Both the projected normal and Rx must be in the opening plane.
581 Assert(std::abs(skeleton_unit[p] * normal_projected_on_opening) <
582 tolerance,
584 Assert(std::abs(skeleton_unit[p] * Rx) < tolerance, ExcInternalError());
585
586 // Now we laterally rotate the pipe segment around its own symmetry axis
587 // that the edge matches the polar axis.
588 const double lateral_angle =
590 normal_projected_on_opening,
591 /*axis=*/skeleton_unit[p]);
592 GridTools::rotate(skeleton_unit[p], lateral_angle, pipe);
593
594 //
595 // Step 5: shift to final position and merge this pipe into the entire
596 // assembly
597 //
598 GridTools::shift(openings[p].first, pipe);
599
600 // Create a manifold object for the mantle of this particular pipe
601 // segment. Since GridGenerator::merge_triangulations() does not copy
602 // manifold objects, but just IDs if requested, we will copy them to
603 // the final triangulation later.
604 manifolds.emplace_back(
605 /*normal_direction=*/normal_projected_on_opening /
606 normal_projected_on_opening.norm(),
607 /*direction=*/skeleton_unit[p],
608 /*point_on_axis=*/openings[p].first,
609 data,
610 tolerance);
611
613 pipe, tria, tria, tolerance_length, /*copy_manifold_ids=*/true);
614 }
615
616 for (unsigned int p = 0; p < n_pipes; ++p)
617 tria.set_manifold(p, manifolds[p]);
618 tria.set_manifold(n_pipes, FlatManifold<3>());
619
620 // Since GridGenerator::merge_triangulations() does not copy boundary IDs
621 // either, we need to set them after the final geometry is created. Luckily,
622 // boundary IDs match with material IDs, so we simply translate them with
623 // the help of manifold IDs to identify openings.
624 for (const auto &cell : tria.active_cell_iterators())
625 for (const auto &face : cell->face_iterators())
626 if (face->at_boundary())
627 {
628 if (face->manifold_id() == numbers::flat_manifold_id ||
629 face->manifold_id() == n_pipes)
630 // opening cross section
631 face->set_boundary_id(cell->material_id());
632 else
633 // cone mantle
634 face->set_boundary_id(n_pipes);
635 }
636 }
637
638#endif
639
640} // namespace GridGenerator
641
642
643// explicit instantiations
644#include "grid/grid_generator_pipe_junction.inst"
645
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const =0
Definition point.h:113
numbers::NumberTraits< Number >::real_type norm() const
unsigned int n_cells() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_NOT_IMPLEMENTED()
Point< 2 > second
Definition grid_out.cc:4633
Point< 2 > first
Definition grid_out.cc:4632
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
std::vector< index_type > data
Definition mpi.cc:746
void hyper_ball_balanced(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1.)
void pipe_junction(Triangulation< dim, spacedim > &tria, const std::vector< std::pair< Point< spacedim >, double > > &openings, const std::pair< Point< spacedim >, double > &bifurcation, const double aspect_ratio=0.5)
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
void rotate(const double angle, Triangulation< dim, spacedim > &triangulation)
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
Tensor< 2, 3, Number > rotation_matrix_3d(const Tensor< 1, 3, Number > &axis, const Number &angle)
Number signed_angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b, const Tensor< 1, spacedim, Number > &axis)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
constexpr double PI
Definition numbers.h:239
constexpr types::manifold_id flat_manifold_id
Definition types.h:336
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)