Reference documentation for deal.II version GIT relicensing-720-g4ca08d40c8 2024-05-22 18:50:02+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\}}\)
Loading...
Searching...
No Matches
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# ifdef DEBUG
294 // Verify user input.
295 Assert(bifurcation.second > 0,
296 ExcMessage("Invalid input: negative radius."));
297 Assert(openings.size() == n_pipes,
298 ExcMessage("Invalid input: only 3 openings allowed."));
299 for (const auto &opening : openings)
300 Assert(opening.second > 0, ExcMessage("Invalid input: negative radius."));
301# endif
302
303 // Each pipe segment will be identified by the index of its opening in the
304 // parameter array. To determine the next and previous entry in the array
305 // for a given index, we create auxiliary functions.
306 const auto cyclic = [](const unsigned int i) -> unsigned int {
307 constexpr unsigned int n_pipes = 3;
308 return (i < (n_pipes - 1)) ? i + 1 : 0;
309 };
310 const auto anticyclic = [](const unsigned int i) -> unsigned int {
311 constexpr unsigned int n_pipes = 3;
312 return (i > 0) ? i - 1 : n_pipes - 1;
313 };
314
315 // Cartesian base represented by unit vectors.
316 const std::array<vector3d, spacedim> directions = {
317 {vector3d({1., 0., 0.}), vector3d({0., 1., 0.}), vector3d({0., 0., 1.})}};
318
319 // The skeleton corresponds to the axis of symmetry in the center of each
320 // pipe segment. Each skeleton vector points from the associated opening to
321 // the common bifurcation point. For convenience, we also compute length and
322 // unit vector of every skeleton vector here.
323 std::array<vector3d, n_pipes> skeleton;
324 for (unsigned int p = 0; p < n_pipes; ++p)
325 skeleton[p] = bifurcation.first - openings[p].first;
326
327 std::array<double, n_pipes> skeleton_length;
328 for (unsigned int p = 0; p < n_pipes; ++p)
329 skeleton_length[p] = skeleton[p].norm();
330
331 // In many assertions that come up below, we will verify the integrity of
332 // the geometry. For this, we introduce a tolerance length which vectors
333 // must exceed to avoid being considered "too short". We relate this length
334 // to the longest pipe segment.
335 const double tolerance_length =
336 tolerance *
337 *std::max_element(skeleton_length.begin(), skeleton_length.end());
338
339 std::array<vector3d, n_pipes> skeleton_unit;
340 for (unsigned int p = 0; p < n_pipes; ++p)
341 {
342 Assert(skeleton_length[p] > tolerance_length,
343 ExcMessage("Invalid input: bifurcation matches opening."));
344 skeleton_unit[p] = skeleton[p] / skeleton_length[p];
345 }
346
347 // To determine the orientation of the pipe segments to each other, we will
348 // construct a plane: starting from the bifurcation point, we will move by
349 // the magnitude one in each of the skeleton directions and span a plane
350 // with the three points we reached.
351 //
352 // The normal vector of this particular plane then describes the edge at
353 // which all pipe segments meet. If we would interpret the bifurcation as a
354 // ball joint, the normal vector would correspond to the polar axis of the
355 // ball.
356 vector3d normal = cross_product_3d(skeleton_unit[1] - skeleton_unit[0],
357 skeleton_unit[2] - skeleton_unit[0]);
358 Assert(normal.norm() > tolerance_length,
359 ExcMessage("Invalid input: all three openings "
360 "are located on one line."));
361 normal /= normal.norm();
362
363 // Projections of all skeleton vectors perpendicular to the normal vector,
364 // or in other words, onto the plane described above.
365 std::array<vector3d, n_pipes> skeleton_plane;
366 for (unsigned int p = 0; p < n_pipes; ++p)
367 {
368 skeleton_plane[p] = skeleton[p] - (skeleton[p] * normal) * normal;
369 Assert(std::abs(skeleton_plane[p] * normal) <
370 tolerance * skeleton_plane[p].norm(),
372 Assert(skeleton_plane[p].norm() > tolerance_length,
373 ExcMessage("Invalid input."));
374 }
375
376 // Create a hyperball domain in 2d that will act as the reference cross
377 // section for each pipe segment.
378 Triangulation<dim - 1, spacedim - 1> tria_base;
380 /*center=*/Point<spacedim - 1>(),
381 /*radius=*/1.);
382
383 // Now move on to actually build the pipe junction geometry!
384 //
385 // For each pipe segment, we create a separate triangulation object which
386 // will be merged with the parameter triangulation in the end.
387 Assert(tria.n_cells() == 0,
388 ExcMessage("The output triangulation object needs to be empty."));
389
390 std::vector<PipeSegment::Manifold<dim, spacedim>> manifolds;
391 for (unsigned int p = 0; p < n_pipes; ++p)
392 {
394
395 //
396 // Step 1: create unit cylinder
397 //
398 // We create a unit cylinder by extrusion from the base cross section.
399 // The number of layers depends on the ratio of the length of the
400 // skeleton and half the minimal radius in the pipe segment. The latter
401 // corresponds to the length in radial direction of the smallest cell in
402 // the base cross section. Further, the aspect ratio of the extruded
403 // cells can be set individually with a function parameter.
404 const unsigned int n_slices =
405 1 + static_cast<unsigned int>(std::ceil(
406 aspect_ratio * skeleton_length[p] /
407 (0.5 * std::min(openings[p].second, bifurcation.second))));
409 n_slices,
410 /*height*/ 1.,
411 pipe);
412
413 // Set all material and manifold indicators on the unit cylinder, simply
414 // because they are easier to handle in this geometry. We will set
415 // boundary indicators at the end of the function. See general
416 // documentation of this function.
417 for (const auto &cell : pipe.active_cell_iterators())
418 {
419 cell->set_material_id(p);
420
421 for (const auto &face : cell->face_iterators())
422 if (face->at_boundary())
423 {
424 const auto center_z = face->center()[2];
425
426 if (std::abs(center_z) < tolerance)
427 {
428 // opening cross section
429 }
430 else if (std::abs(center_z - 1.) < tolerance)
431 {
432 // bifurcation cross section
433 }
434 else
435 {
436 // cone mantle
437 cell->set_all_manifold_ids(n_pipes);
438 face->set_all_manifold_ids(p);
439 }
440 }
441 }
442
443 //
444 // Step 2: transform unit cylinder to pipe segment
445 //
446 // For the given cylinder, we will interpret the base in the xy-plane as
447 // the cross section of the opening, and the top at z=1 as the surface
448 // where all pipe segments meet. On the latter surface, we assign the
449 // section in positive y-direction to face the next (right/cyclic) pipe
450 // segment, and allocate the domain in negative y-direction to border
451 // the previous (left/anticyclic) pipe segment.
452 //
453 // In the end, the transformed pipe segment will look like this:
454 // z z
455 // ^ ^
456 // left | right | /|
457 // anticyclic | cyclic |/ |
458 // /|\ /| |
459 // / | \ / | |
460 // | | | | | |
461 // | | | | | |
462 // ------+----->y ------+----->x
463
464 // Before transforming the unit cylinder however, we compute angle
465 // relations between the skeleton vectors viewed from the bifurcation
466 // point. For this purpose, we interpret the bifurcation as a ball joint
467 // as described above.
468 //
469 // In spherical coordinates, the polar angle describes the kink of the
470 // skeleton vector with respect to the polar axis. If all openings and
471 // the bifurcation are located on a plane, then this angle is pi/2 for
472 // every pipe segment.
473 const double polar_angle =
474 Physics::VectorRelations::angle(skeleton[p], normal);
475 Assert(std::abs(polar_angle) > tolerance &&
476 std::abs(polar_angle - numbers::PI) > tolerance,
477 ExcMessage("Invalid input."));
478
479 // Further, we compute the angles between this pipe segment to the other
480 // two. The angle corresponds to the azimuthal direction if we stick to
481 // the picture of the ball joint.
482 const double azimuth_angle_right =
484 skeleton_plane[cyclic(p)],
485 /*axis=*/normal);
486 Assert(std::abs(azimuth_angle_right) > tolerance,
487 ExcMessage("Invalid input: at least two openings located "
488 "in same direction from bifurcation"));
489
490 const double azimuth_angle_left =
492 skeleton_plane[anticyclic(p)],
493 /*axis=*/-normal);
494 Assert(std::abs(azimuth_angle_left) > tolerance,
495 ExcMessage("Invalid input: at least two openings located "
496 "in same direction from bifurcation"));
497
498 // We compute some trigonometric relations with these angles, and store
499 // them conveniently in a struct to be reused later.
500 PipeSegment::AdditionalData data;
501 data.skeleton_length = skeleton_length[p];
502 data.cosecant_polar = 1. / std::sin(polar_angle);
503 data.cotangent_polar = std::cos(polar_angle) * data.cosecant_polar;
504 data.cotangent_azimuth_half_right = std::cos(.5 * azimuth_angle_right) /
505 std::sin(.5 * azimuth_angle_right);
506 data.cotangent_azimuth_half_left =
507 std::cos(.5 * azimuth_angle_left) / std::sin(.5 * azimuth_angle_left);
508
509 // Now transform the cylinder as described above.
510 const auto pipe_segment_transform =
511 [&](const Point<spacedim> &pt) -> Point<spacedim> {
512 // We transform the cylinder in x- and y-direction to become a
513 // truncated cone, similarly to GridGenerator::truncated_cone().
514 const double r_factor =
515 (bifurcation.second - openings[p].second) * pt[2] +
516 openings[p].second;
517 const double x_new = r_factor * pt[0];
518 const double y_new = r_factor * pt[1];
519
520 // Further, to be able to smoothly merge all pipe segments at the
521 // bifurcation, we also need to transform in z-direction.
522 const double z_factor =
523 PipeSegment::compute_z_expansion(x_new, y_new, data);
524 Assert(z_factor > 0,
525 ExcMessage("Invalid input: at least one pipe segment "
526 "is not long enough in this configuration"));
527 const double z_new = z_factor * pt[2];
528
529 return {x_new, y_new, z_new};
530 };
531 GridTools::transform(pipe_segment_transform, pipe);
532
533 //
534 // Step 3: rotate pipe segment to match skeleton direction
535 //
536 // The symmetry axis of the pipe segment in its current state points in
537 // positive z-direction. We rotate the pipe segment that its symmetry
538 // axis matches the direction of the skeleton vector. For this purpose,
539 // we rotate the pipe segment around the axis that is described by the
540 // cross product of both vectors.
541 const double rotation_angle =
542 Physics::VectorRelations::angle(directions[2], skeleton_unit[p]);
543 const vector3d rotation_axis = [&]() {
544 const vector3d rotation_axis =
545 cross_product_3d(directions[2], skeleton_unit[p]);
546 const double norm = rotation_axis.norm();
547 if (norm < tolerance)
548 return directions[1];
549 else
550 return rotation_axis / norm;
551 }();
552 const Tensor<2, spacedim, double> rotation_matrix =
554 rotation_axis, rotation_angle);
555 GridTools::transform(
556 [&](const Point<spacedim> &pt) { return rotation_matrix * pt; },
557 pipe);
558
559 //
560 // Step 4: rotate laterally to align pipe segments
561 //
562 // On the unit cylinder, we find that the edge on which all pipe
563 // segments meet is parallel to the x-axis. After the transformation to
564 // the pipe segment, we notice that this statement still holds for the
565 // projection of this edge onto the xy-plane, which corresponds to the
566 // cross section of the opening.
567 //
568 // With the latest rotation however, this is no longer the case. We
569 // rotate the unit vector in x-direction in the same fashion, which
570 // gives us the current direction of the projected edge.
571 const vector3d Rx = rotation_matrix * directions[0];
572
573 // To determine how far we need to rotate, we also need to project the
574 // polar axis of the bifurcation ball joint into the same plane.
575 const vector3d normal_projected_on_opening =
576 normal - (normal * skeleton_unit[p]) * skeleton_unit[p];
577
578 // Both the projected normal and Rx must be in the opening plane.
579 Assert(std::abs(skeleton_unit[p] * normal_projected_on_opening) <
580 tolerance,
582 Assert(std::abs(skeleton_unit[p] * Rx) < tolerance, ExcInternalError());
583
584 // Now we laterally rotate the pipe segment around its own symmetry axis
585 // that the edge matches the polar axis.
586 const double lateral_angle =
588 normal_projected_on_opening,
589 /*axis=*/skeleton_unit[p]);
590 GridTools::rotate(skeleton_unit[p], lateral_angle, pipe);
591
592 //
593 // Step 5: shift to final position and merge this pipe into the entire
594 // assembly
595 //
596 GridTools::shift(openings[p].first, pipe);
597
598 // Create a manifold object for the mantle of this particular pipe
599 // segment. Since GridGenerator::merge_triangulations() does not copy
600 // manifold objects, but just IDs if requested, we will copy them to
601 // the final triangulation later.
602 manifolds.emplace_back(
603 /*normal_direction=*/normal_projected_on_opening /
604 normal_projected_on_opening.norm(),
605 /*direction=*/skeleton_unit[p],
606 /*point_on_axis=*/openings[p].first,
607 data,
608 tolerance);
609
611 pipe, tria, tria, tolerance_length, /*copy_manifold_ids=*/true);
612 }
613
614 for (unsigned int p = 0; p < n_pipes; ++p)
615 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
numbers::NumberTraits< Number >::real_type norm() const
unsigned int n_cells() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
#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)
#define DEAL_II_NOT_IMPLEMENTED()
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)
spacedim const Point< spacedim > & p
Definition grid_tools.h:990
void rotate(const double angle, Triangulation< dim, spacedim > &triangulation)
const Triangulation< dim, spacedim > & tria
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
if(marked_vertices.size() !=0) for(auto it
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
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)
static constexpr double PI
Definition numbers.h:259
const types::manifold_id flat_manifold_id
Definition types.h:325
::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 > &)