Reference documentation for deal.II version GIT c14369f203 2023-10-01 07:40: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\}}\)
grid_generator_pipe_junction.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 - 2023 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
9 // Public License as published by the Free Software Foundation; either
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 
18 #include <deal.II/grid/manifold.h>
19 
22 
23 
25 
26 namespace
27 {
34  namespace PipeSegment
35  {
40  struct AdditionalData
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,
59  const AdditionalData &data)
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,
107  const AdditionalData &data,
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 
154  const AdditionalData data;
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,
175  const AdditionalData &data,
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,
186  ExcMessage(
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,
194  ExcMessage(
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>
210  Point<3>
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,
221  ExcMessage(
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 
263 namespace 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,
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, 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(),
372  ExcInternalError());
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>(),
382  /*radius=*/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 =
484  Physics::VectorRelations::signed_angle(skeleton_plane[p],
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 =
492  Physics::VectorRelations::signed_angle(skeleton_plane[p],
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.
501  PipeSegment::AdditionalData data;
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  }();
553  const Tensor<2, spacedim, double> rotation_matrix =
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,
582  ExcInternalError());
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
numbers::NumberTraits< Number >::real_type norm() const
unsigned int n_cells() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 2 > second
Definition: grid_out.cc:4615
Point< 2 > first
Definition: grid_out.cc:4614
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
Expression ceil(const Expression &x)
void hyper_ball_balanced(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1.)
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 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 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)
Definition: grid_tools.cc:2142
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:2132
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
Point< dim > push_forward(const TopoDS_Shape &in_shape, const double u, const double v)
Definition: utilities.cc:835
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 1, dim, Number > pull_back(const Tensor< 1, dim, Number > &v, const Tensor< 2, dim, Number > &F)
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:260
const types::manifold_id flat_manifold_id
Definition: types.h:286
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
const ::Triangulation< dim, spacedim > & tria