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