Reference documentation for deal.II version Git 2276b94619 2020-02-28 10:11:55 +0100
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_generator.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 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 
16 #include <deal.II/distributed/fully_distributed_tria.h>
17 #include <deal.II/distributed/shared_tria.h>
18 #include <deal.II/distributed/tria.h>
19 
20 #include <deal.II/grid/grid_generator.h>
21 #include <deal.II/grid/grid_reordering.h>
22 #include <deal.II/grid/grid_tools.h>
23 #include <deal.II/grid/intergrid_map.h>
24 #include <deal.II/grid/manifold_lib.h>
25 #include <deal.II/grid/tria.h>
26 #include <deal.II/grid/tria_accessor.h>
27 #include <deal.II/grid/tria_iterator.h>
28 
29 #include <cmath>
30 #include <limits>
31 
32 
33 DEAL_II_NAMESPACE_OPEN
34 
35 
36 namespace GridGenerator
37 {
38  namespace Airfoil
39  {
41  // airfoil configuration
42  : airfoil_type("NACA")
43  , naca_id("2412")
44  , joukowski_center(-0.1, 0.14)
45  , airfoil_length(1.0)
46  // far field
47  , height(30.0)
48  , length_b2(15.0)
49  // mesh
50  , incline_factor(0.35)
51  , bias_factor(2.5)
52  , refinements(2)
53  , n_subdivision_x_0(3)
54  , n_subdivision_x_1(2)
55  , n_subdivision_x_2(5)
56  , n_subdivision_y(3)
57  , airfoil_sampling_factor(2)
58  {
59  Assert(
61  ExcMessage(
62  "Mesh is to small to enclose airfoil! Choose larger field or smaller"
63  " chord length!"));
64  Assert(incline_factor < 1.0 && incline_factor >= 0.0,
65  ExcMessage("incline_factor has to be in [0,1)!"));
66  }
67 
68 
69 
70  void
72  {
73  prm.enter_subsection("FarField");
74  {
75  prm.add_parameter(
76  "Height",
77  height,
78  "Mesh height measured from airfoil nose to horizontal boundaries");
79  prm.add_parameter(
80  "LengthB2",
81  length_b2,
82  "Length measured from airfoil leading edge to vertical outlet boundary");
83  prm.add_parameter(
84  "InclineFactor",
86  "Define obliqueness of the vertical mesh around the airfoil");
87  }
88  prm.leave_subsection();
89 
90  prm.enter_subsection("AirfoilType");
91  {
92  prm.add_parameter(
93  "Type",
95  "Type of airfoil geometry, either NACA or Joukowski airfoil",
96  Patterns::Selection("NACA|Joukowski"));
97  }
98  prm.leave_subsection();
99 
100  prm.enter_subsection("NACA");
101  {
102  prm.add_parameter("NacaId", naca_id, "Naca serial number");
103  }
104  prm.leave_subsection();
105 
106  prm.enter_subsection("Joukowski");
107  {
108  prm.add_parameter("Center",
110  "Joukowski circle center coordinates");
111  prm.add_parameter("AirfoilLength",
113  "Joukowski airfoil length leading to trailing edge");
114  }
115  prm.leave_subsection();
116 
117  prm.enter_subsection("Mesh");
118  {
119  prm.add_parameter("Refinements",
120  refinements,
121  "Number of global refinements");
122  prm.add_parameter(
123  "NumberSubdivisionX0",
125  "Number of subdivisions along the airfoil in blocks with material ID 1 and 4");
126  prm.add_parameter(
127  "NumberSubdivisionX1",
129  "Number of subdivisions along the airfoil in blocks with material ID 2 and 5");
130  prm.add_parameter(
131  "NumberSubdivisionX2",
133  "Number of subdivisions in horizontal direction on the right of the trailing edge, i.e., blocks with material ID 3 and 6");
134  prm.add_parameter("NumberSubdivisionY",
136  "Number of subdivisions normal to airfoil");
137  prm.add_parameter(
138  "BiasFactor",
139  bias_factor,
140  "Factor to obtain a finer mesh at the airfoil surface");
141  }
142  prm.leave_subsection();
143  }
144 
145 
146  namespace
147  {
151  class MeshGenerator
152  {
153  public:
154  // IDs of the mesh blocks
155  static const unsigned int id_block_1 = 1;
156  static const unsigned int id_block_2 = 2;
157  static const unsigned int id_block_3 = 3;
158  static const unsigned int id_block_4 = 4;
159  static const unsigned int id_block_5 = 5;
160  static const unsigned int id_block_6 = 6;
161 
165  MeshGenerator(const AdditionalData &data)
166  : refinements(data.refinements)
171  , height(data.height)
172  , length_b2(data.length_b2)
174  , bias_factor(data.bias_factor)
175  , edge_length(1.0)
176  , n_cells_x_0(Utilities::pow(2, refinements) * n_subdivision_x_0)
177  , n_cells_x_1(Utilities::pow(2, refinements) * n_subdivision_x_1)
178  , n_cells_x_2(Utilities::pow(2, refinements) * n_subdivision_x_2)
179  , n_cells_y(Utilities::pow(2, refinements) * n_subdivision_y)
180  , n_points_on_each_side(n_cells_x_0 + n_cells_x_1 + 1)
181  // create points on the airfoil
182  , airfoil_1D(set_airfoil_length(
183  // call either the 'joukowski' or 'naca' static member function
184  data.airfoil_type == "Joukowski" ?
185  joukowski(data.joukowski_center,
186  n_points_on_each_side,
188  (data.airfoil_type == "NACA" ?
189  naca(data.naca_id,
190  n_points_on_each_side,
192  std::array<std::vector<Point<2>>, 2>{
193  {std::vector<Point<2>>{Point<2>(0), Point<2>(1)},
194  std::vector<Point<2>>{
195  Point<2>(0),
196  Point<2>(
197  1)}}} /* dummy vector since we are asserting later*/),
198  data.airfoil_length))
199  , end_b0_x_u(airfoil_1D[0][n_cells_x_0](0))
200  , end_b0_x_l(airfoil_1D[1][n_cells_x_0](0))
201  , nose_x(airfoil_1D[0].front()(0))
202  , tail_x(airfoil_1D[0].back()(0))
203  , tail_y(airfoil_1D[0].back()(1))
204  , center_mesh(0.5 * std::abs(end_b0_x_u + end_b0_x_l))
205  , length_b1_x(tail_x - center_mesh)
206  , gamma(std::atan(height /
207  (edge_length + std::abs(nose_x - center_mesh))))
208  // points on coarse grid
209  // coarse grid has to be symmetric in respect to x-axis to allow
210  // periodic BC and make sure that interpolate() works
211  , A(nose_x - edge_length, 0)
212  , B(nose_x, 0)
213  , C(center_mesh, +std::abs(nose_x - center_mesh) * std::tan(gamma))
214  , D(center_mesh, height)
215  , E(center_mesh, -std::abs(nose_x - center_mesh) * std::tan(gamma))
216  , F(center_mesh, -height)
217  , G(tail_x, height)
218  , H(tail_x, 0)
219  , I(tail_x, -height)
220  , J(tail_x + length_b2, 0)
221  , K(J(0), G(1))
222  , L(J(0), I(1))
223  {
224  Assert(data.airfoil_type == "Joukowski" ||
225  data.airfoil_type == "NACA",
226  ExcMessage("Unknown airfoil type."));
227  }
228 
233  Triangulation<2> & tria_grid,
234  std::vector<GridTools::PeriodicFacePair<
235  typename Triangulation<2>::cell_iterator>> *periodic_faces) const
236  {
237  make_coarse_grid(tria_grid);
238 
239  set_boundary_ids(tria_grid);
240 
241  if (periodic_faces != nullptr)
242  {
244  tria_grid, 5, 4, 1, *periodic_faces);
245  tria_grid.add_periodicity(*periodic_faces);
246  }
247 
248  tria_grid.refine_global(refinements);
249  interpolate(tria_grid);
250  }
251 
257  std::vector<GridTools::PeriodicFacePair<
258  typename Triangulation<2>::cell_iterator>> *periodic_faces) const
259  {
260  (void)parallel_grid;
261  (void)periodic_faces;
262 
263  AssertThrow(false, ExcMessage("Not implemented, yet!")); // TODO [PM]
264  }
265 
266  private:
267  // number of global refinements
268  const unsigned int refinements;
269 
270  // numer of subdivisions of coarse grid in blocks 1 and 4
271  const unsigned int n_subdivision_x_0;
272 
273  // numer of subdivisions of coarse grid in blocks 2 and 5
274  const unsigned int n_subdivision_x_1;
275 
276  // numer of subdivisions of coarse grid in blocks 3 and 6
277  const unsigned int n_subdivision_x_2;
278 
279  // numer of subdivisions of coarse grid in all blocks (normal to
280  // airfoil or in y-direction, respectively)
281  const unsigned int n_subdivision_y;
282 
283  // height of mesh, i.e. length JK or JL and radius of semicircle
284  // (C-Mesh) that arises after interpolation in blocks 1 and 4
285  const double height;
286 
287  // lenght block 3 and 6
288  const double length_b2;
289 
290  // factor to move points G and I horizontal to the right, i.e. make
291  // faces HG and HI inclined instead of vertical
292  const double incline_factor;
293 
294  // bias factor (if factor goes to zero than equal y = x)
295  const double bias_factor;
296 
297  // x-distance between coarse grid vertices A and B, i.e. used only once;
298  const double edge_length;
299 
300  // number of cells (after refining) in block 1 and 4 along airfoil
301  const unsigned int n_cells_x_0;
302 
303  // number of cells (after refining) in block 2 and 5 along airfoil
304  const unsigned int n_cells_x_1;
305 
306  // number of cells (after refining) in block 3 and 6 in x-direction
307  const unsigned int n_cells_x_2;
308 
309  // number of cells (after refining) in all blocks normal to airfoil or
310  // in y-direction, respectively
311  const unsigned int n_cells_y;
312 
313  // number of airfoil points on each side
314  const unsigned int n_points_on_each_side;
315 
316  // vector containing upper/lower airfoil points. First and last point
317  // are identical
318  const std::array<std::vector<Point<2>>, 2> airfoil_1D;
319 
320  // x-coordinate of n-th airfoilpoint where n indicates number of cells
321  // in block 1. end_b0_x_u = end_b0_x_l for symmetric airfoils
322  const double end_b0_x_u;
323 
324  // x-coordinate of n-th airfoilpoint where n indicates number of cells
325  // in block 4. end_b0_x_u = end_b0_x_l for symmetric airfoils
326  const double end_b0_x_l;
327 
328  // x-coordinate of first airfoil point in airfoil_1D[0] and
329  // airfoil_1D[1]
330  const double nose_x;
331 
332  // x-coordinate of last airfoil point in airfoil_1D[0] and airfoil_1D[1]
333  const double tail_x;
334 
335  // y-coordinate of last airfoil point in airfoil_1D[0] and airfoil_1D[1]
336  const double tail_y;
337 
338  // x-coordinate of C,D,E,F indicating ending of blocks 1 and 4 or
339  // beginning of blocks 2 and 5, respectively
340  const double center_mesh;
341 
342  // length of blocks 2 and 5
343  const double length_b1_x;
344 
345  // angle enclosed between faces DAB and FAB
346  const double gamma;
347 
348 
349 
370  const Point<2> A, B, C, D, E, F, G, H, I, J, K, L;
371 
372 
373 
409  static std::array<std::vector<Point<2>>, 2>
410  joukowski(const Point<2> centerpoint,
411  const unsigned int number_points,
412  const unsigned int factor)
413  {
414  std::array<std::vector<Point<2>>, 2> airfoil_1D;
415  const unsigned int total_points = 2 * number_points - 2;
416  const unsigned int n_airfoilpoints = factor * total_points;
417  // joukowski points on the entire airfoil, i.e. upper and lower side
418  const auto jouk_points =
419  joukowski_transform(joukowski_circle(centerpoint, n_airfoilpoints));
420 
421  // vectors to collect airfoil points on either upper or lower side
422  std::vector<Point<2>> upper_points;
423  std::vector<Point<2>> lower_points;
424 
425  {
426  // find point on nose and point on tail
427  unsigned int nose_index = 0;
428  unsigned int tail_index = 0;
429  double nose_x_coordinate = 0;
430  double tail_x_coordinate = 0;
431 
432 
433  // find index in vector to nose point (min) and tail point (max)
434  for (unsigned int i = 0; i < jouk_points.size(); i++)
435  {
436  if (jouk_points[i](0) < nose_x_coordinate)
437  {
438  nose_x_coordinate = jouk_points[i](0);
439  nose_index = i;
440  }
441  if (jouk_points[i](0) > tail_x_coordinate)
442  {
443  tail_x_coordinate = jouk_points[i](0);
444  tail_index = i;
445  }
446  }
447 
448  // copy point on upper side of airfoil
449  for (unsigned int i = tail_index; i < jouk_points.size(); i++)
450  upper_points.emplace_back(jouk_points[i]);
451  for (unsigned int i = 0; i <= nose_index; i++)
452  upper_points.emplace_back(jouk_points[i]);
453  std::reverse(upper_points.begin(), upper_points.end());
454 
455  // copy point on lower side of airfoil
456  lower_points.insert(lower_points.end(),
457  jouk_points.begin() + nose_index,
458  jouk_points.begin() + tail_index + 1);
459  }
460 
461  airfoil_1D[0] = make_points_equidistant(upper_points, number_points);
462  airfoil_1D[1] = make_points_equidistant(lower_points, number_points);
463 
464  // move nose to origin
465  auto move_nose_to_origin = [](std::vector<Point<2>> &vector) {
466  const double nose_x_pos = vector.front()(0);
467  for (auto &i : vector)
468  i(0) -= nose_x_pos;
469  };
470 
471  move_nose_to_origin(airfoil_1D[1]);
472  move_nose_to_origin(airfoil_1D[0]);
473 
474  return airfoil_1D;
475  }
476 
501  static std::vector<Point<2>>
502  joukowski_circle(const Point<2> & center,
503  const unsigned int number_points)
504  {
505  std::vector<Point<2>> circle_points;
506 
507  // Create Circle with number_points - points
508  // unsigned int number_points = 2 * points_per_side - 2;
509 
510  // Calculate radius so that point (x=1|y=0) is enclosed - requirement
511  // for Joukowski transform
512  const double radius = std::sqrt(center(1) * center(1) +
513  (1 - center(0)) * (1 - center(0)));
514  const double radius_test = std::sqrt(
515  center(1) * center(1) + (1 + center(0)) * (1 + center(0)));
516  // Make sure point (x=-1|y=0) is enclosed by the circle
517  (void)radius_test;
518  AssertThrow(
519  radius_test < radius,
520  ExcMessage(
521  "Error creating lower circle: Circle for Joukowski-transform does"
522  " not enclose point zeta = -1! Choose different center "
523  "coordinate."));
524  // Create a full circle with radius 'radius' around Point 'center' of
525  // (number_points) equidistant points.
526  const double theta = 2 * numbers::PI / number_points;
527  // first point is leading edge then counterclockwise
528  for (unsigned int i = 0; i < number_points; i++)
529  circle_points.emplace_back(center[0] - radius * cos(i * theta),
530  center[1] - radius * sin(i * theta));
531 
532  return circle_points;
533  }
534 
542  static std::vector<Point<2>>
543  joukowski_transform(const std::vector<Point<2>> &circle_points)
544  {
545  std::vector<Point<2>> joukowski_points(circle_points.size());
546 
547  // transform each point
548  for (unsigned int i = 0; i < circle_points.size(); i++)
549  {
550  const double chi = circle_points[i](0);
551  const double eta = circle_points[i](1);
552  const std::complex<double> zeta(chi, eta);
553  const std::complex<double> z = zeta + 1. / zeta;
554 
555  joukowski_points[i] = {real(z), imag(z)};
556  }
557  return joukowski_points;
558  }
559 
576  static std::array<std::vector<Point<2>>, 2>
577  naca(const std::string &serialnumber,
578  const unsigned int number_points,
579  const unsigned int factor)
580  {
581  // number of non_equidistant airfoilpoints among which will be
582  // interpolated
583  const unsigned int n_airfoilpoints = factor * number_points;
584 
585  // create equidistant airfoil points for upper and lower side
586  return {{make_points_equidistant(
587  naca_create_points(serialnumber, n_airfoilpoints, true),
588  number_points),
589  make_points_equidistant(
590  naca_create_points(serialnumber, n_airfoilpoints, false),
591  number_points)}};
592  }
593 
605  static std::vector<Point<2>>
606  naca_create_points(const std::string &serialnumber,
607  const unsigned int number_points,
608  const bool is_upper)
609  {
610  Assert(serialnumber.length() == 4,
611  ExcMessage("This NACA-serial number is not implemented!"));
612 
613  return naca_create_points_4_digits(serialnumber,
614  number_points,
615  is_upper);
616  }
617 
632  static std::vector<Point<2>>
633  naca_create_points_4_digits(const std::string &serialnumber,
634  const unsigned int number_points,
635  const bool is_upper)
636  {
637  // conversion string (char * ) to int
638  const unsigned int digit_0 = (serialnumber[0] - '0');
639  const unsigned int digit_1 = (serialnumber[1] - '0');
640  const unsigned int digit_2 = (serialnumber[2] - '0');
641  const unsigned int digit_3 = (serialnumber[3] - '0');
642 
643  const unsigned int digit_23 = 10 * digit_2 + digit_3;
644 
645  // maximum thickness in percentage of the cord
646  const double t = static_cast<double>(digit_23) / 100.0;
647 
648  std::vector<Point<2>> naca_points;
649 
650  if (digit_0 == 0 && digit_1 == 0) // is symmetric
651  for (unsigned int i = 0; i < number_points; i++)
652  {
653  const double x = i * 1 / (1.0 * number_points - 1);
654  const double y_t =
655  5 * t *
656  (0.2969 * std::pow(x, 0.5) - 0.126 * x -
657  0.3516 * std::pow(x, 2) + 0.2843 * std::pow(x, 3) -
658  0.1036 * std::pow(x, 4)); // half thickness at a position x
659 
660  if (is_upper)
661  naca_points.emplace_back(x, +y_t);
662  else
663  naca_points.emplace_back(x, -y_t);
664  }
665  else // is asymmetric
666  for (unsigned int i = 0; i < number_points; i++)
667  {
668  const double m = 1.0 * digit_0 / 100; // max. chamber
669  const double p = 1.0 * digit_1 / 10; // location of max. chamber
670  const double x = i * 1 / (1.0 * number_points - 1);
671 
672  const double y_c =
673  (x <= p) ? m / std::pow(p, 2) * (2 * p * x - std::pow(x, 2)) :
674  m / std::pow(1 - p, 2) *
675  ((1 - 2 * p) + 2 * p * x - std::pow(x, 2));
676 
677  const double dy_c = (x <= p) ?
678  2 * m / std::pow(p, 2) * (p - x) :
679  2 * m / std::pow(1 - p, 2) * (p - x);
680 
681  const double y_t =
682  5 * t *
683  (0.2969 * std::pow(x, 0.5) - 0.126 * x -
684  0.3516 * std::pow(x, 2) + 0.2843 * std::pow(x, 3) -
685  0.1036 * std::pow(x, 4)); // half thicknes at a position x
686 
687  const double theta = std::atan(dy_c);
688 
689  if (is_upper)
690  naca_points.emplace_back(x - y_t * std::sin(theta),
691  y_c + y_t * std::cos(theta));
692  else
693  naca_points.emplace_back(x + y_t * std::sin(theta),
694  y_c - y_t * std::cos(theta));
695  }
696 
697  return naca_points;
698  }
699 
700 
701 
710  static std::array<std::vector<Point<2>>, 2>
711  set_airfoil_length(const std::array<std::vector<Point<2>>, 2> &input,
712  const double desired_len)
713  {
714  std::array<std::vector<Point<2>>, 2> output;
715  output[0] = set_airfoil_length(input[0], desired_len);
716  output[1] = set_airfoil_length(input[1], desired_len);
717 
718  return output;
719  }
720 
728  static std::vector<Point<2>>
729  set_airfoil_length(const std::vector<Point<2>> &input,
730  const double desired_len)
731  {
732  std::vector<Point<2>> output = input;
733 
734  const double scale =
735  desired_len / input.front().distance(input.back());
736 
737  for (auto &x : output)
738  x *= scale;
739 
740  return output;
741  }
742 
753  static std::vector<Point<2>>
754  make_points_equidistant(
755  const std::vector<Point<2>> &non_equidistant_points,
756  const unsigned int number_points)
757  {
758  const unsigned int n_points =
759  non_equidistant_points
760  .size(); // number provided airfoilpoints to interpolate
761 
762  // calculate arclength
763  std::vector<double> arclength_L(non_equidistant_points.size(), 0);
764  for (unsigned int i = 0; i < non_equidistant_points.size() - 1; i++)
765  arclength_L[i + 1] =
766  arclength_L[i] +
767  non_equidistant_points[i + 1].distance(non_equidistant_points[i]);
768 
769 
770  const auto airfoil_length =
771  arclength_L.back(); // arclength upper or lower side
772  const auto deltaX = airfoil_length / (number_points - 1);
773 
774  // Create equidistant points: keep the first (and last) point
775  // unchanged
776  std::vector<Point<2>> equidist(
777  number_points); // number_points is required points on each side for
778  // mesh
779  equidist[0] = non_equidistant_points[0];
780  equidist[number_points - 1] = non_equidistant_points[n_points - 1];
781 
782 
783  // loop over all subsections
784  for (unsigned int j = 0, i = 1; j < n_points - 1; j++)
785  {
786  // get reference left and right end of this section
787  const auto Lj = arclength_L[j];
788  const auto Ljp = arclength_L[j + 1];
789 
790  while (Lj <= i * deltaX && i * deltaX <= Ljp &&
791  i < number_points - 1)
792  {
793  equidist[i] = Point<2>((i * deltaX - Lj) / (Ljp - Lj) *
794  (non_equidistant_points[j + 1] -
795  non_equidistant_points[j]) +
796  non_equidistant_points[j]);
797  ++i;
798  }
799  }
800  return equidist;
801  }
802 
803 
804 
811  void make_coarse_grid(Triangulation<2> &tria) const
812  {
813  // create vector of serial triangulations for each block and
814  // temporary storage for merging them
815  std::vector<Triangulation<2>> trias(10);
816 
817  // helper function to create a subdivided quadrilateral
818  auto make = [](Triangulation<2> & tria,
819  const std::vector<Point<2>> & corner_vertices,
820  const std::vector<unsigned int> &repetitions,
821  const unsigned int material_id) {
822  // create subdivided rectangle with corner points (-1,-1)
823  // and (+1, +1). It serves as reference system
825  repetitions,
826  {-1, -1},
827  {+1, +1});
828 
829  // move all vertices to the correct position
830  for (auto it = tria.begin_vertex(); it < tria.end_vertex(); ++it)
831  {
832  auto & point = it->vertex();
833  const double xi = point(0);
834  const double eta = point(1);
835 
836  // bilinear mapping
837  point = 0.25 * ((1 - xi) * (1 - eta) * corner_vertices[0] +
838  (1 + xi) * (1 - eta) * corner_vertices[1] +
839  (1 - xi) * (1 + eta) * corner_vertices[2] +
840  (1 + xi) * (1 + eta) * corner_vertices[3]);
841  }
842 
843  // set material id of block
844  for (auto cell : tria.active_cell_iterators())
845  cell->set_material_id(material_id);
846  };
847 
848  // create a subdivided quadrilateral for each block (see last number
849  // of block id)
850  make(trias[0],
851  {A, B, D, C},
852  {n_subdivision_y, n_subdivision_x_0},
853  id_block_1);
854  make(trias[1],
855  {F, E, A, B},
856  {n_subdivision_y, n_subdivision_x_0},
857  id_block_4);
858  make(trias[2],
859  {C, H, D, G},
860  {n_subdivision_x_1, n_subdivision_y},
861  id_block_2);
862  make(trias[3],
863  {F, I, E, H},
864  {n_subdivision_x_1, n_subdivision_y},
865  id_block_5);
866  make(trias[4],
867  {H, J, G, K},
868  {n_subdivision_x_2, n_subdivision_y},
869  id_block_3);
870  make(trias[5],
871  {I, L, H, J},
872  {n_subdivision_x_2, n_subdivision_y},
873  id_block_6);
874 
875 
876  // merge triangulation (warning: do not change the order here since
877  // this might change the face ids)
878  GridGenerator::merge_triangulations(trias[0], trias[1], trias[6]);
879  GridGenerator::merge_triangulations(trias[2], trias[3], trias[7]);
880  GridGenerator::merge_triangulations(trias[4], trias[5], trias[8]);
881  GridGenerator::merge_triangulations(trias[6], trias[7], trias[9]);
882  GridGenerator::merge_triangulations(trias[8], trias[9], tria);
883  }
884 
885  /*
886  * Loop over all (cells and) boundary faces of a given triangulation
887  * and set the boundary_ids depending on the material_id of the cell and
888  * the face number. The resulting boundary_ids are:
889  * - 0: inlet
890  * - 1: outlet
891  * - 2: upper airfoil surface (aka. suction side)
892  * - 3, lower airfoil surface (aka. pressure side),
893  * - 4: upper far-field side
894  * - 5: lower far-field side
895  */
896  static void set_boundary_ids(Triangulation<2> &tria)
897  {
898  for (auto cell : tria.active_cell_iterators())
899  for (unsigned int f : GeometryInfo<2>::face_indices())
900  {
901  if (cell->face(f)->at_boundary() == false)
902  continue;
903 
904  const auto mid = cell->material_id();
905 
906  if ((mid == id_block_1 && f == 0) ||
907  (mid == id_block_4 && f == 0))
908  cell->face(f)->set_boundary_id(0); // inlet
909  else if ((mid == id_block_3 && f == 0) ||
910  (mid == id_block_6 && f == 2))
911  cell->face(f)->set_boundary_id(1); // outlet
912  else if ((mid == id_block_1 && f == 1) ||
913  (mid == id_block_2 && f == 1))
914  cell->face(f)->set_boundary_id(2); // upper airfoil side
915  else if ((mid == id_block_4 && f == 1) ||
916  (mid == id_block_5 && f == 3))
917  cell->face(f)->set_boundary_id(3); // lower airfoil side
918  else if ((mid == id_block_2 && f == 0) ||
919  (mid == id_block_3 && f == 2))
920  cell->face(f)->set_boundary_id(4); // upper far-field side
921  else if ((mid == id_block_5 && f == 2) ||
922  (mid == id_block_6 && f == 0))
923  cell->face(f)->set_boundary_id(5); // lower far-field side
924  else
925  Assert(false, ExcIndexRange(mid, id_block_1, id_block_6));
926  }
927  }
928 
929  /*
930  * Interpolate all vertices of the given triangulation onto the airfoil
931  * geometry, depending on the material_id of the block.
932  * Due to symmetry of coarse grid in respect to
933  * x-axis (by definition of points A-L), blocks 1&4, 2&4 and 3&6 can be
934  * interpolated with the same geometric computations Consider a
935  * bias_factor and incline_factor during interpolation to obtain a more
936  * dense mesh next to airfoil geometry and recieve an inclined boundary
937  * between block 2&3 and 5&6, respectively
938  */
939  void interpolate(Triangulation<2> &tria) const
940  {
941  // array storing the information if a vertex was processed
942  std::vector<bool> vertex_processed(tria.n_vertices(), false);
943 
944  // rotation matrix for clockwise rotation of block 1 by angle gamma
945  Tensor<2, 2, double> rotation_matrix_1, rotation_matrix_2;
946 
947  rotation_matrix_1[0][0] = +std::cos(-gamma);
948  rotation_matrix_1[0][1] = -std::sin(-gamma);
949  rotation_matrix_1[1][0] = +std::sin(-gamma);
950  rotation_matrix_1[1][1] = +std::cos(-gamma);
951 
952  rotation_matrix_2 = transpose(rotation_matrix_1);
953 
954  // horizontal offset in order to place coarse-grid node A in the
955  // origin
956  const Point<2, double> horizontal_offset(A(0), 0.0);
957 
958  // Move block 1 so that face BC coincides the x-axis
959  const Point<2, double> trapeze_offset(0.0,
960  std::sin(gamma) * edge_length);
961 
962  // loop over vertices of all cells
963  for (auto &cell : tria)
964  for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell;
965  ++v)
966  {
967  // vertex has been already processed: nothing to do
968  if (vertex_processed[cell.vertex_index(v)])
969  continue;
970 
971  // mark vertex as processed
972  vertex_processed[cell.vertex_index(v)] = true;
973 
974  auto &node = cell.vertex(v);
975 
976  // distinguish blocks
977  if (cell.material_id() == id_block_1 ||
978  cell.material_id() == id_block_4) // block 1 and 4
979  {
980  // step 1: rotate block 1 clockwise by gamma and move block
981  // 1 so that A(0) is on y-axis so that faces AD and BC are
982  // horizontal. This simplifies the computation of the
983  // required indices for interpolation (all x-nodes are
984  // positiv) Move trapeze to be in first quadrant by addig
985  // trapeze_offset
986  Point<2, double> node_;
987  if (cell.material_id() == id_block_1)
988  {
989  node_ = Point<2, double>(rotation_matrix_1 *
990  (node - horizontal_offset) +
991  trapeze_offset);
992  }
993  // step 1: rotate block 4 counterclockwise and move down so
994  // that trapeze is located in fourth quadrant (subtracting
995  // trapeze_offset)
996  else if (cell.material_id() == id_block_4)
997  {
998  node_ = Point<2, double>(rotation_matrix_2 *
999  (node - horizontal_offset) -
1000  trapeze_offset);
1001  }
1002  // step 2: compute indices ix and iy and interpolate
1003  // trapezoid to a rectangle of length pi/2.
1004  {
1005  const double trapeze_height =
1006  std::sin(gamma) * edge_length;
1007  const double L = height / std::sin(gamma);
1008  const double l_a = std::cos(gamma) * edge_length;
1009  const double l_b = trapeze_height * std::tan(gamma);
1010  const double x1 = std::abs(node_(1)) / std::tan(gamma);
1011  const double x2 = L - l_a - l_b;
1012  const double x3 = std::abs(node_(1)) * std::tan(gamma);
1013  const double Dx = x1 + x2 + x3;
1014  const double deltax =
1015  (trapeze_height - std::abs(node_(1))) / std::tan(gamma);
1016  const double dx = Dx / n_cells_x_0;
1017  const double dy = trapeze_height / n_cells_y;
1018  const int ix =
1019  static_cast<int>(std::round((node_(0) - deltax) / dx));
1020  const int iy =
1021  static_cast<int>(std::round(std::abs(node_(1)) / dy));
1022 
1023  node_(0) = numbers::PI / 2 * (1.0 * ix) / n_cells_x_0;
1024  node_(1) = height * (1.0 * iy) / n_cells_y;
1025  }
1026 
1027  // step 3: Interpolation between semicircle (of C-Mesh) and
1028  // airfoil contour
1029  {
1030  const double dx = numbers::PI / 2 / n_cells_x_0;
1031  const double dy = height / n_cells_y;
1032  const int ix =
1033  static_cast<int>(std::round(node_(0) / dx));
1034  const int iy =
1035  static_cast<int>(std::round(node_(1) / dy));
1036  const double alpha =
1037  bias_alpha(1 - (1.0 * iy) / n_cells_y);
1038  const double theta = node_(0);
1039  const Point<2> p(-height * std::cos(theta) + center_mesh,
1040  ((cell.material_id() == id_block_1) ?
1041  (height) :
1042  (-height)) *
1043  std::sin(theta));
1044  node =
1045  airfoil_1D[(
1046  (cell.material_id() == id_block_1) ? (0) : (1))][ix] *
1047  alpha +
1048  p * (1 - alpha);
1049  }
1050  }
1051  else if (cell.material_id() == id_block_2 ||
1052  cell.material_id() == id_block_5) // block 2 and 5
1053  {
1054  // geometric parameters and indices for interpolation
1055  Assert(
1056  (std::abs(D(1) - C(1)) == std::abs(F(1) - E(1))) &&
1057  (std::abs(C(1)) == std::abs(E(1))) &&
1058  (std::abs(G(1)) == std::abs(I(1))),
1059  ExcMessage(
1060  "Points D,C,G and E,F,I are not defined symmetric to "
1061  "x-axis, which is required to interpolate block 2"
1062  " and 5 with same geometric computations."));
1063  const double l_y = D(1) - C(1);
1064  const double l_h = D(1) - l_y;
1065  const double by = -l_h / length_b1_x * (node(0) - H(0));
1066  const double dy = (height - by) / n_cells_y;
1067  const int iy = static_cast<int>(
1068  std::round((std::abs(node(1)) - by) / dy));
1069  const double dx = length_b1_x / n_cells_x_1;
1070  const int ix = static_cast<int>(
1071  std::round(std::abs(node(0) - center_mesh) / dx));
1072 
1073  const double alpha = bias_alpha(1 - (1.0 * iy) / n_cells_y);
1074  // define points on upper/lower horizontal far field side,
1075  // i.e. face DG or FI. Incline factor to move points G and I
1076  // to the right by distance incline_facor*lenght_b2
1077  const Point<2> p(ix * dx + center_mesh +
1078  incline_factor * length_b2 * ix /
1079  n_cells_x_1,
1080  ((cell.material_id() == id_block_2) ?
1081  (height) :
1082  (-height)));
1083  // interpolate between y = height and upper airfoil points
1084  // (block2) or y = -height and lower airfoil points (block5)
1085  node = airfoil_1D[(
1086  (cell.material_id() == id_block_2) ? (0) : (1))]
1087  [n_cells_x_0 + ix] *
1088  alpha +
1089  p * (1 - alpha);
1090  }
1091  else if (cell.material_id() == id_block_3 ||
1092  cell.material_id() == id_block_6) // block 3 and 6
1093  {
1094  // compute indices ix and iy
1095  const double dx = length_b2 / n_cells_x_2;
1096  const double dy = height / n_cells_y;
1097  const int ix = static_cast<int>(
1098  std::round(std::abs(node(0) - H(0)) / dx));
1099  const int iy =
1100  static_cast<int>(std::round(std::abs(node(1)) / dy));
1101 
1102  const double alpha_y = bias_alpha(1 - 1.0 * iy / n_cells_y);
1103  const double alpha_x =
1104  bias_alpha(1 - (static_cast<double>(ix)) / n_cells_x_2);
1105  // define on upper/lower horizontal far field side at y =
1106  // +/- height, i.e. face GK or IL incline factor to move
1107  // points G and H to the right
1108  const Point<2> p1(J(0) - (1 - incline_factor) * length_b2 *
1109  (alpha_x),
1110  ((cell.material_id() == id_block_3) ?
1111  (height) :
1112  (-height)));
1113  // define points on HJ but use tail_y as y-coordinate, in
1114  // case last airfoil point has y =/= 0
1115  const Point<2> p2(J(0) - alpha_x * length_b2, tail_y);
1116  node = p1 * (1 - alpha_y) + p2 * alpha_y;
1117  }
1118  else
1119  {
1120  Assert(false,
1121  ExcIndexRange(cell.material_id(),
1122  id_block_1,
1123  id_block_6));
1124  }
1125  }
1126  }
1127 
1128 
1129  /*
1130  * This function returns a bias factor 'alpha' which is used to make the
1131  * mesh more tight in close distance of the airfoil.
1132  * It is a bijective function mapping from [0,1] onto [0,1] where values
1133  * near 1 are made tighter.
1134  */
1135  double
1136  bias_alpha(double alpha) const
1137  {
1138  return std::tanh(bias_factor * alpha) / std::tanh(bias_factor);
1139  }
1140  };
1141  } // namespace
1142 
1143 
1144 
1145  void internal_create_triangulation(
1146  Triangulation<2, 2> & tria,
1147  std::vector<GridTools::PeriodicFacePair<
1148  typename Triangulation<2, 2>::cell_iterator>> *periodic_faces,
1149  const AdditionalData & additional_data)
1150  {
1151  MeshGenerator mesh_generator(additional_data);
1152  // Cast the the triangulation to the right type so that the right
1153  // specialization of the function create_triangulation is picked up.
1154  if (auto parallel_tria =
1156  &tria))
1157  mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1158  else if (auto parallel_tria = dynamic_cast<
1160  &tria))
1161  mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1162  else
1163  mesh_generator.create_triangulation(tria, periodic_faces);
1164  }
1165 
1166  template <>
1168  {
1169  Assert(false, ExcMessage("Airfoils only exist for 2D and 3D!"));
1170  }
1171 
1172 
1173 
1174  template <>
1176  std::vector<GridTools::PeriodicFacePair<
1178  const AdditionalData &)
1179  {
1180  Assert(false, ExcMessage("Airfoils only exist for 2D and 3D!"));
1181  }
1182 
1183 
1184 
1185  template <>
1187  const AdditionalData &additional_data)
1188  {
1189  internal_create_triangulation(tria, nullptr, additional_data);
1190  }
1191 
1192 
1193 
1194  template <>
1195  void create_triangulation(
1196  Triangulation<2, 2> & tria,
1197  std::vector<GridTools::PeriodicFacePair<
1198  typename Triangulation<2, 2>::cell_iterator>> &periodic_faces,
1199  const AdditionalData & additional_data)
1200  {
1201  internal_create_triangulation(tria, &periodic_faces, additional_data);
1202  }
1203 
1204 
1205 
1206  template <>
1207  void create_triangulation(
1208  Triangulation<3, 3> & tria,
1209  std::vector<GridTools::PeriodicFacePair<
1210  typename Triangulation<3, 3>::cell_iterator>> &periodic_faces,
1211  const AdditionalData & additional_data)
1212  {
1213  Assert(false, ExcMessage("3D airfoils are not implemented yet!"));
1214  (void)tria;
1215  (void)additional_data;
1216  (void)periodic_faces;
1217  }
1218  } // namespace Airfoil
1219 
1220 
1221  namespace
1222  {
1227  template <int dim, int spacedim>
1228  void
1229  colorize_hyper_rectangle(Triangulation<dim, spacedim> &tria)
1230  {
1231  // there is nothing to do in 1d
1232  if (dim > 1)
1233  {
1234  // there is only one cell, so
1235  // simple task
1236  const typename Triangulation<dim, spacedim>::cell_iterator cell =
1237  tria.begin();
1238  for (auto f : GeometryInfo<dim>::face_indices())
1239  cell->face(f)->set_boundary_id(f);
1240  }
1241  }
1242 
1243 
1244 
1245  template <int spacedim>
1246  void colorize_subdivided_hyper_rectangle(Triangulation<1, spacedim> &tria,
1247  const Point<spacedim> &,
1248  const Point<spacedim> &,
1249  const double)
1250  {
1251  for (typename Triangulation<1, spacedim>::cell_iterator cell =
1252  tria.begin();
1253  cell != tria.end();
1254  ++cell)
1255  if (cell->center()(0) > 0)
1256  cell->set_material_id(1);
1257  // boundary indicators are set to
1258  // 0 (left) and 1 (right) by default.
1259  }
1260 
1261 
1262 
1263  template <int dim, int spacedim>
1264  void
1265  colorize_subdivided_hyper_rectangle(Triangulation<dim, spacedim> &tria,
1266  const Point<spacedim> & p1,
1267  const Point<spacedim> & p2,
1268  const double epsilon)
1269  {
1270  // run through all faces and check
1271  // if one of their center coordinates matches
1272  // one of the corner points. Comparisons
1273  // are made using an epsilon which
1274  // should be smaller than the smallest cell
1275  // diameter.
1276 
1278  tria.begin_face(),
1279  endface =
1280  tria.end_face();
1281  for (; face != endface; ++face)
1282  if (face->at_boundary())
1283  if (face->boundary_id() == 0)
1284  {
1285  const Point<spacedim> center(face->center());
1286 
1287  if (std::abs(center(0) - p1[0]) < epsilon)
1288  face->set_boundary_id(0);
1289  else if (std::abs(center(0) - p2[0]) < epsilon)
1290  face->set_boundary_id(1);
1291  else if (dim > 1 && std::abs(center(1) - p1[1]) < epsilon)
1292  face->set_boundary_id(2);
1293  else if (dim > 1 && std::abs(center(1) - p2[1]) < epsilon)
1294  face->set_boundary_id(3);
1295  else if (dim > 2 && std::abs(center(2) - p1[2]) < epsilon)
1296  face->set_boundary_id(4);
1297  else if (dim > 2 && std::abs(center(2) - p2[2]) < epsilon)
1298  face->set_boundary_id(5);
1299  else
1300  // triangulation says it
1301  // is on the boundary,
1302  // but we could not find
1303  // on which boundary.
1304  Assert(false, ExcInternalError());
1305  }
1306 
1307  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
1308  tria.begin();
1309  cell != tria.end();
1310  ++cell)
1311  {
1312  char id = 0;
1313  for (unsigned int d = 0; d < dim; ++d)
1314  if (cell->center()(d) > 0)
1315  id += (1 << d);
1316  cell->set_material_id(id);
1317  }
1318  }
1319 
1320 
1325  void colorize_hyper_shell(Triangulation<2> &tria,
1326  const Point<2> &,
1327  const double,
1328  const double)
1329  {
1330  // In spite of receiving geometrical
1331  // data, we do this only based on
1332  // topology.
1333 
1334  // For the mesh based on cube,
1335  // this is highly irregular
1336  for (Triangulation<2>::cell_iterator cell = tria.begin();
1337  cell != tria.end();
1338  ++cell)
1339  {
1340  Assert(cell->face(2)->at_boundary(), ExcInternalError());
1341  cell->face(2)->set_all_boundary_ids(1);
1342  }
1343  }
1344 
1345 
1350  void colorize_hyper_shell(Triangulation<3> &tria,
1351  const Point<3> &,
1352  const double,
1353  const double)
1354  {
1355  // the following uses a good amount
1356  // of knowledge about the
1357  // orientation of cells. this is
1358  // probably not good style...
1359  if (tria.n_cells() == 6)
1360  {
1361  Triangulation<3>::cell_iterator cell = tria.begin();
1362 
1363  Assert(cell->face(4)->at_boundary(), ExcInternalError());
1364  cell->face(4)->set_all_boundary_ids(1);
1365 
1366  ++cell;
1367  Assert(cell->face(2)->at_boundary(), ExcInternalError());
1368  cell->face(2)->set_all_boundary_ids(1);
1369 
1370  ++cell;
1371  Assert(cell->face(2)->at_boundary(), ExcInternalError());
1372  cell->face(2)->set_all_boundary_ids(1);
1373 
1374  ++cell;
1375  Assert(cell->face(0)->at_boundary(), ExcInternalError());
1376  cell->face(0)->set_all_boundary_ids(1);
1377 
1378  ++cell;
1379  Assert(cell->face(2)->at_boundary(), ExcInternalError());
1380  cell->face(2)->set_all_boundary_ids(1);
1381 
1382  ++cell;
1383  Assert(cell->face(0)->at_boundary(), ExcInternalError());
1384  cell->face(0)->set_all_boundary_ids(1);
1385  }
1386  else if (tria.n_cells() == 12)
1387  {
1388  // again use some internal
1389  // knowledge
1390  for (Triangulation<3>::cell_iterator cell = tria.begin();
1391  cell != tria.end();
1392  ++cell)
1393  {
1394  Assert(cell->face(5)->at_boundary(), ExcInternalError());
1395  cell->face(5)->set_all_boundary_ids(1);
1396  }
1397  }
1398  else if (tria.n_cells() == 96)
1399  {
1400  // the 96-cell hypershell is
1401  // based on a once refined
1402  // 12-cell mesh. consequently,
1403  // since the outer faces all
1404  // are face_no==5 above, so
1405  // they are here (unless they
1406  // are in the interior). Use
1407  // this to assign boundary
1408  // indicators, but also make
1409  // sure that we encounter
1410  // exactly 48 such faces
1411  unsigned int count = 0;
1412  for (Triangulation<3>::cell_iterator cell = tria.begin();
1413  cell != tria.end();
1414  ++cell)
1415  if (cell->face(5)->at_boundary())
1416  {
1417  cell->face(5)->set_all_boundary_ids(1);
1418  ++count;
1419  }
1420  Assert(count == 48, ExcInternalError());
1421  }
1422  else
1423  Assert(false, ExcNotImplemented());
1424  }
1425 
1426 
1427 
1433  void colorize_quarter_hyper_shell(Triangulation<3> &tria,
1434  const Point<3> & center,
1435  const double inner_radius,
1436  const double outer_radius)
1437  {
1438  if (tria.n_cells() != 3)
1439  AssertThrow(false, ExcNotImplemented());
1440 
1441  double middle = (outer_radius - inner_radius) / 2e0 + inner_radius;
1442  double eps = 1e-3 * middle;
1443  Triangulation<3>::cell_iterator cell = tria.begin();
1444 
1445  for (; cell != tria.end(); ++cell)
1446  for (unsigned int f : GeometryInfo<3>::face_indices())
1447  {
1448  if (!cell->face(f)->at_boundary())
1449  continue;
1450 
1451  double radius = cell->face(f)->center().norm() - center.norm();
1452  if (std::fabs(cell->face(f)->center()(0)) <
1453  eps) // x = 0 set boundary 2
1454  {
1455  cell->face(f)->set_boundary_id(2);
1456  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1457  ++j)
1458  if (cell->face(f)->line(j)->at_boundary())
1459  if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1460  cell->face(f)->line(j)->vertex(1).norm()) >
1461  eps)
1462  cell->face(f)->line(j)->set_boundary_id(2);
1463  }
1464  else if (std::fabs(cell->face(f)->center()(1)) <
1465  eps) // y = 0 set boundary 3
1466  {
1467  cell->face(f)->set_boundary_id(3);
1468  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1469  ++j)
1470  if (cell->face(f)->line(j)->at_boundary())
1471  if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1472  cell->face(f)->line(j)->vertex(1).norm()) >
1473  eps)
1474  cell->face(f)->line(j)->set_boundary_id(3);
1475  }
1476  else if (std::fabs(cell->face(f)->center()(2)) <
1477  eps) // z = 0 set boundary 4
1478  {
1479  cell->face(f)->set_boundary_id(4);
1480  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1481  ++j)
1482  if (cell->face(f)->line(j)->at_boundary())
1483  if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1484  cell->face(f)->line(j)->vertex(1).norm()) >
1485  eps)
1486  cell->face(f)->line(j)->set_boundary_id(4);
1487  }
1488  else if (radius < middle) // inner radius set boundary 0
1489  {
1490  cell->face(f)->set_boundary_id(0);
1491  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1492  ++j)
1493  if (cell->face(f)->line(j)->at_boundary())
1494  if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1495  cell->face(f)->line(j)->vertex(1).norm()) <
1496  eps)
1497  cell->face(f)->line(j)->set_boundary_id(0);
1498  }
1499  else if (radius > middle) // outer radius set boundary 1
1500  {
1501  cell->face(f)->set_boundary_id(1);
1502  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1503  ++j)
1504  if (cell->face(f)->line(j)->at_boundary())
1505  if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1506  cell->face(f)->line(j)->vertex(1).norm()) <
1507  eps)
1508  cell->face(f)->line(j)->set_boundary_id(1);
1509  }
1510  else
1511  Assert(false, ExcInternalError());
1512  }
1513  }
1514 
1515  } // namespace
1516 
1517 
1518  template <int dim, int spacedim>
1519  void
1521  const Point<dim> & p_1,
1522  const Point<dim> & p_2,
1523  const bool colorize)
1524  {
1525  // First, extend dimensions from dim to spacedim and
1526  // normalize such that p1 is lower in all coordinate
1527  // directions. Additional entries will be 0.
1528  Point<spacedim> p1, p2;
1529  for (unsigned int i = 0; i < dim; ++i)
1530  {
1531  p1(i) = std::min(p_1(i), p_2(i));
1532  p2(i) = std::max(p_1(i), p_2(i));
1533  }
1534 
1535  std::vector<Point<spacedim>> vertices(GeometryInfo<dim>::vertices_per_cell);
1536  switch (dim)
1537  {
1538  case 1:
1539  vertices[0] = p1;
1540  vertices[1] = p2;
1541  break;
1542  case 2:
1543  vertices[0] = vertices[1] = p1;
1544  vertices[2] = vertices[3] = p2;
1545 
1546  vertices[1](0) = p2(0);
1547  vertices[2](0) = p1(0);
1548  break;
1549  case 3:
1550  vertices[0] = vertices[1] = vertices[2] = vertices[3] = p1;
1551  vertices[4] = vertices[5] = vertices[6] = vertices[7] = p2;
1552 
1553  vertices[1](0) = p2(0);
1554  vertices[2](1) = p2(1);
1555  vertices[3](0) = p2(0);
1556  vertices[3](1) = p2(1);
1557 
1558  vertices[4](0) = p1(0);
1559  vertices[4](1) = p1(1);
1560  vertices[5](1) = p1(1);
1561  vertices[6](0) = p1(0);
1562 
1563  break;
1564  default:
1565  Assert(false, ExcNotImplemented());
1566  }
1567 
1568  // Prepare cell data
1569  std::vector<CellData<dim>> cells(1);
1570  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
1571  cells[0].vertices[i] = i;
1572  cells[0].material_id = 0;
1573 
1574  tria.create_triangulation(vertices, cells, SubCellData());
1575 
1576  // Assign boundary indicators
1577  if (colorize)
1578  colorize_hyper_rectangle(tria);
1579  }
1580 
1581 
1582  template <int dim, int spacedim>
1583  void
1585  const double left,
1586  const double right,
1587  const bool colorize)
1588  {
1589  Assert(left < right,
1590  ExcMessage("Invalid left-to-right bounds of hypercube"));
1591 
1592  Point<dim> p1, p2;
1593  for (unsigned int i = 0; i < dim; ++i)
1594  {
1595  p1(i) = left;
1596  p2(i) = right;
1597  }
1598  hyper_rectangle(tria, p1, p2, colorize);
1599  }
1600 
1601  template <int dim>
1602  void
1603  simplex(Triangulation<dim> &tria, const std::vector<Point<dim>> &vertices)
1604  {
1605  AssertDimension(vertices.size(), dim + 1);
1606  Assert(dim > 1, ExcNotImplemented());
1607  Assert(dim < 4, ExcNotImplemented());
1608 
1609 #ifdef DEBUG
1610  Tensor<2, dim> vector_matrix;
1611  for (unsigned int d = 0; d < dim; ++d)
1612  for (unsigned int c = 1; c <= dim; ++c)
1613  vector_matrix[c - 1][d] = vertices[c](d) - vertices[0](d);
1614  Assert(determinant(vector_matrix) > 0.,
1615  ExcMessage("Vertices of simplex must form a right handed system"));
1616 #endif
1617 
1618  // Set up the vertices by first copying into points.
1619  std::vector<Point<dim>> points = vertices;
1620  Point<dim> center;
1621  // Compute the edge midpoints and add up everything to compute the
1622  // center point.
1623  for (unsigned int i = 0; i <= dim; ++i)
1624  {
1625  points.push_back(0.5 * (points[i] + points[(i + 1) % (dim + 1)]));
1626  center += points[i];
1627  }
1628  if (dim > 2)
1629  {
1630  // In 3D, we have some more edges to deal with
1631  for (unsigned int i = 1; i < dim; ++i)
1632  points.push_back(0.5 * (points[i - 1] + points[i + 1]));
1633  // And we need face midpoints
1634  for (unsigned int i = 0; i <= dim; ++i)
1635  points.push_back(1. / 3. *
1636  (points[i] + points[(i + 1) % (dim + 1)] +
1637  points[(i + 2) % (dim + 1)]));
1638  }
1639  points.push_back((1. / (dim + 1)) * center);
1640 
1641  std::vector<CellData<dim>> cells(dim + 1);
1642  switch (dim)
1643  {
1644  case 2:
1645  AssertDimension(points.size(), 7);
1646  cells[0].vertices[0] = 0;
1647  cells[0].vertices[1] = 3;
1648  cells[0].vertices[2] = 5;
1649  cells[0].vertices[3] = 6;
1650  cells[0].material_id = 0;
1651 
1652  cells[1].vertices[0] = 3;
1653  cells[1].vertices[1] = 1;
1654  cells[1].vertices[2] = 6;
1655  cells[1].vertices[3] = 4;
1656  cells[1].material_id = 0;
1657 
1658  cells[2].vertices[0] = 5;
1659  cells[2].vertices[1] = 6;
1660  cells[2].vertices[2] = 2;
1661  cells[2].vertices[3] = 4;
1662  cells[2].material_id = 0;
1663  break;
1664  case 3:
1665  AssertDimension(points.size(), 15);
1666  cells[0].vertices[0] = 0;
1667  cells[0].vertices[1] = 4;
1668  cells[0].vertices[2] = 8;
1669  cells[0].vertices[3] = 10;
1670  cells[0].vertices[4] = 7;
1671  cells[0].vertices[5] = 13;
1672  cells[0].vertices[6] = 12;
1673  cells[0].vertices[7] = 14;
1674  cells[0].material_id = 0;
1675 
1676  cells[1].vertices[0] = 4;
1677  cells[1].vertices[1] = 1;
1678  cells[1].vertices[2] = 10;
1679  cells[1].vertices[3] = 5;
1680  cells[1].vertices[4] = 13;
1681  cells[1].vertices[5] = 9;
1682  cells[1].vertices[6] = 14;
1683  cells[1].vertices[7] = 11;
1684  cells[1].material_id = 0;
1685 
1686  cells[2].vertices[0] = 8;
1687  cells[2].vertices[1] = 10;
1688  cells[2].vertices[2] = 2;
1689  cells[2].vertices[3] = 5;
1690  cells[2].vertices[4] = 12;
1691  cells[2].vertices[5] = 14;
1692  cells[2].vertices[6] = 6;
1693  cells[2].vertices[7] = 11;
1694  cells[2].material_id = 0;
1695 
1696  cells[3].vertices[0] = 7;
1697  cells[3].vertices[1] = 13;
1698  cells[3].vertices[2] = 12;
1699  cells[3].vertices[3] = 14;
1700  cells[3].vertices[4] = 3;
1701  cells[3].vertices[5] = 9;
1702  cells[3].vertices[6] = 6;
1703  cells[3].vertices[7] = 11;
1704  cells[3].material_id = 0;
1705  break;
1706  default:
1707  Assert(false, ExcNotImplemented());
1708  }
1709  tria.create_triangulation(points, cells, SubCellData());
1710  }
1711 
1712 
1713  void moebius(Triangulation<3> & tria,
1714  const unsigned int n_cells,
1715  const unsigned int n_rotations,
1716  const double R,
1717  const double r)
1718  {
1719  const unsigned int dim = 3;
1720  Assert(n_cells > 4,
1721  ExcMessage(
1722  "More than 4 cells are needed to create a moebius grid."));
1723  Assert(r > 0 && R > 0,
1724  ExcMessage("Outer and inner radius must be positive."));
1725  Assert(R > r,
1726  ExcMessage("Outer radius must be greater than inner radius."));
1727 
1728 
1729  std::vector<Point<dim>> vertices(4 * n_cells);
1730  double beta_step = n_rotations * numbers::PI / 2.0 / n_cells;
1731  double alpha_step = 2.0 * numbers::PI / n_cells;
1732 
1733  for (unsigned int i = 0; i < n_cells; ++i)
1734  for (unsigned int j = 0; j < 4; ++j)
1735  {
1736  vertices[4 * i + j][0] =
1737  R * std::cos(i * alpha_step) +
1738  r * std::cos(i * beta_step + j * numbers::PI / 2.0) *
1739  std::cos(i * alpha_step);
1740  vertices[4 * i + j][1] =
1741  R * std::sin(i * alpha_step) +
1742  r * std::cos(i * beta_step + j * numbers::PI / 2.0) *
1743  std::sin(i * alpha_step);
1744  vertices[4 * i + j][2] =
1745  r * std::sin(i * beta_step + j * numbers::PI / 2.0);
1746  }
1747 
1748  unsigned int offset = 0;
1749 
1750  std::vector<CellData<dim>> cells(n_cells);
1751  for (unsigned int i = 0; i < n_cells; ++i)
1752  {
1753  for (unsigned int j = 0; j < 2; ++j)
1754  {
1755  cells[i].vertices[0 + 4 * j] = offset + 0 + 4 * j;
1756  cells[i].vertices[1 + 4 * j] = offset + 3 + 4 * j;
1757  cells[i].vertices[2 + 4 * j] = offset + 2 + 4 * j;
1758  cells[i].vertices[3 + 4 * j] = offset + 1 + 4 * j;
1759  }
1760  offset += 4;
1761  cells[i].material_id = 0;
1762  }
1763 
1764  // now correct the last four vertices
1765  cells[n_cells - 1].vertices[4] = (0 + n_rotations) % 4;
1766  cells[n_cells - 1].vertices[5] = (3 + n_rotations) % 4;
1767  cells[n_cells - 1].vertices[6] = (2 + n_rotations) % 4;
1768  cells[n_cells - 1].vertices[7] = (1 + n_rotations) % 4;
1769 
1771  tria.create_triangulation_compatibility(vertices, cells, SubCellData());
1772  }
1773 
1774 
1775 
1776  template <>
1777  void torus<2, 3>(Triangulation<2, 3> &tria,
1778  const double R,
1779  const double r,
1780  const unsigned int,
1781  const double)
1782  {
1783  Assert(R > r,
1784  ExcMessage("Outer radius R must be greater than the inner "
1785  "radius r."));
1786  Assert(r > 0.0, ExcMessage("The inner radius r must be positive."));
1787 
1788  const unsigned int dim = 2;
1789  const unsigned int spacedim = 3;
1790  std::vector<Point<spacedim>> vertices(16);
1791 
1792  vertices[0] = Point<spacedim>(R - r, 0, 0);
1793  vertices[1] = Point<spacedim>(R, -r, 0);
1794  vertices[2] = Point<spacedim>(R + r, 0, 0);
1795  vertices[3] = Point<spacedim>(R, r, 0);
1796  vertices[4] = Point<spacedim>(0, 0, R - r);
1797  vertices[5] = Point<spacedim>(0, -r, R);
1798  vertices[6] = Point<spacedim>(0, 0, R + r);
1799  vertices[7] = Point<spacedim>(0, r, R);
1800  vertices[8] = Point<spacedim>(-(R - r), 0, 0);
1801  vertices[9] = Point<spacedim>(-R, -r, 0);
1802  vertices[10] = Point<spacedim>(-(R + r), 0, 0);
1803  vertices[11] = Point<spacedim>(-R, r, 0);
1804  vertices[12] = Point<spacedim>(0, 0, -(R - r));
1805  vertices[13] = Point<spacedim>(0, -r, -R);
1806  vertices[14] = Point<spacedim>(0, 0, -(R + r));
1807  vertices[15] = Point<spacedim>(0, r, -R);
1808 
1809  std::vector<CellData<dim>> cells(16);
1810  // Right Hand Orientation
1811  cells[0].vertices[0] = 0;
1812  cells[0].vertices[1] = 4;
1813  cells[0].vertices[2] = 7;
1814  cells[0].vertices[3] = 3;
1815  cells[0].material_id = 0;
1816 
1817  cells[1].vertices[0] = 1;
1818  cells[1].vertices[1] = 5;
1819  cells[1].vertices[2] = 4;
1820  cells[1].vertices[3] = 0;
1821  cells[1].material_id = 0;
1822 
1823  cells[2].vertices[0] = 2;
1824  cells[2].vertices[1] = 6;
1825  cells[2].vertices[2] = 5;
1826  cells[2].vertices[3] = 1;
1827  cells[2].material_id = 0;
1828 
1829  cells[3].vertices[0] = 3;
1830  cells[3].vertices[1] = 7;
1831  cells[3].vertices[2] = 6;
1832  cells[3].vertices[3] = 2;
1833  cells[3].material_id = 0;
1834 
1835  cells[4].vertices[0] = 4;
1836  cells[4].vertices[1] = 8;
1837  cells[4].vertices[2] = 11;
1838  cells[4].vertices[3] = 7;
1839  cells[4].material_id = 0;
1840 
1841  cells[5].vertices[0] = 5;
1842  cells[5].vertices[1] = 9;
1843  cells[5].vertices[2] = 8;
1844  cells[5].vertices[3] = 4;
1845  cells[5].material_id = 0;
1846 
1847  cells[6].vertices[0] = 6;
1848  cells[6].vertices[1] = 10;
1849  cells[6].vertices[2] = 9;
1850  cells[6].vertices[3] = 5;
1851  cells[6].material_id = 0;
1852 
1853  cells[7].vertices[0] = 7;
1854  cells[7].vertices[1] = 11;
1855  cells[7].vertices[2] = 10;
1856  cells[7].vertices[3] = 6;
1857  cells[7].material_id = 0;
1858 
1859  cells[8].vertices[0] = 8;
1860  cells[8].vertices[1] = 12;
1861  cells[8].vertices[2] = 15;
1862  cells[8].vertices[3] = 11;
1863  cells[8].material_id = 0;
1864 
1865  cells[9].vertices[0] = 9;
1866  cells[9].vertices[1] = 13;
1867  cells[9].vertices[2] = 12;
1868  cells[9].vertices[3] = 8;
1869  cells[9].material_id = 0;
1870 
1871  cells[10].vertices[0] = 10;
1872  cells[10].vertices[1] = 14;
1873  cells[10].vertices[2] = 13;
1874  cells[10].vertices[3] = 9;
1875  cells[10].material_id = 0;
1876 
1877  cells[11].vertices[0] = 11;
1878  cells[11].vertices[1] = 15;
1879  cells[11].vertices[2] = 14;
1880  cells[11].vertices[3] = 10;
1881  cells[11].material_id = 0;
1882 
1883  cells[12].vertices[0] = 12;
1884  cells[12].vertices[1] = 0;
1885  cells[12].vertices[2] = 3;
1886  cells[12].vertices[3] = 15;
1887  cells[12].material_id = 0;
1888 
1889  cells[13].vertices[0] = 13;
1890  cells[13].vertices[1] = 1;
1891  cells[13].vertices[2] = 0;
1892  cells[13].vertices[3] = 12;
1893  cells[13].material_id = 0;
1894 
1895  cells[14].vertices[0] = 14;
1896  cells[14].vertices[1] = 2;
1897  cells[14].vertices[2] = 1;
1898  cells[14].vertices[3] = 13;
1899  cells[14].material_id = 0;
1900 
1901  cells[15].vertices[0] = 15;
1902  cells[15].vertices[1] = 3;
1903  cells[15].vertices[2] = 2;
1904  cells[15].vertices[3] = 14;
1905  cells[15].material_id = 0;
1906 
1907  // Must call this to be able to create a
1908  // correct triangulation in dealii, read
1909  // GridReordering<> doc
1911  tria.create_triangulation_compatibility(vertices, cells, SubCellData());
1912 
1913  tria.set_all_manifold_ids(0);
1914  tria.set_manifold(0, TorusManifold<2>(R, r));
1915  }
1916 
1917 
1918 
1919  template <>
1920  void torus<3, 3>(Triangulation<3, 3> &tria,
1921  const double R,
1922  const double r,
1923  const unsigned int n_cells_toroidal,
1924  const double phi)
1925  {
1926  Assert(R > r,
1927  ExcMessage("Outer radius R must be greater than the inner "
1928  "radius r."));
1929  Assert(r > 0.0, ExcMessage("The inner radius r must be positive."));
1930  Assert(n_cells_toroidal > 2,
1931  ExcMessage("Number of cells in toroidal direction has "
1932  "to be at least 3."));
1933  AssertThrow(phi > 0.0 && phi < 2.0 * numbers::PI + 1.0e-15,
1934  ExcMessage("Invalid angle phi specified."));
1935 
1936  // the first 8 vertices are in the x-y-plane
1937  Point<3> const p = Point<3>(R, 0.0, 0.0);
1938  double const a = 1. / (1 + std::sqrt(2.0));
1939  // A value of 1 indicates "open" torus with angle < 2*pi, which
1940  // means that we need an additional layer of vertices
1941  const unsigned int additional_layer =
1942  (phi < 2.0 * numbers::PI - 1.0e-15) ?
1943  1 :
1944  0; // torus is closed (angle of 2*pi)
1945  const unsigned int n_point_layers_toroidal =
1946  n_cells_toroidal + additional_layer;
1947  std::vector<Point<3>> vertices(8 * n_point_layers_toroidal);
1948  vertices[0] = p + Point<3>(-1, -1, 0) * (r / std::sqrt(2.0)),
1949  vertices[1] = p + Point<3>(+1, -1, 0) * (r / std::sqrt(2.0)),
1950  vertices[2] = p + Point<3>(-1, -1, 0) * (r / std::sqrt(2.0) * a),
1951  vertices[3] = p + Point<3>(+1, -1, 0) * (r / std::sqrt(2.0) * a),
1952  vertices[4] = p + Point<3>(-1, +1, 0) * (r / std::sqrt(2.0) * a),
1953  vertices[5] = p + Point<3>(+1, +1, 0) * (r / std::sqrt(2.0) * a),
1954  vertices[6] = p + Point<3>(-1, +1, 0) * (r / std::sqrt(2.0)),
1955  vertices[7] = p + Point<3>(+1, +1, 0) * (r / std::sqrt(2.0));
1956 
1957  // create remaining vertices by rotating around negative y-axis (the
1958  // direction is to ensure positive cell measures)
1959  double const phi_cell = phi / n_cells_toroidal;
1960  for (unsigned int c = 1; c < n_point_layers_toroidal; ++c)
1961  {
1962  for (unsigned int v = 0; v < 8; ++v)
1963  {
1964  double const r_2d = vertices[v][0];
1965  vertices[8 * c + v][0] = r_2d * std::cos(phi_cell * c);
1966  vertices[8 * c + v][1] = vertices[v][1];
1967  vertices[8 * c + v][2] = r_2d * std::sin(phi_cell * c);
1968  }
1969  }
1970 
1971  // cell connectivity
1972  std::vector<CellData<3>> cells(5 * n_cells_toroidal);
1973  for (unsigned int c = 0; c < n_cells_toroidal; ++c)
1974  {
1975  for (unsigned int j = 0; j < 2; ++j)
1976  {
1977  const unsigned int offset =
1978  (8 * (c + j)) % (8 * n_point_layers_toroidal);
1979 
1980  // cell 0 in x-y-plane
1981  cells[5 * c].vertices[0 + j * 4] = offset + 0;
1982  cells[5 * c].vertices[1 + j * 4] = offset + 1;
1983  cells[5 * c].vertices[2 + j * 4] = offset + 2;
1984  cells[5 * c].vertices[3 + j * 4] = offset + 3;
1985  // cell 1 in x-y-plane (cell on torus centerline)
1986  cells[5 * c + 1].vertices[0 + j * 4] = offset + 2;
1987  cells[5 * c + 1].vertices[1 + j * 4] = offset + 3;
1988  cells[5 * c + 1].vertices[2 + j * 4] = offset + 4;
1989  cells[5 * c + 1].vertices[3 + j * 4] = offset + 5;
1990  // cell 2 in x-y-plane
1991  cells[5 * c + 2].vertices[0 + j * 4] = offset + 4;
1992  cells[5 * c + 2].vertices[1 + j * 4] = offset + 5;
1993  cells[5 * c + 2].vertices[2 + j * 4] = offset + 6;
1994  cells[5 * c + 2].vertices[3 + j * 4] = offset + 7;
1995  // cell 3 in x-y-plane
1996  cells[5 * c + 3].vertices[0 + j * 4] = offset + 0;
1997  cells[5 * c + 3].vertices[1 + j * 4] = offset + 2;
1998  cells[5 * c + 3].vertices[2 + j * 4] = offset + 6;
1999  cells[5 * c + 3].vertices[3 + j * 4] = offset + 4;
2000  // cell 4 in x-y-plane
2001  cells[5 * c + 4].vertices[0 + j * 4] = offset + 3;
2002  cells[5 * c + 4].vertices[1 + j * 4] = offset + 1;
2003  cells[5 * c + 4].vertices[2 + j * 4] = offset + 5;
2004  cells[5 * c + 4].vertices[3 + j * 4] = offset + 7;
2005  }
2006 
2007  cells[5 * c].material_id = 0;
2008  // mark cell on torus centerline
2009  cells[5 * c + 1].material_id = 1;
2010  cells[5 * c + 2].material_id = 0;
2011  cells[5 * c + 3].material_id = 0;
2012  cells[5 * c + 4].material_id = 0;
2013  }
2014 
2015  tria.create_triangulation(vertices, cells, SubCellData());
2016 
2017  tria.reset_all_manifolds();
2018  tria.set_all_manifold_ids(0);
2019 
2020  for (auto &cell : tria.cell_iterators())
2021  {
2022  // identify faces on torus surface and set manifold to 1
2023  for (unsigned int f : GeometryInfo<3>::face_indices())
2024  {
2025  // faces 4 and 5 are those with normal vector aligned with torus
2026  // centerline
2027  if (cell->face(f)->at_boundary() && f != 4 && f != 5)
2028  {
2029  cell->face(f)->set_all_manifold_ids(1);
2030  }
2031  }
2032 
2033  // set manifold id to 2 for those cells that are on the torus centerline
2034  if (cell->material_id() == 1)
2035  {
2036  cell->set_all_manifold_ids(2);
2037  // reset to 0
2038  cell->set_material_id(0);
2039  }
2040  }
2041 
2042  tria.set_manifold(1, TorusManifold<3>(R, r));
2043  tria.set_manifold(2,
2044  CylindricalManifold<3>(Tensor<1, 3>({0., 1., 0.}),
2045  Point<3>()));
2047  transfinite.initialize(tria);
2048  tria.set_manifold(0, transfinite);
2049  }
2050 
2051 
2052 
2053  template <int dim, int spacedim>
2054  void
2056  const std::vector<Point<spacedim>> &vertices,
2057  const bool colorize)
2058  {
2060  ExcMessage("Wrong number of vertices."));
2061 
2062  // First create a hyper_rectangle and then deform it.
2063  hyper_cube(tria, 0, 1, colorize);
2064 
2066  tria.begin_active();
2067  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
2068  cell->vertex(i) = vertices[i];
2069 
2070  // Check that the order of the vertices makes sense, i.e., the volume of the
2071  // cell is positive.
2072  Assert(GridTools::volume(tria) > 0.,
2073  ExcMessage(
2074  "The volume of the cell is not greater than zero. "
2075  "This could be due to the wrong ordering of the vertices."));
2076  }
2077 
2078 
2079 
2080  template <>
2082  const Point<3> (&/*corners*/)[3],
2083  const bool /*colorize*/)
2084  {
2085  Assert(false, ExcNotImplemented());
2086  }
2087 
2088  template <>
2090  const Point<1> (&/*corners*/)[1],
2091  const bool /*colorize*/)
2092  {
2093  Assert(false, ExcNotImplemented());
2094  }
2095 
2096  // Implementation for 2D only
2097  template <>
2098  void parallelogram(Triangulation<2> &tria,
2099  const Point<2> (&corners)[2],
2100  const bool colorize)
2101  {
2102  Point<2> origin;
2103  std::array<Tensor<1, 2>, 2> edges;
2104  edges[0] = corners[0];
2105  edges[1] = corners[1];
2106  std::vector<unsigned int> subdivisions;
2107  subdivided_parallelepiped<2, 2>(
2108  tria, origin, edges, subdivisions, colorize);
2109  }
2110 
2111 
2112 
2113  template <int dim>
2114  void
2116  const Point<dim> (&corners)[dim],
2117  const bool colorize)
2118  {
2119  unsigned int n_subdivisions[dim];
2120  for (unsigned int i = 0; i < dim; ++i)
2121  n_subdivisions[i] = 1;
2122 
2123  // and call the function below
2124  subdivided_parallelepiped(tria, n_subdivisions, corners, colorize);
2125  }
2126 
2127  template <int dim>
2128  void
2130  const unsigned int n_subdivisions,
2131  const Point<dim> (&corners)[dim],
2132  const bool colorize)
2133  {
2134  // Equalize number of subdivisions in each dim-direction, their
2135  // validity will be checked later
2136  unsigned int n_subdivisions_[dim];
2137  for (unsigned int i = 0; i < dim; ++i)
2138  n_subdivisions_[i] = n_subdivisions;
2139 
2140  // and call the function below
2141  subdivided_parallelepiped(tria, n_subdivisions_, corners, colorize);
2142  }
2143 
2144  template <int dim>
2145  void
2147 #ifndef _MSC_VER
2148  const unsigned int (&n_subdivisions)[dim],
2149 #else
2150  const unsigned int *n_subdivisions,
2151 #endif
2152  const Point<dim> (&corners)[dim],
2153  const bool colorize)
2154  {
2155  Point<dim> origin;
2156  std::vector<unsigned int> subdivisions;
2157  std::array<Tensor<1, dim>, dim> edges;
2158  for (unsigned int i = 0; i < dim; ++i)
2159  {
2160  subdivisions.push_back(n_subdivisions[i]);
2161  edges[i] = corners[i];
2162  }
2163 
2164  subdivided_parallelepiped<dim, dim>(
2165  tria, origin, edges, subdivisions, colorize);
2166  }
2167 
2168  // Parallelepiped implementation in 1d, 2d, and 3d. @note The
2169  // implementation in 1d is similar to hyper_rectangle(), and in 2d is
2170  // similar to parallelogram().
2171  //
2172  // The GridReordering::reorder_grid is made use of towards the end of
2173  // this function. Thus the triangulation is explicitly constructed for
2174  // all dim here since it is slightly different in that respect
2175  // (cf. hyper_rectangle(), parallelogram()).
2176  template <int dim, int spacedim>
2177  void
2179  const Point<spacedim> & origin,
2180  const std::array<Tensor<1, spacedim>, dim> &edges,
2181  const std::vector<unsigned int> &subdivisions,
2182  const bool colorize)
2183  {
2184  std::vector<unsigned int> compute_subdivisions = subdivisions;
2185  if (compute_subdivisions.size() == 0)
2186  {
2187  compute_subdivisions.resize(dim, 1);
2188  }
2189 
2190  Assert(compute_subdivisions.size() == dim,
2191  ExcMessage("One subdivision must be provided for each dimension."));
2192  // check subdivisions
2193  for (unsigned int i = 0; i < dim; ++i)
2194  {
2195  Assert(compute_subdivisions[i] > 0,
2196  ExcInvalidRepetitions(subdivisions[i]));
2197  Assert(
2198  edges[i].norm() > 0,
2199  ExcMessage(
2200  "Edges in subdivided_parallelepiped() must not be degenerate."));
2201  }
2202 
2203  /*
2204  * Verify that the edge points to the right in 1D, vectors are oriented in
2205  * a counter clockwise direction in 2D, or form a right handed system in
2206  * 3D.
2207  */
2208  bool twisted_data = false;
2209  switch (dim)
2210  {
2211  case 1:
2212  {
2213  twisted_data = (edges[0][0] < 0);
2214  break;
2215  }
2216  case 2:
2217  {
2218  if (spacedim == 2) // this check does not make sense otherwise
2219  {
2220  const double plane_normal =
2221  edges[0][0] * edges[1][1] - edges[0][1] * edges[1][0];
2222  twisted_data = (plane_normal < 0.0);
2223  }
2224  break;
2225  }
2226  case 3:
2227  {
2228  // Check that the first two vectors are not linear combinations to
2229  // avoid zero division later on.
2230  Assert(std::abs(edges[0] * edges[1] /
2231  (edges[0].norm() * edges[1].norm()) -
2232  1.0) > 1.0e-15,
2233  ExcMessage(
2234  "Edges in subdivided_parallelepiped() must point in"
2235  " different directions."));
2236  const Tensor<1, spacedim> plane_normal =
2237  cross_product_3d(edges[0], edges[1]);
2238 
2239  /*
2240  * Ensure that edges 1, 2, and 3 form a right-handed set of
2241  * vectors. This works by applying the definition of the dot product
2242  *
2243  * cos(theta) = dot(x, y)/(norm(x)*norm(y))
2244  *
2245  * and then, since the normal vector and third edge should both
2246  * point away from the plane formed by the first two edges, the
2247  * angle between them must be between 0 and pi/2; hence we just need
2248  *
2249  * 0 < dot(x, y).
2250  */
2251  twisted_data = (plane_normal * edges[2] < 0.0);
2252  break;
2253  }
2254  default:
2255  Assert(false, ExcInternalError());
2256  }
2257  (void)twisted_data; // make the static analyzer happy
2258  Assert(
2259  !twisted_data,
2261  "The triangulation you are trying to create will consist of cells"
2262  " with negative measures. This is usually the result of input data"
2263  " that does not define a right-handed coordinate system. The usual"
2264  " fix for this is to ensure that in 1D the given point is to the"
2265  " right of the origin (or the given edge tensor is positive), in 2D"
2266  " that the two edges (and their cross product) obey the right-hand"
2267  " rule (which may usually be done by switching the order of the"
2268  " points or edge tensors), or in 3D that the edges form a"
2269  " right-handed coordinate system (which may also be accomplished by"
2270  " switching the order of the first two points or edge tensors)."));
2271 
2272  // Check corners do not overlap (unique)
2273  for (unsigned int i = 0; i < dim; ++i)
2274  for (unsigned int j = i + 1; j < dim; ++j)
2275  Assert((edges[i] != edges[j]),
2276  ExcMessage(
2277  "Degenerate edges of subdivided_parallelepiped encountered."));
2278 
2279  // Create a list of points
2280  std::vector<Point<spacedim>> points;
2281 
2282  switch (dim)
2283  {
2284  case 1:
2285  for (unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2286  points.push_back(origin + edges[0] / compute_subdivisions[0] * x);
2287  break;
2288 
2289  case 2:
2290  for (unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2291  for (unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2292  points.push_back(origin + edges[0] / compute_subdivisions[0] * x +
2293  edges[1] / compute_subdivisions[1] * y);
2294  break;
2295 
2296  case 3:
2297  for (unsigned int z = 0; z <= compute_subdivisions[2]; ++z)
2298  for (unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2299  for (unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2300  points.push_back(origin +
2301  edges[0] / compute_subdivisions[0] * x +
2302  edges[1] / compute_subdivisions[1] * y +
2303  edges[2] / compute_subdivisions[2] * z);
2304  break;
2305 
2306  default:
2307  Assert(false, ExcNotImplemented());
2308  }
2309 
2310  // Prepare cell data
2311  unsigned int n_cells = 1;
2312  for (unsigned int i = 0; i < dim; ++i)
2313  n_cells *= compute_subdivisions[i];
2314  std::vector<CellData<dim>> cells(n_cells);
2315 
2316  // Create fixed ordering of
2317  switch (dim)
2318  {
2319  case 1:
2320  for (unsigned int x = 0; x < compute_subdivisions[0]; ++x)
2321  {
2322  cells[x].vertices[0] = x;
2323  cells[x].vertices[1] = x + 1;
2324 
2325  // wipe material id
2326  cells[x].material_id = 0;
2327  }
2328  break;
2329 
2330  case 2:
2331  {
2332  // Shorthand
2333  const unsigned int n_dy = compute_subdivisions[1];
2334  const unsigned int n_dx = compute_subdivisions[0];
2335 
2336  for (unsigned int y = 0; y < n_dy; ++y)
2337  for (unsigned int x = 0; x < n_dx; ++x)
2338  {
2339  const unsigned int c = y * n_dx + x;
2340  cells[c].vertices[0] = y * (n_dx + 1) + x;
2341  cells[c].vertices[1] = y * (n_dx + 1) + x + 1;
2342  cells[c].vertices[2] = (y + 1) * (n_dx + 1) + x;
2343  cells[c].vertices[3] = (y + 1) * (n_dx + 1) + x + 1;
2344 
2345  // wipe material id
2346  cells[c].material_id = 0;
2347  }
2348  }
2349  break;
2350 
2351  case 3:
2352  {
2353  // Shorthand
2354  const unsigned int n_dz = compute_subdivisions[2];
2355  const unsigned int n_dy = compute_subdivisions[1];
2356  const unsigned int n_dx = compute_subdivisions[0];
2357 
2358  for (unsigned int z = 0; z < n_dz; ++z)
2359  for (unsigned int y = 0; y < n_dy; ++y)
2360  for (unsigned int x = 0; x < n_dx; ++x)
2361  {
2362  const unsigned int c = z * n_dy * n_dx + y * n_dx + x;
2363 
2364  cells[c].vertices[0] =
2365  z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2366  cells[c].vertices[1] =
2367  z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x + 1;
2368  cells[c].vertices[2] =
2369  z * (n_dy + 1) * (n_dx + 1) + (y + 1) * (n_dx + 1) + x;
2370  cells[c].vertices[3] = z * (n_dy + 1) * (n_dx + 1) +
2371  (y + 1) * (n_dx + 1) + x + 1;
2372  cells[c].vertices[4] =
2373  (z + 1) * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2374  cells[c].vertices[5] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2375  y * (n_dx + 1) + x + 1;
2376  cells[c].vertices[6] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2377  (y + 1) * (n_dx + 1) + x;
2378  cells[c].vertices[7] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2379  (y + 1) * (n_dx + 1) + x + 1;
2380 
2381  // wipe material id
2382  cells[c].material_id = 0;
2383  }
2384  break;
2385  }
2386 
2387  default:
2388  Assert(false, ExcNotImplemented());
2389  }
2390 
2391  // Create triangulation
2392  // reorder the cells to ensure that they satisfy the convention for
2393  // edge and face directions
2395  tria.create_triangulation(points, cells, SubCellData());
2396 
2397  // Finally assign boundary indicators according to hyper_rectangle
2398  if (colorize)
2399  {
2401  tria.begin_active(),
2402  endc = tria.end();
2403  for (; cell != endc; ++cell)
2404  {
2405  for (unsigned int face = 0;
2406  face < GeometryInfo<dim>::faces_per_cell;
2407  ++face)
2408  {
2409  if (cell->face(face)->at_boundary())
2410  cell->face(face)->set_boundary_id(face);
2411  }
2412  }
2413  }
2414  }
2415 
2416 
2417  template <int dim, int spacedim>
2418  void
2420  const unsigned int repetitions,
2421  const double left,
2422  const double right,
2423  const bool colorize)
2424  {
2425  Assert(repetitions >= 1, ExcInvalidRepetitions(repetitions));
2426  Assert(left < right,
2427  ExcMessage("Invalid left-to-right bounds of hypercube"));
2428 
2429  Point<dim> p0, p1;
2430  for (unsigned int i = 0; i < dim; ++i)
2431  {
2432  p0[i] = left;
2433  p1[i] = right;
2434  }
2435 
2436  std::vector<unsigned int> reps(dim, repetitions);
2437  subdivided_hyper_rectangle(tria, reps, p0, p1, colorize);
2438  }
2439 
2440 
2441 
2442  template <int dim, int spacedim>
2443  void
2445  const std::vector<unsigned int> &repetitions,
2446  const Point<dim> & p_1,
2447  const Point<dim> & p_2,
2448  const bool colorize)
2449  {
2450  Assert(repetitions.size() == dim, ExcInvalidRepetitionsDimension(dim));
2451 
2452  // First, extend dimensions from dim to spacedim and
2453  // normalize such that p1 is lower in all coordinate
2454  // directions. Additional entries will be 0.
2455  Point<spacedim> p1, p2;
2456  for (unsigned int i = 0; i < dim; ++i)
2457  {
2458  p1(i) = std::min(p_1(i), p_2(i));
2459  p2(i) = std::max(p_1(i), p_2(i));
2460  }
2461 
2462  // calculate deltas and validate input
2463  std::vector<Point<spacedim>> delta(dim);
2464  for (unsigned int i = 0; i < dim; ++i)
2465  {
2466  Assert(repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
2467 
2468  delta[i][i] = (p2[i] - p1[i]) / repetitions[i];
2469  Assert(
2470  delta[i][i] > 0.0,
2471  ExcMessage(
2472  "The first dim entries of coordinates of p1 and p2 need to be different."));
2473  }
2474 
2475  // then generate the points
2476  std::vector<Point<spacedim>> points;
2477  switch (dim)
2478  {
2479  case 1:
2480  for (unsigned int x = 0; x <= repetitions[0]; ++x)
2481  points.push_back(p1 + x * delta[0]);
2482  break;
2483 
2484  case 2:
2485  for (unsigned int y = 0; y <= repetitions[1]; ++y)
2486  for (unsigned int x = 0; x <= repetitions[0]; ++x)
2487  points.push_back(p1 + x * delta[0] + y * delta[1]);
2488  break;
2489 
2490  case 3:
2491  for (unsigned int z = 0; z <= repetitions[2]; ++z)
2492  for (unsigned int y = 0; y <= repetitions[1]; ++y)
2493  for (unsigned int x = 0; x <= repetitions[0]; ++x)
2494  points.push_back(p1 + x * delta[0] + y * delta[1] +
2495  z * delta[2]);
2496  break;
2497 
2498  default:
2499  Assert(false, ExcNotImplemented());
2500  }
2501 
2502  // next create the cells
2503  std::vector<CellData<dim>> cells;
2504  switch (dim)
2505  {
2506  case 1:
2507  {
2508  cells.resize(repetitions[0]);
2509  for (unsigned int x = 0; x < repetitions[0]; ++x)
2510  {
2511  cells[x].vertices[0] = x;
2512  cells[x].vertices[1] = x + 1;
2513  cells[x].material_id = 0;
2514  }
2515  break;
2516  }
2517 
2518  case 2:
2519  {
2520  cells.resize(repetitions[1] * repetitions[0]);
2521  for (unsigned int y = 0; y < repetitions[1]; ++y)
2522  for (unsigned int x = 0; x < repetitions[0]; ++x)
2523  {
2524  const unsigned int c = x + y * repetitions[0];
2525  cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
2526  cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2527  cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2528  cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
2529  cells[c].material_id = 0;
2530  }
2531  break;
2532  }
2533 
2534  case 3:
2535  {
2536  const unsigned int n_x = (repetitions[0] + 1);
2537  const unsigned int n_xy =
2538  (repetitions[0] + 1) * (repetitions[1] + 1);
2539 
2540  cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
2541  for (unsigned int z = 0; z < repetitions[2]; ++z)
2542  for (unsigned int y = 0; y < repetitions[1]; ++y)
2543  for (unsigned int x = 0; x < repetitions[0]; ++x)
2544  {
2545  const unsigned int c = x + y * repetitions[0] +
2546  z * repetitions[0] * repetitions[1];
2547  cells[c].vertices[0] = z * n_xy + y * n_x + x;
2548  cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2549  cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2550  cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2551  cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2552  cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2553  cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2554  cells[c].vertices[7] =
2555  (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2556  cells[c].material_id = 0;
2557  }
2558  break;
2559  }
2560 
2561  default:
2562  Assert(false, ExcNotImplemented());
2563  }
2564 
2565  tria.create_triangulation(points, cells, SubCellData());
2566 
2567  if (colorize)
2568  {
2569  // to colorize, run through all
2570  // faces of all cells and set
2571  // boundary indicator to the
2572  // correct value if it was 0.
2573 
2574  // use a large epsilon to
2575  // compare numbers to avoid
2576  // roundoff problems.
2577  double epsilon = 10;
2578  for (unsigned int i = 0; i < dim; ++i)
2579  epsilon = std::min(epsilon, 0.01 * delta[i][i]);
2580  Assert(epsilon > 0,
2581  ExcMessage(
2582  "The distance between corner points must be positive."))
2583 
2584  // actual code is external since
2585  // 1-D is different from 2/3D.
2586  colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
2587  }
2588  }
2589 
2590 
2591 
2592  template <int dim>
2593  void
2595  const std::vector<std::vector<double>> &step_sz,
2596  const Point<dim> & p_1,
2597  const Point<dim> & p_2,
2598  const bool colorize)
2599  {
2600  Assert(step_sz.size() == dim, ExcInvalidRepetitionsDimension(dim));
2601 
2602  // First, normalize input such that
2603  // p1 is lower in all coordinate
2604  // directions and check the consistency of
2605  // step sizes, i.e. that they all
2606  // add up to the sizes specified by
2607  // p_1 and p_2
2608  Point<dim> p1(p_1);
2609  Point<dim> p2(p_2);
2610  std::vector<std::vector<double>> step_sizes(step_sz);
2611 
2612  for (unsigned int i = 0; i < dim; ++i)
2613  {
2614  if (p1(i) > p2(i))
2615  {
2616  std::swap(p1(i), p2(i));
2617  std::reverse(step_sizes[i].begin(), step_sizes[i].end());
2618  }
2619 
2620  double x = 0;
2621  for (unsigned int j = 0; j < step_sizes.at(i).size(); j++)
2622  x += step_sizes[i][j];
2623  Assert(std::fabs(x - (p2(i) - p1(i))) <= 1e-12 * std::fabs(x),
2624  ExcMessage(
2625  "The sequence of step sizes in coordinate direction " +
2627  " must be equal to the distance of the two given "
2628  "points in this coordinate direction."));
2629  }
2630 
2631 
2632  // then generate the necessary
2633  // points
2634  std::vector<Point<dim>> points;
2635  switch (dim)
2636  {
2637  case 1:
2638  {
2639  double x = 0;
2640  for (unsigned int i = 0;; ++i)
2641  {
2642  points.push_back(Point<dim>(p1[0] + x));
2643 
2644  // form partial sums. in
2645  // the last run through
2646  // avoid accessing
2647  // non-existent values
2648  // and exit early instead
2649  if (i == step_sizes[0].size())
2650  break;
2651 
2652  x += step_sizes[0][i];
2653  }
2654  break;
2655  }
2656 
2657  case 2:
2658  {
2659  double y = 0;
2660  for (unsigned int j = 0;; ++j)
2661  {
2662  double x = 0;
2663  for (unsigned int i = 0;; ++i)
2664  {
2665  points.push_back(Point<dim>(p1[0] + x, p1[1] + y));
2666  if (i == step_sizes[0].size())
2667  break;
2668 
2669  x += step_sizes[0][i];
2670  }
2671 
2672  if (j == step_sizes[1].size())
2673  break;
2674 
2675  y += step_sizes[1][j];
2676  }
2677  break;
2678  }
2679  case 3:
2680  {
2681  double z = 0;
2682  for (unsigned int k = 0;; ++k)
2683  {
2684  double y = 0;
2685  for (unsigned int j = 0;; ++j)
2686  {
2687  double x = 0;
2688  for (unsigned int i = 0;; ++i)
2689  {
2690  points.push_back(
2691  Point<dim>(p1[0] + x, p1[1] + y, p1[2] + z));
2692  if (i == step_sizes[0].size())
2693  break;
2694 
2695  x += step_sizes[0][i];
2696  }
2697 
2698  if (j == step_sizes[1].size())
2699  break;
2700 
2701  y += step_sizes[1][j];
2702  }
2703 
2704  if (k == step_sizes[2].size())
2705  break;
2706 
2707  z += step_sizes[2][k];
2708  }
2709  break;
2710  }
2711 
2712  default:
2713  Assert(false, ExcNotImplemented());
2714  }
2715 
2716  // next create the cells
2717  // Prepare cell data
2718  std::vector<CellData<dim>> cells;
2719  switch (dim)
2720  {
2721  case 1:
2722  {
2723  cells.resize(step_sizes[0].size());
2724  for (unsigned int x = 0; x < step_sizes[0].size(); ++x)
2725  {
2726  cells[x].vertices[0] = x;
2727  cells[x].vertices[1] = x + 1;
2728  cells[x].material_id = 0;
2729  }
2730  break;
2731  }
2732 
2733  case 2:
2734  {
2735  cells.resize(step_sizes[1].size() * step_sizes[0].size());
2736  for (unsigned int y = 0; y < step_sizes[1].size(); ++y)
2737  for (unsigned int x = 0; x < step_sizes[0].size(); ++x)
2738  {
2739  const unsigned int c = x + y * step_sizes[0].size();
2740  cells[c].vertices[0] = y * (step_sizes[0].size() + 1) + x;
2741  cells[c].vertices[1] = y * (step_sizes[0].size() + 1) + x + 1;
2742  cells[c].vertices[2] =
2743  (y + 1) * (step_sizes[0].size() + 1) + x;
2744  cells[c].vertices[3] =
2745  (y + 1) * (step_sizes[0].size() + 1) + x + 1;
2746  cells[c].material_id = 0;
2747  }
2748  break;
2749  }
2750 
2751  case 3:
2752  {
2753  const unsigned int n_x = (step_sizes[0].size() + 1);
2754  const unsigned int n_xy =
2755  (step_sizes[0].size() + 1) * (step_sizes[1].size() + 1);
2756 
2757  cells.resize(step_sizes[2].size() * step_sizes[1].size() *
2758  step_sizes[0].size());
2759  for (unsigned int z = 0; z < step_sizes[2].size(); ++z)
2760  for (unsigned int y = 0; y < step_sizes[1].size(); ++y)
2761  for (unsigned int x = 0; x < step_sizes[0].size(); ++x)
2762  {
2763  const unsigned int c =
2764  x + y * step_sizes[0].size() +
2765  z * step_sizes[0].size() * step_sizes[1].size();
2766  cells[c].vertices[0] = z * n_xy + y * n_x + x;
2767  cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2768  cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2769  cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2770  cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2771  cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2772  cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2773  cells[c].vertices[7] =
2774  (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2775  cells[c].material_id = 0;
2776  }
2777  break;
2778  }
2779 
2780  default:
2781  Assert(false, ExcNotImplemented());
2782  }
2783 
2784  tria.create_triangulation(points, cells, SubCellData());
2785 
2786  if (colorize)
2787  {
2788  // to colorize, run through all
2789  // faces of all cells and set
2790  // boundary indicator to the
2791  // correct value if it was 0.
2792 
2793  // use a large epsilon to
2794  // compare numbers to avoid
2795  // roundoff problems.
2796  double min_size =
2797  *std::min_element(step_sizes[0].begin(), step_sizes[0].end());
2798  for (unsigned int i = 1; i < dim; ++i)
2799  min_size = std::min(min_size,
2800  *std::min_element(step_sizes[i].begin(),
2801  step_sizes[i].end()));
2802  const double epsilon = 0.01 * min_size;
2803 
2804  // actual code is external since
2805  // 1-D is different from 2/3D.
2806  colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
2807  }
2808  }
2809 
2810 
2811 
2812  template <>
2813  void
2815  const std::vector<std::vector<double>> &spacing,
2816  const Point<1> & p,
2817  const Table<1, types::material_id> &material_id,
2818  const bool colorize)
2819  {
2820  Assert(spacing.size() == 1, ExcInvalidRepetitionsDimension(1));
2821 
2822  const unsigned int n_cells = material_id.size(0);
2823 
2824  Assert(spacing[0].size() == n_cells, ExcInvalidRepetitionsDimension(1));
2825 
2826  double delta = std::numeric_limits<double>::max();
2827  for (unsigned int i = 0; i < n_cells; i++)
2828  {
2829  Assert(spacing[0][i] >= 0, ExcInvalidRepetitions(-1));
2830  delta = std::min(delta, spacing[0][i]);
2831  }
2832 
2833  // generate the necessary points
2834  std::vector<Point<1>> points;
2835  double ax = p[0];
2836  for (unsigned int x = 0; x <= n_cells; ++x)
2837  {
2838  points.emplace_back(ax);
2839  if (x < n_cells)
2840  ax += spacing[0][x];
2841  }
2842  // create the cells
2843  unsigned int n_val_cells = 0;
2844  for (unsigned int i = 0; i < n_cells; i++)
2845  if (material_id[i] != numbers::invalid_material_id)
2846  n_val_cells++;
2847 
2848  std::vector<CellData<1>> cells(n_val_cells);
2849  unsigned int id = 0;
2850  for (unsigned int x = 0; x < n_cells; ++x)
2851  if (material_id[x] != numbers::invalid_material_id)
2852  {
2853  cells[id].vertices[0] = x;
2854  cells[id].vertices[1] = x + 1;
2855  cells[id].material_id = material_id[x];
2856  id++;
2857  }
2858  // create triangulation
2859  SubCellData t;
2860  GridTools::delete_unused_vertices(points, cells, t);
2861 
2862  tria.create_triangulation(points, cells, t);
2863 
2864  // set boundary indicator
2865  if (colorize)
2866  Assert(false, ExcNotImplemented());
2867  }
2868 
2869 
2870  template <>
2871  void
2873  const std::vector<std::vector<double>> &spacing,
2874  const Point<2> & p,
2875  const Table<2, types::material_id> &material_id,
2876  const bool colorize)
2877  {
2878  Assert(spacing.size() == 2, ExcInvalidRepetitionsDimension(2));
2879 
2880  std::vector<unsigned int> repetitions(2);
2881  unsigned int n_cells = 1;
2882  double delta = std::numeric_limits<double>::max();
2883  for (unsigned int i = 0; i < 2; i++)
2884  {
2885  repetitions[i] = spacing[i].size();
2886  n_cells *= repetitions[i];
2887  for (unsigned int j = 0; j < repetitions[i]; j++)
2888  {
2889  Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
2890  delta = std::min(delta, spacing[i][j]);
2891  }
2892  Assert(material_id.size(i) == repetitions[i],
2894  }
2895 
2896  // generate the necessary points
2897  std::vector<Point<2>> points;
2898  double ay = p[1];
2899  for (unsigned int y = 0; y <= repetitions[1]; ++y)
2900  {
2901  double ax = p[0];
2902  for (unsigned int x = 0; x <= repetitions[0]; ++x)
2903  {
2904  points.emplace_back(ax, ay);
2905  if (x < repetitions[0])
2906  ax += spacing[0][x];
2907  }
2908  if (y < repetitions[1])
2909  ay += spacing[1][y];
2910  }
2911 
2912  // create the cells
2913  unsigned int n_val_cells = 0;
2914  for (unsigned int i = 0; i < material_id.size(0); i++)
2915  for (unsigned int j = 0; j < material_id.size(1); j++)
2916  if (material_id[i][j] != numbers::invalid_material_id)
2917  n_val_cells++;
2918 
2919  std::vector<CellData<2>> cells(n_val_cells);
2920  unsigned int id = 0;
2921  for (unsigned int y = 0; y < repetitions[1]; ++y)
2922  for (unsigned int x = 0; x < repetitions[0]; ++x)
2923  if (material_id[x][y] != numbers::invalid_material_id)
2924  {
2925  cells[id].vertices[0] = y * (repetitions[0] + 1) + x;
2926  cells[id].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2927  cells[id].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2928  cells[id].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
2929  cells[id].material_id = material_id[x][y];
2930  id++;
2931  }
2932 
2933  // create triangulation
2934  SubCellData t;
2935  GridTools::delete_unused_vertices(points, cells, t);
2936 
2937  tria.create_triangulation(points, cells, t);
2938 
2939  // set boundary indicator
2940  if (colorize)
2941  {
2942  double eps = 0.01 * delta;
2943  Triangulation<2>::cell_iterator cell = tria.begin(), endc = tria.end();
2944  for (; cell != endc; ++cell)
2945  {
2946  Point<2> cell_center = cell->center();
2947  for (unsigned int f : GeometryInfo<2>::face_indices())
2948  if (cell->face(f)->boundary_id() == 0)
2949  {
2950  Point<2> face_center = cell->face(f)->center();
2951  for (unsigned int i = 0; i < 2; ++i)
2952  {
2953  if (face_center[i] < cell_center[i] - eps)
2954  cell->face(f)->set_boundary_id(i * 2);
2955  if (face_center[i] > cell_center[i] + eps)
2956  cell->face(f)->set_boundary_id(i * 2 + 1);
2957  }
2958  }
2959  }
2960  }
2961  }
2962 
2963 
2964  template <>
2965  void
2967  const std::vector<std::vector<double>> &spacing,
2968  const Point<3> & p,
2969  const Table<3, types::material_id> &material_id,
2970  const bool colorize)
2971  {
2972  const unsigned int dim = 3;
2973 
2974  Assert(spacing.size() == dim, ExcInvalidRepetitionsDimension(dim));
2975 
2976  std::vector<unsigned int> repetitions(dim);
2977  unsigned int n_cells = 1;
2978  double delta = std::numeric_limits<double>::max();
2979  for (unsigned int i = 0; i < dim; i++)
2980  {
2981  repetitions[i] = spacing[i].size();
2982  n_cells *= repetitions[i];
2983  for (unsigned int j = 0; j < repetitions[i]; j++)
2984  {
2985  Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
2986  delta = std::min(delta, spacing[i][j]);
2987  }
2988  Assert(material_id.size(i) == repetitions[i],
2990  }
2991 
2992  // generate the necessary points
2993  std::vector<Point<dim>> points;
2994  double az = p[2];
2995  for (unsigned int z = 0; z <= repetitions[2]; ++z)
2996  {
2997  double ay = p[1];
2998  for (unsigned int y = 0; y <= repetitions[1]; ++y)
2999  {
3000  double ax = p[0];
3001  for (unsigned int x = 0; x <= repetitions[0]; ++x)
3002  {
3003  points.emplace_back(ax, ay, az);
3004  if (x < repetitions[0])
3005  ax += spacing[0][x];
3006  }
3007  if (y < repetitions[1])
3008  ay += spacing[1][y];
3009  }
3010  if (z < repetitions[2])
3011  az += spacing[2][z];
3012  }
3013 
3014  // create the cells
3015  unsigned int n_val_cells = 0;
3016  for (unsigned int i = 0; i < material_id.size(0); i++)
3017  for (unsigned int j = 0; j < material_id.size(1); j++)
3018  for (unsigned int k = 0; k < material_id.size(2); k++)
3019  if (material_id[i][j][k] != numbers::invalid_material_id)
3020  n_val_cells++;
3021 
3022  std::vector<CellData<dim>> cells(n_val_cells);
3023  unsigned int id = 0;
3024  const unsigned int n_x = (repetitions[0] + 1);
3025  const unsigned int n_xy = (repetitions[0] + 1) * (repetitions[1] + 1);
3026  for (unsigned int z = 0; z < repetitions[2]; ++z)
3027  for (unsigned int y = 0; y < repetitions[1]; ++y)
3028  for (unsigned int x = 0; x < repetitions[0]; ++x)
3029  if (material_id[x][y][z] != numbers::invalid_material_id)
3030  {
3031  cells[id].vertices[0] = z * n_xy + y * n_x + x;
3032  cells[id].vertices[1] = z * n_xy + y * n_x + x + 1;
3033  cells[id].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3034  cells[id].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3035  cells[id].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3036  cells[id].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3037  cells[id].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3038  cells[id].vertices[7] = (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3039  cells[id].material_id = material_id[x][y][z];
3040  id++;
3041  }
3042 
3043  // create triangulation
3044  SubCellData t;
3045  GridTools::delete_unused_vertices(points, cells, t);
3046 
3047  tria.create_triangulation(points, cells, t);
3048 
3049  // set boundary indicator
3050  if (colorize)
3051  {
3052  double eps = 0.01 * delta;
3054  endc = tria.end();
3055  for (; cell != endc; ++cell)
3056  {
3057  Point<dim> cell_center = cell->center();
3058  for (auto f : GeometryInfo<dim>::face_indices())
3059  if (cell->face(f)->boundary_id() == 0)
3060  {
3061  Point<dim> face_center = cell->face(f)->center();
3062  for (unsigned int i = 0; i < dim; ++i)
3063  {
3064  if (face_center[i] < cell_center[i] - eps)
3065  cell->face(f)->set_boundary_id(i * 2);
3066  if (face_center[i] > cell_center[i] + eps)
3067  cell->face(f)->set_boundary_id(i * 2 + 1);
3068  }
3069  }
3070  }
3071  }
3072  }
3073 
3074  template <int dim, int spacedim>
3075  void
3077  const std::vector<unsigned int> &holes)
3078  {
3079  AssertDimension(holes.size(), dim);
3080  // The corner points of the first cell. If there is a desire at
3081  // some point to change the geometry of the cells, they can be
3082  // made an argument to the function.
3083 
3084  Point<spacedim> p1;
3085  Point<spacedim> p2;
3086  for (unsigned int d = 0; d < dim; ++d)
3087  p2(d) = 1.;
3088 
3089  // then check that all repetitions
3090  // are >= 1, and calculate deltas
3091  // convert repetitions from double
3092  // to int by taking the ceiling.
3093  std::vector<Point<spacedim>> delta(dim);
3094  unsigned int repetitions[dim];
3095  for (unsigned int i = 0; i < dim; ++i)
3096  {
3097  Assert(holes[i] >= 1,
3098  ExcMessage("At least one hole needed in each direction"));
3099  repetitions[i] = 2 * holes[i] + 1;
3100  delta[i][i] = (p2[i] - p1[i]);
3101  }
3102 
3103  // then generate the necessary
3104  // points
3105  std::vector<Point<spacedim>> points;
3106  switch (dim)
3107  {
3108  case 1:
3109  for (unsigned int x = 0; x <= repetitions[0]; ++x)
3110  points.push_back(p1 + x * delta[0]);
3111  break;
3112 
3113  case 2:
3114  for (unsigned int y = 0; y <= repetitions[1]; ++y)
3115  for (unsigned int x = 0; x <= repetitions[0]; ++x)
3116  points.push_back(p1 + x * delta[0] + y * delta[1]);
3117  break;
3118 
3119  case 3:
3120  for (unsigned int z = 0; z <= repetitions[2]; ++z)
3121  for (unsigned int y = 0; y <= repetitions[1]; ++y)
3122  for (unsigned int x = 0; x <= repetitions[0]; ++x)
3123  points.push_back(p1 + x * delta[0] + y * delta[1] +
3124  z * delta[2]);
3125  break;
3126 
3127  default:
3128  Assert(false, ExcNotImplemented());
3129  }
3130 
3131  // next create the cells
3132  // Prepare cell data
3133  std::vector<CellData<dim>> cells;
3134  switch (dim)
3135  {
3136  case 2:
3137  {
3138  cells.resize(repetitions[1] * repetitions[0] - holes[1] * holes[0]);
3139  unsigned int c = 0;
3140  for (unsigned int y = 0; y < repetitions[1]; ++y)
3141  for (unsigned int x = 0; x < repetitions[0]; ++x)
3142  {
3143  if ((x % 2 == 1) && (y % 2 == 1))
3144  continue;
3145  Assert(c < cells.size(), ExcInternalError());
3146  cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
3147  cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
3148  cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
3149  cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
3150  cells[c].material_id = 0;
3151  ++c;
3152  }
3153  break;
3154  }
3155 
3156  case 3:
3157  {
3158  const unsigned int n_x = (repetitions[0] + 1);
3159  const unsigned int n_xy =
3160  (repetitions[0] + 1) * (repetitions[1] + 1);
3161 
3162  cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
3163 
3164  unsigned int c = 0;
3165  for (unsigned int z = 0; z < repetitions[2]; ++z)
3166  for (unsigned int y = 0; y < repetitions[1]; ++y)
3167  for (unsigned int x = 0; x < repetitions[0]; ++x)
3168  {
3169  Assert(c < cells.size(), ExcInternalError());
3170  cells[c].vertices[0] = z * n_xy + y * n_x + x;
3171  cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
3172  cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3173  cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3174  cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3175  cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3176  cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3177  cells[c].vertices[7] =
3178  (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3179  cells[c].material_id = 0;
3180  ++c;
3181  }
3182  break;
3183  }
3184 
3185  default:
3186  Assert(false, ExcNotImplemented());
3187  }
3188 
3189  tria.create_triangulation(points, cells, SubCellData());
3190  }
3191 
3192 
3193 
3194  template <>
3195  void plate_with_a_hole(Triangulation<1> & /*tria*/,
3196  const double /*inner_radius*/,
3197  const double /*outer_radius*/,
3198  const double /*pad_bottom*/,
3199  const double /*pad_top*/,
3200  const double /*pad_left*/,
3201  const double /*pad_right*/,
3202  const Point<1> /*center*/,
3203  const types::manifold_id /*polar_manifold_id*/,
3204  const types::manifold_id /*tfi_manifold_id*/,
3205  const double /*L*/,
3206  const unsigned int /*n_slices*/,
3207  const bool /*colorize*/)
3208  {
3209  Assert(false, ExcNotImplemented());
3210  }
3211 
3212 
3213 
3214  template <>
3215  void channel_with_cylinder(Triangulation<1> & /*tria*/,
3216  const double /*shell_region_width*/,
3217  const unsigned int /*n_shells*/,
3218  const double /*skewness*/,
3219  const bool /*colorize*/)
3220  {
3221  Assert(false, ExcNotImplemented());
3222  }
3223 
3224 
3225 
3226  namespace internal
3227  {
3228  // helper function to check if point is in 2D box
3229  bool inline point_in_2d_box(const Point<2> &p,
3230  const Point<2> &c,
3231  const double radius)
3232  {
3233  return (std::abs(p[0] - c[0]) < radius) &&
3234  (std::abs(p[1] - c[1]) < radius);
3235  }
3236 
3237 
3238 
3239  // Find the minimal distance between two vertices. This is useful for
3240  // computing a tolerance for merging vertices in
3241  // GridTools::merge_triangulations.
3242  template <int dim, int spacedim>
3243  double
3244  minimal_vertex_distance(const Triangulation<dim, spacedim> &triangulation)
3245  {
3246  double length = std::numeric_limits<double>::max();
3247  for (const auto &cell : triangulation.active_cell_iterators())
3248  for (unsigned int n = 0; n < GeometryInfo<dim>::lines_per_cell; ++n)
3249  length = std::min(length, cell->line(n)->diameter());
3250  return length;
3251  }
3252  } // namespace internal
3253 
3254 
3255 
3256  template <>
3257  void plate_with_a_hole(Triangulation<2> & tria,
3258  const double inner_radius,
3259  const double outer_radius,
3260  const double pad_bottom,
3261  const double pad_top,
3262  const double pad_left,
3263  const double pad_right,
3264  const Point<2> new_center,
3265  const types::manifold_id polar_manifold_id,
3266  const types::manifold_id tfi_manifold_id,
3267  const double L,
3268  const unsigned int /*n_slices*/,
3269  const bool colorize)
3270  {
3271  const bool with_padding =
3272  pad_bottom > 0 || pad_top > 0 || pad_left > 0 || pad_right > 0;
3273 
3274  Assert(pad_bottom >= 0., ExcMessage("Negative bottom padding."));
3275  Assert(pad_top >= 0., ExcMessage("Negative top padding."));
3276  Assert(pad_left >= 0., ExcMessage("Negative left padding."));
3277  Assert(pad_right >= 0., ExcMessage("Negative right padding."));
3278 
3279  const Point<2> center;
3280 
3281  auto min_line_length = [](const Triangulation<2> &tria) -> double {
3282  double length = std::numeric_limits<double>::max();
3283  for (const auto &cell : tria.active_cell_iterators())
3284  for (unsigned int n = 0; n < GeometryInfo<2>::lines_per_cell; ++n)
3285  length = std::min(length, cell->line(n)->diameter());
3286  return length;
3287  };
3288 
3289  // start by setting up the cylinder triangulation
3290  Triangulation<2> cylinder_tria_maybe;
3291  Triangulation<2> &cylinder_tria = with_padding ? cylinder_tria_maybe : tria;
3293  inner_radius,
3294  outer_radius,
3295  L,
3296  /*repetitions*/ 1,
3297  colorize);
3298 
3299  // we will deal with face manifold ids after we merge triangulations
3300  for (const auto &cell : cylinder_tria.active_cell_iterators())
3301  cell->set_manifold_id(tfi_manifold_id);
3302 
3303  const Point<2> bl(-outer_radius - pad_left, -outer_radius - pad_bottom);
3304  const Point<2> tr(outer_radius + pad_right, outer_radius + pad_top);
3305  if (with_padding)
3306  {
3307  // hyper_cube_with_cylindrical_hole will have 2 cells along
3308  // each face, so the element size is outer_radius
3309 
3310  auto add_sizes = [](std::vector<double> &step_sizes,
3311  const double padding,
3312  const double h) -> void {
3313  // use std::round instead of std::ceil to improve aspect ratio
3314  // in case padding is only slightly larger than h.
3315  const auto rounded =
3316  static_cast<unsigned int>(std::round(padding / h));
3317  // in case padding is much smaller than h, make sure we
3318  // have at least 1 element
3319  const unsigned int num = (padding > 0. && rounded == 0) ? 1 : rounded;
3320  for (unsigned int i = 0; i < num; ++i)
3321  step_sizes.push_back(padding / num);
3322  };
3323 
3324  std::vector<std::vector<double>> step_sizes(2);
3325  // x-coord
3326  // left:
3327  add_sizes(step_sizes[0], pad_left, outer_radius);
3328  // center
3329  step_sizes[0].push_back(outer_radius);
3330  step_sizes[0].push_back(outer_radius);
3331  // right
3332  add_sizes(step_sizes[0], pad_right, outer_radius);
3333  // y-coord
3334  // bottom
3335  add_sizes(step_sizes[1], pad_bottom, outer_radius);
3336  // center
3337  step_sizes[1].push_back(outer_radius);
3338  step_sizes[1].push_back(outer_radius);
3339  // top
3340  add_sizes(step_sizes[1], pad_top, outer_radius);
3341 
3342  // now create bulk
3343  Triangulation<2> bulk_tria;
3345  bulk_tria, step_sizes, bl, tr, colorize);
3346 
3347  // now remove cells reserved from the cylindrical hole
3348  std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3349  for (const auto &cell : bulk_tria.active_cell_iterators())
3350  if (internal::point_in_2d_box(cell->center(), center, outer_radius))
3351  cells_to_remove.insert(cell);
3352 
3353  Triangulation<2> tria_without_cylinder;
3355  bulk_tria, cells_to_remove, tria_without_cylinder);
3356 
3357  const double tolerance =
3358  std::min(min_line_length(tria_without_cylinder),
3359  min_line_length(cylinder_tria)) /
3360  2.0;
3361 
3362  GridGenerator::merge_triangulations(tria_without_cylinder,
3363  cylinder_tria,
3364  tria,
3365  tolerance);
3366  }
3367 
3368  // now set manifold ids:
3369  for (const auto &cell : tria.active_cell_iterators())
3370  {
3371  // set all non-boundary manifold ids on the cells that came from the
3372  // grid around the cylinder to the new TFI manifold id.
3373  if (cell->manifold_id() == tfi_manifold_id)
3374  {
3375  for (unsigned int face_n = 0;
3376  face_n < GeometryInfo<2>::faces_per_cell;
3377  ++face_n)
3378  {
3379  const auto &face = cell->face(face_n);
3380  if (face->at_boundary() &&
3381  internal::point_in_2d_box(face->center(),
3382  center,
3383  outer_radius))
3384  face->set_manifold_id(polar_manifold_id);
3385  else
3386  face->set_manifold_id(tfi_manifold_id);
3387  }
3388  }
3389  else
3390  {
3391  // ensure that all other manifold ids (including the faces
3392  // opposite the cylinder) are set to the flat id
3393  cell->set_all_manifold_ids(numbers::flat_manifold_id);
3394  }
3395  }
3396 
3397  static constexpr double tol =
3398  std::numeric_limits<double>::epsilon() * 10000;
3399  if (colorize)
3400  for (const auto &cell : tria.active_cell_iterators())
3401  for (unsigned int face_n = 0; face_n < GeometryInfo<2>::faces_per_cell;
3402  ++face_n)
3403  {
3404  const auto face = cell->face(face_n);
3405  if (face->at_boundary())
3406  {
3407  const Point<2> center = face->center();
3408  // left side
3409  if (std::abs(center[0] - bl[0]) < tol * std::abs(bl[0]))
3410  face->set_boundary_id(0);
3411  // right side
3412  else if (std::abs(center[0] - tr[0]) < tol * std::abs(tr[0]))
3413  face->set_boundary_id(1);
3414  // bottom
3415  else if (std::abs(center[1] - bl[1]) < tol * std::abs(bl[1]))
3416  face->set_boundary_id(2);
3417  // top
3418  else if (std::abs(center[1] - tr[1]) < tol * std::abs(tr[1]))
3419  face->set_boundary_id(3);
3420  // cylinder boundary
3421  else
3422  {
3423  Assert(cell->manifold_id() == tfi_manifold_id,
3424  ExcInternalError());
3425  face->set_boundary_id(4);
3426  }
3427  }
3428  }
3429 
3430  // move to the new center
3431  GridTools::shift(new_center, tria);
3432 
3433  PolarManifold<2> polar_manifold(new_center);
3434  tria.set_manifold(polar_manifold_id, polar_manifold);
3435  TransfiniteInterpolationManifold<2> inner_manifold;
3436  inner_manifold.initialize(tria);
3437  tria.set_manifold(tfi_manifold_id, inner_manifold);
3438  }
3439 
3440 
3441 
3442  template <>
3443  void plate_with_a_hole(Triangulation<3> & tria,
3444  const double inner_radius,
3445  const double outer_radius,
3446  const double pad_bottom,
3447  const double pad_top,
3448  const double pad_left,
3449  const double pad_right,
3450  const Point<3> new_center,
3451  const types::manifold_id polar_manifold_id,
3452  const types::manifold_id tfi_manifold_id,
3453  const double L,
3454  const unsigned int n_slices,
3455  const bool colorize)
3456  {
3457  Triangulation<2> tria_2;
3458  plate_with_a_hole(tria_2,
3459  inner_radius,
3460  outer_radius,
3461  pad_bottom,
3462  pad_top,
3463  pad_left,
3464  pad_right,
3465  Point<2>(new_center[0], new_center[1]),
3466  polar_manifold_id,
3467  tfi_manifold_id,
3468  L,
3469  n_slices,
3470  colorize);
3471 
3472  // extrude to 3D
3473  extrude_triangulation(tria_2, n_slices, L, tria, true);
3474 
3475  // shift in Z direction to match specified center
3476  GridTools::shift(Point<3>(0, 0, new_center[2] - L / 2.), tria);
3477 
3478  // set up the new manifolds
3479  const Tensor<1, 3> direction{{0.0, 0.0, 1.0}};
3480  const CylindricalManifold<3> cylindrical_manifold(
3481  direction,
3482  /*axial_point*/ new_center);
3483  TransfiniteInterpolationManifold<3> inner_manifold;
3484  inner_manifold.initialize(tria);
3485  tria.set_manifold(polar_manifold_id, cylindrical_manifold);
3486  tria.set_manifold(tfi_manifold_id, inner_manifold);
3487  }
3488 
3489 
3490 
3491  template <>
3493  const double shell_region_width,
3494  const unsigned int n_shells,
3495  const double skewness,
3496  const bool colorize)
3497  {
3498  Assert(0.0 <= shell_region_width && shell_region_width < 0.05,
3499  ExcMessage("The width of the shell region must be less than 0.05 "
3500  "(and preferably close to 0.03)"));
3501  const types::manifold_id polar_manifold_id = 0;
3502  const types::manifold_id tfi_manifold_id = 1;
3503 
3504  // We begin by setting up a grid that is 4 by 22 cells. While not
3505  // squares, these have pretty good aspect ratios.
3506  Triangulation<2> bulk_tria;
3508  {22u, 4u},
3509  Point<2>(0.0, 0.0),
3510  Point<2>(2.2, 0.41));
3511  // bulk_tria now looks like this:
3512  //
3513  // +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
3514  // | | | | | | | | | | | | | | | | | | | | | | |
3515  // +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
3516  // | |XX|XX| | | | | | | | | | | | | | | | | | | |
3517  // +--+--O--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
3518  // | |XX|XX| | | | | | | | | | | | | | | | | | | |
3519  // +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
3520  // | | | | | | | | | | | | | | | | | | | | | | |
3521  // +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
3522  //
3523  // Note that these cells are not quite squares: they are all 0.1 by
3524  // 0.1025.
3525  //
3526  // The next step is to remove the cells marked with XXs: we will place
3527  // the grid around the cylinder there later. The next loop does two
3528  // things:
3529  // 1. Determines which cells need to be removed from the Triangulation
3530  // (i.e., find the cells marked with XX in the picture).
3531  // 2. Finds the location of the vertex marked with 'O' and uses that to
3532  // calculate the shift vector for aligning cylinder_tria with
3533  // tria_without_cylinder.
3534  std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3535  Tensor<1, 2> cylinder_triangulation_offset;
3536  for (const auto &cell : bulk_tria.active_cell_iterators())
3537  {
3538  if ((cell->center() - Point<2>(0.2, 0.2)).norm() < 0.15)
3539  cells_to_remove.insert(cell);
3540 
3541  if (cylinder_triangulation_offset == Tensor<1, 2>())
3542  {
3543  for (unsigned int vertex_n = 0;
3544  vertex_n < GeometryInfo<2>::vertices_per_cell;
3545  ++vertex_n)
3546  if (cell->vertex(vertex_n) == Point<2>())
3547  {
3548  // cylinder_tria is centered at zero, so we need to
3549  // shift it up and to the right by two cells:
3550  cylinder_triangulation_offset =
3551  2.0 * (cell->vertex(3) - Point<2>());
3552  break;
3553  }
3554  }
3555  }
3556  Triangulation<2> tria_without_cylinder;
3558  bulk_tria, cells_to_remove, tria_without_cylinder);
3559 
3560  // set up the cylinder triangulation. Note that this function sets the
3561  // manifold ids of the interior boundary cells to 0
3562  // (polar_manifold_id).
3563  Triangulation<2> cylinder_tria;
3565  0.05 + shell_region_width,
3566  0.41 / 4.0);
3567  // The bulk cells are not quite squares, so we need to move the left
3568  // and right sides of cylinder_tria inwards so that it fits in
3569  // bulk_tria:
3570  for (const auto &cell : cylinder_tria.active_cell_iterators())
3571  for (unsigned int vertex_n = 0;
3572  vertex_n < GeometryInfo<2>::vertices_per_cell;
3573  ++vertex_n)
3574  {
3575  if (std::abs(cell->vertex(vertex_n)[0] - -0.41 / 4.0) < 1e-10)
3576  cell->vertex(vertex_n)[0] = -0.1;
3577  else if (std::abs(cell->vertex(vertex_n)[0] - 0.41 / 4.0) < 1e-10)
3578  cell->vertex(vertex_n)[0] = 0.1;
3579  }
3580 
3581  // Assign interior manifold ids to be the TFI id.
3582  for (const auto &cell : cylinder_tria.active_cell_iterators())
3583  {
3584  cell->set_manifold_id(tfi_manifold_id);
3585  for (unsigned int face_n = 0; face_n < GeometryInfo<2>::faces_per_cell;
3586  ++face_n)
3587  if (!cell->face(face_n)->at_boundary())
3588  cell->face(face_n)->set_manifold_id(tfi_manifold_id);
3589  }
3590  if (0.0 < shell_region_width)
3591  {
3592  Assert(0 < n_shells,
3593  ExcMessage("If the shell region has positive width then "
3594  "there must be at least one shell."));
3595  Triangulation<2> shell_tria;
3597  Point<2>(),
3598  0.05,
3599  0.05 + shell_region_width,
3600  n_shells,
3601  skewness,
3602  8);
3603 
3604  // Make the tolerance as large as possible since these cells can
3605  // be quite close together
3606  const double vertex_tolerance =
3607  std::min(internal::minimal_vertex_distance(shell_tria),
3608  internal::minimal_vertex_distance(cylinder_tria)) *
3609  0.5;
3610 
3611  shell_tria.set_all_manifold_ids(polar_manifold_id);
3612  Triangulation<2> temp;
3614  shell_tria, cylinder_tria, temp, vertex_tolerance, true);
3615  cylinder_tria = std::move(temp);
3616  }
3617  GridTools::shift(cylinder_triangulation_offset, cylinder_tria);
3618 
3619  // Compute the tolerance again, since the shells may be very close to
3620  // each-other:
3621  const double vertex_tolerance =
3622  std::min(internal::minimal_vertex_distance(tria_without_cylinder),
3623  internal::minimal_vertex_distance(cylinder_tria)) /
3624  10;
3626  tria_without_cylinder, cylinder_tria, tria, vertex_tolerance, true);
3627 
3628  // Move the vertices in the middle of the faces of cylinder_tria slightly
3629  // to give a better mesh quality. We have to balance the quality of these
3630  // cells with the quality of the outer cells (initially rectangles). For
3631  // constant radial distance, we would place them at the distance 0.1 *
3632  // sqrt(2.) from the center. In case the shell region width is more than
3633  // 0.1/6., we choose to place them at 0.1 * 4./3. from the center, which
3634  // ensures that the shortest edge of the outer cells is 2./3. of the
3635  // original length. If the shell region width is less, we make the edge
3636  // length of the inner part and outer part (in the shorter x direction)
3637  // the same.
3638  {
3639  const double shift =
3640  std::min(0.125 + shell_region_width * 0.5, 0.1 * 4. / 3.);
3641  for (const auto &cell : tria.active_cell_iterators())
3642  for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v)
3643  if (cell->vertex(v).distance(Point<2>(0.1, 0.205)) < 1e-10)
3644  cell->vertex(v) = Point<2>(0.2 - shift, 0.205);
3645  else if (cell->vertex(v).distance(Point<2>(0.3, 0.205)) < 1e-10)
3646  cell->vertex(v) = Point<2>(0.2 + shift, 0.205);
3647  else if (cell->vertex(v).distance(Point<2>(0.2, 0.1025)) < 1e-10)
3648  cell->vertex(v) = Point<2>(0.2, 0.2 - shift);
3649  else if (cell->vertex(v).distance(Point<2>(0.2, 0.3075)) < 1e-10)
3650  cell->vertex(v) = Point<2>(0.2, 0.2 + shift);
3651  }
3652 
3653  // Ensure that all manifold ids on a polar cell really are set to the
3654  // polar manifold id:
3655  for (const auto &cell : tria.active_cell_iterators())
3656  if (cell->manifold_id() == polar_manifold_id)
3657  cell->set_all_manifold_ids(polar_manifold_id);
3658 
3659  // Ensure that all other manifold ids (including the interior faces
3660  // opposite the cylinder) are set to the flat manifold id:
3661  for (const auto &cell : tria.active_cell_iterators())
3662  if (cell->manifold_id() != polar_manifold_id &&
3663  cell->manifold_id() != tfi_manifold_id)
3664  cell->set_all_manifold_ids(numbers::flat_manifold_id);
3665 
3666  // We need to calculate the current center so that we can move it later:
3667  // to start get a unique list of (points to) vertices on the cylinder
3668  std::vector<Point<2> *> cylinder_pointers;
3669  for (const auto &face : tria.active_face_iterators())
3670  if (face->manifold_id() == polar_manifold_id)
3671  {
3672  cylinder_pointers.push_back(&face->vertex(0));
3673  cylinder_pointers.push_back(&face->vertex(1));
3674  }
3675  // de-duplicate
3676  std::sort(cylinder_pointers.begin(), cylinder_pointers.end());
3677  cylinder_pointers.erase(std::unique(cylinder_pointers.begin(),
3678  cylinder_pointers.end()),
3679  cylinder_pointers.end());
3680 
3681  // find the current center...
3682  Point<2> center;
3683  for (const Point<2> *const ptr : cylinder_pointers)
3684  center += *ptr / double(cylinder_pointers.size());
3685 
3686  // and recenter at (0.2, 0.2)
3687  for (Point<2> *const ptr : cylinder_pointers)
3688  *ptr += Point<2>(0.2, 0.2) - center;
3689 
3690  // attach manifolds
3691  PolarManifold<2> polar_manifold(Point<2>(0.2, 0.2));
3692  tria.set_manifold(polar_manifold_id, polar_manifold);
3693  TransfiniteInterpolationManifold<2> inner_manifold;
3694  inner_manifold.initialize(tria);
3695  tria.set_manifold(tfi_manifold_id, inner_manifold);
3696 
3697  if (colorize)
3698  for (const auto &face : tria.active_face_iterators())
3699  if (face->at_boundary())
3700  {
3701  const Point<2> center = face->center();
3702  // left side
3703  if (std::abs(center[0] - 0.0) < 1e-10)
3704  face->set_boundary_id(0);
3705  // right side
3706  else if (std::abs(center[0] - 2.2) < 1e-10)
3707  face->set_boundary_id(1);
3708  // cylinder boundary
3709  else if (face->manifold_id() == polar_manifold_id)
3710  face->set_boundary_id(2);
3711  // sides of channel
3712  else
3713  {
3714  Assert(std::abs(center[1] - 0.00) < 1.0e-10 ||
3715  std::abs(center[1] - 0.41) < 1.0e-10,
3716  ExcInternalError());
3717  face->set_boundary_id(3);
3718  }
3719  }
3720  }
3721 
3722 
3723 
3724  template <>
3726  const double shell_region_width,
3727  const unsigned int n_shells,
3728  const double skewness,
3729  const bool colorize)
3730  {
3731  Triangulation<2> tria_2;
3733  tria_2, shell_region_width, n_shells, skewness, colorize);
3734  extrude_triangulation(tria_2, 5, 0.41, tria, true);
3735 
3736  // set up the new 3D manifolds
3737  const types::manifold_id cylindrical_manifold_id = 0;
3738  const types::manifold_id tfi_manifold_id = 1;
3739  const PolarManifold<2> *const m_ptr =
3740  dynamic_cast<const PolarManifold<2> *>(
3741  &tria_2.get_manifold(cylindrical_manifold_id));
3742  Assert(m_ptr != nullptr, ExcInternalError());
3743  const Point<3> axial_point(m_ptr->center[0], m_ptr->center[1], 0.0);
3744  const Tensor<1, 3> direction{{0.0, 0.0, 1.0}};
3745 
3746  const CylindricalManifold<3> cylindrical_manifold(direction, axial_point);
3747  TransfiniteInterpolationManifold<3> inner_manifold;
3748  inner_manifold.initialize(tria);
3749  tria.set_manifold(cylindrical_manifold_id, cylindrical_manifold);
3750  tria.set_manifold(tfi_manifold_id, inner_manifold);
3751 
3752  // From extrude_triangulation: since the maximum boundary id of tria_2 was
3753  // 3, the bottom boundary id is 4 and the top is 5: both are walls, so set
3754  // them to 3
3755  if (colorize)
3756  for (const auto &face : tria.active_face_iterators())
3757  if (face->boundary_id() == 4 || face->boundary_id() == 5)
3758  face->set_boundary_id(3);
3759  }
3760 
3761 
3762 
3763  template <int dim, int spacedim>
3764  void
3766  const std::vector<unsigned int> &sizes,
3767  const bool colorize)
3768  {
3770  Assert(dim > 1, ExcNotImplemented());
3771  Assert(dim < 4, ExcNotImplemented());
3772 
3773  // If there is a desire at some point to change the geometry of
3774  // the cells, this tensor can be made an argument to the function.
3775  Tensor<1, dim> dimensions;
3776  for (unsigned int d = 0; d < dim; ++d)
3777  dimensions[d] = 1.;
3778 
3779  std::vector<Point<spacedim>> points;
3780  unsigned int n_cells = 1;
3781  for (unsigned int i : GeometryInfo<dim>::face_indices())
3782  n_cells += sizes[i];
3783 
3784  std::vector<CellData<dim>> cells(n_cells);
3785  // Vertices of the center cell
3786  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
3787  {
3788  Point<spacedim> p;
3789  for (unsigned int d = 0; d < dim; ++d)
3790  p(d) = 0.5 * dimensions[d] *
3793  points.push_back(p);
3794  cells[0].vertices[i] = i;
3795  }
3796  cells[0].material_id = 0;
3797 
3798  // The index of the first cell of the leg.
3799  unsigned int cell_index = 1;
3800  // The legs of the cross
3801  for (const unsigned int face : GeometryInfo<dim>::face_indices())
3802  {
3803  const unsigned int oface = GeometryInfo<dim>::opposite_face[face];
3804  const unsigned int dir = GeometryInfo<dim>::unit_normal_direction[face];
3805 
3806  // We are moving in the direction of face
3807  for (unsigned int j = 0; j < sizes[face]; ++j, ++cell_index)
3808  {
3809  const unsigned int last_cell = (j == 0) ? 0U : (cell_index - 1);
3810 
3811  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
3812  ++v)
3813  {
3814  const unsigned int cellv =
3816  const unsigned int ocellv =
3818  // First the vertices which already exist
3819  cells[cell_index].vertices[ocellv] =
3820  cells[last_cell].vertices[cellv];
3821 
3822  // Now the new vertices
3823  cells[cell_index].vertices[cellv] = points.size();
3824 
3825  Point<spacedim> p = points[cells[cell_index].vertices[ocellv]];
3827  dimensions[dir];
3828  points.push_back(p);
3829  }
3830  cells[cell_index].material_id = (colorize) ? (face + 1U) : 0U;
3831  }
3832  }
3833  tria.create_triangulation(points, cells, SubCellData());
3834  }
3835 
3836 
3837  template <>
3838  void
3839  hyper_cube_slit(Triangulation<1> &, const double, const double, const bool)
3840  {
3841  Assert(false, ExcNotImplemented());
3842  }
3843 
3844 
3845 
3846  template <>
3848  const double,
3849  const double,
3850  const double,
3851  const bool)
3852  {
3853  Assert(false, ExcNotImplemented());
3854  }
3855 
3856 
3857 
3858  template <>
3859  void hyper_L(Triangulation<1> &, const double, const double, const bool)
3860  {
3861  Assert(false, ExcNotImplemented());
3862  }
3863 
3864 
3865 
3866  template <>
3867  void
3868  hyper_ball(Triangulation<1> &, const Point<1> &, const double, const bool)
3869  {
3870  Assert(false, ExcNotImplemented());
3871  }
3872 
3873 
3874 
3875  template <>
3876  void cylinder(Triangulation<1> &, const double, const double)
3877  {
3878  Assert(false, ExcNotImplemented());
3879  }
3880 
3881 
3882 
3883  template <>
3884  void
3885  truncated_cone(Triangulation<1> &, const double, const double, const double)
3886  {
3887  Assert(false, ExcNotImplemented());
3888  }
3889 
3890 
3891 
3892  template <>
3894  const Point<1> &,
3895  const double,
3896  const double,
3897  const unsigned int,
3898  const bool)
3899  {
3900  Assert(false, ExcNotImplemented());
3901  }
3902 
3903  template <>
3905  const double,
3906  const double,
3907  const double,
3908  const unsigned int,
3909  const unsigned int)
3910  {
3911  Assert(false, ExcNotImplemented());
3912  }
3913 
3914 
3915  template <>
3916  void quarter_hyper_ball(Triangulation<1> &, const Point<1> &, const double)
3917  {
3918  Assert(false, ExcNotImplemented());
3919  }
3920 
3921 
3922  template <>
3923  void half_hyper_ball(Triangulation<1> &, const Point<1> &, const double)
3924  {
3925  Assert(false, ExcNotImplemented());
3926  }
3927 
3928 
3929  template <>
3931  const Point<1> &,
3932  const double,
3933  const double,
3934  const unsigned int,
3935  const bool)
3936  {
3937  Assert(false, ExcNotImplemented());
3938  }
3939 
3940  template <>
3942  const Point<1> &,
3943  const double,
3944  const double,
3945  const unsigned int,
3946  const bool)
3947  {
3948  Assert(false, ExcNotImplemented());
3949  }
3950 
3951  template <>
3953  const double left,
3954  const double right,
3955  const double thickness,
3956  const bool colorize)
3957  {
3958  Assert(left < right,
3959  ExcMessage("Invalid left-to-right bounds of enclosed hypercube"));
3960 
3961  std::vector<Point<2>> vertices(16);
3962  double coords[4];
3963  coords[0] = left - thickness;
3964  coords[1] = left;
3965  coords[2] = right;
3966  coords[3] = right + thickness;
3967 
3968  unsigned int k = 0;
3969  for (const double y : coords)
3970  for (const double x : coords)
3971  vertices[k++] = Point<2>(x, y);
3972 
3973  const types::material_id materials[9] = {5, 4, 6, 1, 0, 2, 9, 8, 10};
3974 
3975  std::vector<CellData<2>> cells(9);
3976  k = 0;
3977  for (unsigned int i0 = 0; i0 < 3; ++i0)
3978  for (unsigned int i1 = 0; i1 < 3; ++i1)
3979  {
3980  cells[k].vertices[0] = i1 + 4 * i0;
3981  cells[k].vertices[1] = i1 + 4 * i0 + 1;
3982  cells[k].vertices[2] = i1 + 4 * i0 + 4;
3983  cells[k].vertices[3] = i1 + 4 * i0 + 5;
3984  if (colorize)
3985  cells[k].material_id = materials[k];
3986  ++k;
3987  }
3988  tria.create_triangulation(vertices,
3989  cells,
3990  SubCellData()); // no boundary information
3991  }
3992 
3993 
3994 
3995  // Implementation for 2D only
3996  template <>
3997  void hyper_cube_slit(Triangulation<2> &tria,
3998  const double left,
3999  const double right,
4000  const bool colorize)
4001  {
4002  const double rl2 = (right + left) / 2;
4003  const Point<2> vertices[10] = {Point<2>(left, left),
4004  Point<2>(rl2, left),
4005  Point<2>(rl2, rl2),
4006  Point<2>(left, rl2),
4007  Point<2>(right, left),
4008  Point<2>(right, rl2),
4009  Point<2>(rl2, right),
4010  Point<2>(left, right),
4011  Point<2>(right, right),
4012  Point<2>(rl2, left)};
4013  const int cell_vertices[4][4] = {{0, 1, 3, 2},
4014  {9, 4, 2, 5},
4015  {3, 2, 7, 6},
4016  {2, 5, 6, 8}};
4017  std::vector<CellData<2>> cells(4, CellData<2>());
4018  for (unsigned int i = 0; i < 4; ++i)
4019  {
4020  for (unsigned int j = 0; j < 4; ++j)
4021  cells[i].vertices[j] = cell_vertices[i][j];
4022  cells[i].material_id = 0;
4023  }
4024  tria.create_triangulation(std::vector<Point<2>>(std::begin(vertices),
4025  std::end(vertices)),
4026  cells,
4027  SubCellData()); // no boundary information
4028 
4029  if (colorize)
4030  {
4031  Triangulation<2>::cell_iterator cell = tria.begin();
4032  cell->face(1)->set_boundary_id(1);
4033  ++cell;
4034  cell->face(0)->set_boundary_id(2);
4035  }
4036  }
4037 
4038 
4039 
4040  template <>
4041  void truncated_cone(Triangulation<2> &triangulation,
4042  const double radius_0,
4043  const double radius_1,
4044  const double half_length)
4045  {
4046  Point<2> vertices_tmp[4];
4047 
4048  vertices_tmp[0] = Point<2>(-half_length, -radius_0);
4049  vertices_tmp[1] = Point<2>(half_length, -radius_1);
4050  vertices_tmp[2] = Point<2>(-half_length, radius_0);
4051  vertices_tmp[3] = Point<2>(half_length, radius_1);
4052 
4053  const std::vector<Point<2>> vertices(std::begin(vertices_tmp),
4054  std::end(vertices_tmp));
4055  unsigned int cell_vertices[1][GeometryInfo<2>::vertices_per_cell];
4056 
4057  for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
4058  cell_vertices[0][i] = i;
4059 
4060  std::vector<CellData<2>> cells(1, CellData<2>());
4061 
4062  for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
4063  cells[0].vertices[i] = cell_vertices[0][i];
4064 
4065  cells[0].material_id = 0;
4066  triangulation.create_triangulation(vertices, cells, SubCellData());
4067 
4068  Triangulation<2>::cell_iterator cell = triangulation.begin();
4069 
4070  cell->face(0)->set_boundary_id(1);
4071  cell->face(1)->set_boundary_id(2);
4072 
4073  for (unsigned int i = 2; i < 4; ++i)
4074  cell->face(i)->set_boundary_id(0);
4075  }
4076 
4077 
4078 
4079  // Implementation for 2D only
4080  template <>
4081  void hyper_L(Triangulation<2> &tria,
4082  const double a,
4083  const double b,
4084  const bool colorize)
4085  {
4086  const Point<2> vertices[8] = {Point<2>(a, a),
4087  Point<2>((a + b) / 2, a),
4088  Point<2>(b, a),
4089  Point<2>(a, (a + b) / 2),
4090  Point<2>((a + b) / 2, (a + b) / 2),
4091  Point<2>(b, (a + b) / 2),
4092  Point<2>(a, b),
4093  Point<2>((a + b) / 2, b)};
4094  const int cell_vertices[3][4] = {{0, 1, 3, 4}, {1, 2, 4, 5}, {3, 4, 6, 7}};
4095 
4096  std::vector<CellData<2>> cells(3, CellData<2>());
4097 
4098  for (unsigned int i = 0; i < 3; ++i)
4099  {
4100  for (unsigned int j = 0; j < 4; ++j)
4101  cells[i].vertices[j] = cell_vertices[i][j];
4102  cells[i].material_id = 0;
4103  }
4104 
4105  tria.create_triangulation(std::vector<Point<2>>(std::begin(vertices),
4106  std::end(vertices)),
4107  cells,
4108  SubCellData());
4109 
4110  if (colorize)
4111  {
4112  Triangulation<2>::cell_iterator cell = tria.begin();
4113 
4114  cell->face(0)->set_boundary_id(0);
4115  cell->face(2)->set_boundary_id(1);
4116  cell++;
4117 
4118  cell->face(1)->set_boundary_id(2);
4119  cell->face(2)->set_boundary_id(1);
4120  cell->face(3)->set_boundary_id(3);
4121  cell++;
4122 
4123  cell->face(0)->set_boundary_id(0);
4124  cell->face(1)->set_boundary_id(4);
4125  cell->face(3)->set_boundary_id(5);
4126  }
4127  }
4128 
4129 
4130 
4131  template <int dim, int spacedim>
4132  void
4134  const std::vector<unsigned int> &repetitions,
4135  const Point<dim> & bottom_left,
4136  const Point<dim> & top_right,
4137  const std::vector<int> & n_cells_to_remove)
4138  {
4139  Assert(dim > 1, ExcNotImplemented());
4140  // Check the consistency of the dimensions provided.
4141  AssertDimension(repetitions.size(), dim);
4142  AssertDimension(n_cells_to_remove.size(), dim);
4143  for (unsigned int d = 0; d < dim; ++d)
4144  {
4145  Assert(std::fabs(n_cells_to_remove[d]) <= repetitions[d],
4146  ExcMessage("Attempting to cut away too many cells."));
4147  }
4148  // Create the domain to be cut
4149  Triangulation<dim, spacedim> rectangle;
4151  repetitions,
4152  bottom_left,
4153  top_right);
4154  // compute the vertex of the cut step, we will cut according to the
4155  // location of the cartesian coordinates of the cell centers
4156  std::array<double, dim> h;
4157  Point<dim> cut_step;
4158  for (unsigned int d = 0; d < dim; ++d)
4159  {
4160  // mesh spacing in each direction in cartesian coordinates
4161  h[d] = (top_right[d] - bottom_left[d]) / repetitions[d];
4162  // left to right, bottom to top, front to back
4163  if (n_cells_to_remove[d] >= 0)
4164  {
4165  // cartesian coordinates of vertex location
4166  cut_step[d] =
4167  h[d] * std::fabs(n_cells_to_remove[d]) + bottom_left[d];
4168  }
4169  // right to left, top to bottom, back to front
4170  else
4171  {
4172  cut_step[d] = top_right[d] - h[d] * std::fabs(n_cells_to_remove[d]);
4173  }
4174  }
4175 
4176 
4177  // compute cells to remove
4178  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>
4179  cells_to_remove;
4180  std::copy_if(
4181  rectangle.active_cell_iterators().begin(),
4182  rectangle.active_cell_iterators().end(),
4183  std::inserter(cells_to_remove, cells_to_remove.end()),
4184  [&](
4186  -> bool {
4187  for (unsigned int d = 0; d < dim; ++d)
4188  if ((n_cells_to_remove[d] > 0 && cell->center()[d] >= cut_step[d]) ||
4189  (n_cells_to_remove[d] < 0 && cell->center()[d] <= cut_step[d]))
4190  return false;
4191 
4192  return true;
4193  });
4194 
4196  cells_to_remove,
4197  tria);
4198  }
4199 
4200 
4201 
4202  // Implementation for 2D only
4203  template <>
4204  void hyper_ball(Triangulation<2> &tria,
4205  const Point<2> & p,
4206  const double radius,
4207  const bool internal_manifolds)
4208  {
4209  // equilibrate cell sizes at
4210  // transition from the inner part
4211  // to the radial cells
4212  const double a = 1. / (1 + std::sqrt(2.0));
4213  const Point<2> vertices[8] = {
4214  p + Point<2>(-1, -1) * (radius / std::sqrt(2.0)),
4215  p + Point<2>(+1, -1) * (radius / std::sqrt(2.0)),
4216  p + Point<2>(-1, -1) * (radius / std::sqrt(2.0) * a),
4217  p + Point<2>(+1, -1) * (radius / std::sqrt(2.0) * a),
4218  p + Point<2>(-1, +1) * (radius / std::sqrt(2.0) * a),
4219  p + Point<2>(+1, +1) * (radius / std::sqrt(2.0) * a),
4220  p + Point<2>(-1, +1) * (radius / std::sqrt(2.0)),
4221  p + Point<2>(+1, +1) * (radius / std::sqrt(2.0))};
4222 
4223  const int cell_vertices[5][4] = {
4224  {0, 1, 2, 3}, {0, 2, 6, 4}, {2, 3, 4, 5}, {1, 7, 3, 5}, {6, 4, 7, 5}};
4225 
4226  std::vector<CellData<2>> cells(5, CellData<2>());
4227 
4228  for (unsigned int i = 0; i < 5; ++i)
4229  {
4230  for (unsigned int j = 0; j < 4; ++j)
4231  cells[i].vertices[j] = cell_vertices[i][j];
4232  cells[i].material_id = 0;
4233  cells[i].manifold_id = i == 2 ? 1 : numbers::flat_manifold_id;
4234  }
4235 
4236  tria.create_triangulation(std::vector<Point<2>>(std::begin(vertices),
4237  std::end(vertices)),
4238  cells,
4239  SubCellData()); // no boundary information
4241  tria.set_manifold(0, SphericalManifold<2>(p));
4242  if (internal_manifolds)
4243  tria.set_manifold(1, SphericalManifold<2>(p));
4244  }
4245 
4246 
4247 
4248  template <>
4249  void hyper_shell(Triangulation<2> & tria,
4250  const Point<2> & center,
4251  const double inner_radius,
4252  const double outer_radius,
4253  const unsigned int n_cells,
4254  const bool colorize)
4255  {
4256  Assert((inner_radius > 0) && (inner_radius < outer_radius),
4257  ExcInvalidRadii());
4258 
4259  const double pi = numbers::PI;
4260 
4261  // determine the number of cells
4262  // for the grid. if not provided by
4263  // the user determine it such that
4264  // the length of each cell on the
4265  // median (in the middle between
4266  // the two circles) is equal to its
4267  // radial extent (which is the
4268  // difference between the two
4269  // radii)
4270  const unsigned int N =
4271  (n_cells == 0 ? static_cast<unsigned int>(
4272  std::ceil((2 * pi * (outer_radius + inner_radius) / 2) /
4273  (outer_radius - inner_radius))) :
4274  n_cells);
4275 
4276  // set up N vertices on the
4277  // outer and N vertices on
4278  // the inner circle. the
4279  // first N ones are on the
4280  // outer one, and all are
4281  // numbered counter-clockwise
4282  std::vector<Point<2>> vertices(2 * N);
4283  for (unsigned int i = 0; i < N; ++i)
4284  {
4285  vertices[i] =
4286  Point<2>(std::cos(2 * pi * i / N), std::sin(2 * pi * i / N)) *
4287  outer_radius;
4288  vertices[i + N] = vertices[i] * (inner_radius / outer_radius);
4289 
4290  vertices[i] += center;
4291  vertices[i + N] += center;
4292  }
4293 
4294  std::vector<CellData<2>> cells(N, CellData<2>());
4295 
4296  for (unsigned int i = 0; i < N; ++i)
4297  {
4298  cells[i].vertices[0] = i;
4299  cells[i].vertices[1] = (i + 1) % N;
4300  cells[i].vertices[2] = N + i;
4301  cells[i].vertices[3] = N + ((i + 1) % N);
4302 
4303  cells[i].material_id = 0;
4304  }
4305 
4306  tria.create_triangulation(vertices, cells, SubCellData());
4307 
4308  if (colorize)
4309  colorize_hyper_shell(tria, center, inner_radius, outer_radius);
4310 
4311  tria.set_all_manifold_ids(0);
4312  tria.set_manifold(0, SphericalManifold<2>(center));
4313  }
4314 
4315 
4316 
4317  template <int dim>
4318  void
4320  const Point<dim> & inner_center,
4321  const Point<dim> & outer_center,
4322  const double inner_radius,
4323  const double outer_radius,
4324  const unsigned int n_cells)
4325  {
4327  tria, outer_center, inner_radius, outer_radius, n_cells, true);
4328 
4329  // check the consistency of the dimensions provided
4330  Assert(
4331  outer_radius - inner_radius > outer_center.distance(inner_center),
4333  "The inner radius is greater than or equal to the outer radius plus eccentricity."));
4334 
4335  // shift nodes along the inner boundary according to the position of
4336  // inner_circle
4337  std::set<Point<dim> *> vertices_to_move;
4338 
4339  for (const auto &face : tria.active_face_iterators())
4340  if (face->boundary_id() == 0)
4341  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
4342  vertices_to_move.insert(&face->vertex(v));
4343 
4344  const auto shift = inner_center - outer_center;
4345  for (const auto &p : vertices_to_move)
4346  (*p) += shift;
4347 
4348  // the original hyper_shell function assigns the same manifold id
4349  // to all cells and faces. Set all manifolds ids to a different
4350  // value (2), then use boundary ids to assign different manifolds to
4351  // the inner (0) and outer manifolds (1). Use a transfinite manifold
4352  // for all faces and cells aside from the boundaries.
4353  tria.set_all_manifold_ids(2);
4355 
4356  SphericalManifold<dim> inner_manifold(inner_center);
4357  SphericalManifold<dim> outer_manifold(outer_center);
4358 
4360  transfinite.initialize(tria);
4361 
4362  tria.set_manifold(0, inner_manifold);
4363  tria.set_manifold(1, outer_manifold);
4364  tria.set_manifold(2, transfinite);
4365  }
4366 
4367 
4368 
4369  // Implementation for 2D only
4370  template <>
4371  void cylinder(Triangulation<2> &tria,
4372  const double radius,
4373  const double half_length)
4374  {
4375  Point<2> p1(-half_length, -radius);
4376  Point<2> p2(half_length, radius);
4377 
4378  hyper_rectangle(tria, p1, p2, true);
4379 
4382  while (f != end)
4383  {
4384  switch (f->boundary_id())
4385  {
4386  case 0:
4387  f->set_boundary_id(1);
4388  break;
4389  case 1:
4390  f->set_boundary_id(2);
4391  break;
4392  default:
4393  f->set_boundary_id(0);
4394  break;
4395  }
4396  ++f;
4397  }
4398  }
4399 
4400 
4401 
4402  // Implementation for 2D only
4403  template <>
4405  const double,
4406  const double,
4407  const double,
4408  const unsigned int,
4409  const unsigned int)
4410  {
4411  Assert(false, ExcNotImplemented());
4412  }
4413 
4414 
4415  template <>
4417  const Point<2> & p,
4418  const double radius)
4419  {
4420  const unsigned int dim = 2;
4421 
4422  // equilibrate cell sizes at
4423  // transition from the inner part
4424  // to the radial cells
4425  const Point<dim> vertices[7] = {p + Point<dim>(0, 0) * radius,
4426  p + Point<dim>(+1, 0) * radius,
4427  p + Point<dim>(+1, 0) * (radius / 2),
4428  p + Point<dim>(0, +1) * (radius / 2),
4429  p + Point<dim>(+1, +1) *
4430  (radius / (2 * std::sqrt(2.0))),
4431  p + Point<dim>(0, +1) * radius,
4432  p + Point<dim>(+1, +1) *
4433  (radius / std::sqrt(2.0))};
4434 
4435  const int cell_vertices[3][4] = {{0, 2, 3, 4}, {1, 6, 2, 4}, {5, 3, 6, 4}};
4436 
4437  std::vector<CellData<dim>> cells(3, CellData<dim>());
4438 
4439  for (unsigned int i = 0; i < 3; ++i)
4440  {
4441  for (unsigned int j = 0; j < 4; ++j)
4442  cells[i].vertices[j] = cell_vertices[i][j];
4443  cells[i].material_id = 0;
4444  }
4445 
4446  tria.create_triangulation(std::vector<Point<dim>>(std::begin(vertices),
4447  std::end(vertices)),
4448  cells,
4449  SubCellData()); // no boundary information
4450 
4453 
4455 
4456  while (cell != end)
4457  {
4458  for (unsigned int i : GeometryInfo<dim>::face_indices())
4459  {
4460  if (cell->face(i)->boundary_id() ==
4462  continue;
4463 
4464  // If one the components is the same as the respective
4465  // component of the center, then this is part of the plane
4466  if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius ||
4467  cell->face(i)->center()(1) < p(1) + 1.e-5 * radius)
4468  {
4469  cell->face(i)->set_boundary_id(1);
4470  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
4471  }
4472  }
4473  ++cell;
4474  }
4475  tria.set_manifold(0, SphericalManifold<2>(p));
4476  }
4477 
4478 
4479  template <>
4480  void half_hyper_ball(Triangulation<2> &tria,
4481  const Point<2> & p,
4482  const double radius)
4483  {
4484  // equilibrate cell sizes at
4485  // transition from the inner part
4486  // to the radial cells
4487  const double a = 1. / (1 + std::sqrt(2.0));
4488  const Point<2> vertices[8] = {
4489  p + Point<2>(0, -1) * radius,
4490  p + Point<2>(+1, -1) * (radius / std::sqrt(2.0)),
4491  p + Point<2>(0, -1) * (radius / std::sqrt(2.0) * a),
4492  p + Point<2>(+1, -1) * (radius / std::sqrt(2.0) * a),
4493  p + Point<2>(0, +1) * (radius / std::sqrt(2.0) * a),
4494  p + Point<2>(+1, +1) * (radius / std::sqrt(2.0) * a),
4495  p + Point<2>(0, +1) * radius,
4496  p + Point<2>(+1, +1) * (radius / std::sqrt(2.0))};
4497 
4498  const int cell_vertices[5][4] = {{0, 1, 2, 3},
4499  {2, 3, 4, 5},
4500  {1, 7, 3, 5},
4501  {6, 4, 7, 5}};
4502 
4503  std::vector<CellData<2>> cells(4, CellData<2>());
4504 
4505  for (unsigned int i = 0; i < 4; ++i)
4506  {
4507  for (unsigned int j = 0; j < 4; ++j)
4508  cells[i].vertices[j] = cell_vertices[i][j];
4509  cells[i].material_id = 0;
4510  }
4511 
4512  tria.create_triangulation(std::vector<Point<2>>(std::begin(vertices),
4513  std::end(vertices)),
4514  cells,
4515  SubCellData()); // no boundary information
4516 
4517  Triangulation<2>::cell_iterator cell = tria.begin();
4518  Triangulation<2>::cell_iterator end = tria.end();
4519 
4521 
4522  while (cell != end)
4523  {
4524  for (unsigned int i : GeometryInfo<2>::face_indices())
4525  {
4526  if (cell->face(i)->boundary_id() ==
4528  continue;
4529 
4530  // If x is zero, then this is part of the plane
4531  if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius)
4532  {
4533  cell->face(i)->set_boundary_id(1);
4534  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
4535  }
4536  }
4537  ++cell;
4538  }
4539  tria.set_manifold(0, SphericalManifold<2>(p));
4540  }
4541 
4542 
4543 
4544  // Implementation for 2D only
4545  template <>
4546  void half_hyper_shell(Triangulation<2> & tria,
4547  const Point<2> & center,
4548  const double inner_radius,
4549  const double outer_radius,
4550  const unsigned int n_cells,
4551  const bool colorize)
4552  {
4553  Assert((inner_radius > 0) && (inner_radius < outer_radius),
4554  ExcInvalidRadii());
4555 
4556  const double pi = numbers::PI;
4557  // determine the number of cells
4558  // for the grid. if not provided by
4559  // the user determine it such that
4560  // the length of each cell on the
4561  // median (in the middle between
4562  // the two circles) is equal to its
4563  // radial extent (which is the
4564  // difference between the two
4565  // radii)
4566  const unsigned int N =
4567  (n_cells == 0 ? static_cast<unsigned int>(
4568  std::ceil((pi * (outer_radius + inner_radius) / 2) /
4569  (outer_radius - inner_radius))) :
4570  n_cells);
4571 
4572  // set up N+1 vertices on the
4573  // outer and N+1 vertices on
4574  // the inner circle. the
4575  // first N+1 ones are on the
4576  // outer one, and all are
4577  // numbered counter-clockwise
4578  std::vector<Point<2>> vertices(2 * (N + 1));
4579  for (unsigned int i = 0; i <= N; ++i)
4580  {
4581  // enforce that the x-coordinates
4582  // of the first and last point of
4583  // each half-circle are exactly
4584  // zero (contrary to what we may
4585  // compute using the imprecise
4586  // value of pi)
4587  vertices[i] =
4588  Point<2>(((i == 0) || (i == N) ? 0 : std::cos(pi * i / N - pi / 2)),
4589  std::sin(pi * i / N - pi / 2)) *
4590  outer_radius;
4591  vertices[i + N + 1] = vertices[i] * (inner_radius / outer_radius);
4592 
4593  vertices[i] += center;
4594  vertices[i + N + 1] += center;
4595  }
4596 
4597 
4598  std::vector<CellData<2>> cells(N, CellData<2>());
4599 
4600  for (unsigned int i = 0; i < N; ++i)
4601  {
4602  cells[i].vertices[0] = i;
4603  cells[i].vertices[1] = (i + 1) % (N + 1);
4604  cells[i].vertices[2] = N + 1 + i;
4605  cells[i].vertices[3] = N + 1 + ((i + 1) % (N + 1));
4606 
4607  cells[i].material_id = 0;
4608  }
4609 
4610  tria.create_triangulation(vertices, cells, SubCellData());
4611 
4612  if (colorize)
4613  {
4614  Triangulation<2>::cell_iterator cell = tria.begin();
4615  for (; cell != tria.end(); ++cell)
4616  {
4617  cell->face(2)->set_boundary_id(1);
4618  }
4619  tria.begin()->face(0)->set_boundary_id(3);
4620 
4621  tria.last()->face(1)->set_boundary_id(2);
4622  }
4623  tria.set_all_manifold_ids(0);
4624  tria.set_manifold(0, SphericalManifold<2>(center));
4625  }
4626 
4627 
4628  template <>
4630  const Point<2> & center,
4631  const double inner_radius,
4632  const double outer_radius,
4633  const unsigned int n_cells,
4634  const bool colorize)
4635  {
4636  Assert((inner_radius > 0) && (inner_radius < outer_radius),
4637  ExcInvalidRadii());
4638 
4639  const double pi = numbers::PI;
4640  // determine the number of cells
4641  // for the grid. if not provided by
4642  // the user determine it such that
4643  // the length of each cell on the
4644  // median (in the middle between
4645  // the two circles) is equal to its
4646  // radial extent (which is the
4647  // difference between the two
4648  // radii)
4649  const unsigned int N =
4650  (n_cells == 0 ? static_cast<unsigned int>(
4651  std::ceil((pi * (outer_radius + inner_radius) / 4) /
4652  (outer_radius - inner_radius))) :
4653  n_cells);
4654 
4655  // set up N+1 vertices on the
4656  // outer and N+1 vertices on
4657  // the inner circle. the
4658  // first N+1 ones are on the
4659  // outer one, and all are
4660  // numbered counter-clockwise
4661  std::vector<Point<2>> vertices(2 * (N + 1));
4662  for (unsigned int i = 0; i <= N; ++i)
4663  {
4664  // enforce that the x-coordinates
4665  // of the last point is exactly
4666  // zero (contrary to what we may
4667  // compute using the imprecise
4668  // value of pi)
4669  vertices[i] = Point<2>(((i == N) ? 0 : std::cos(pi * i / N / 2)),
4670  std::sin(pi * i / N / 2)) *
4671  outer_radius;
4672  vertices[i + N + 1] = vertices[i] * (inner_radius / outer_radius);
4673 
4674  vertices[i] += center;
4675  vertices[i + N + 1] += center;
4676  }
4677 
4678 
4679  std::vector<CellData<2>> cells(N, CellData<2>());
4680 
4681  for (unsigned int i = 0; i < N; ++i)
4682  {
4683  cells[i].vertices[0] = i;
4684  cells[i].vertices[1] = (i + 1) % (N + 1);
4685  cells[i].vertices[2] = N + 1 + i;
4686  cells[i].vertices[3] = N + 1 + ((i + 1) % (N + 1));
4687 
4688  cells[i].material_id = 0;
4689  }
4690 
4691  tria.create_triangulation(vertices, cells, SubCellData());
4692 
4693  if (colorize)
4694  {
4695  Triangulation<2>::cell_iterator cell = tria.begin();
4696  for (; cell != tria.end(); ++cell)
4697  {
4698  cell->face(2)->set_boundary_id(1);
4699  }
4700  tria.begin()->face(0)->set_boundary_id(3);
4701 
4702  tria.last()->face(1)->set_boundary_id(2);
4703  }
4704 
4705  tria.set_all_manifold_ids(0);
4706  tria.set_manifold(0, SphericalManifold<2>(center));
4707  }
4708 
4709 
4710 
4711  // Implementation for 3D only
4712  template <>
4713  void hyper_cube_slit(Triangulation<3> &tria,
4714  const double left,
4715  const double right,
4716  const bool colorize)
4717  {
4718  const double rl2 = (right + left) / 2;
4719  const double len = (right - left) / 2.;
4720 
4721  const Point<3> vertices[20] = {
4722  Point<3>(left, left, -len / 2.), Point<3>(rl2, left, -len / 2.),
4723  Point<3>(rl2, rl2, -len / 2.), Point<3>(left, rl2, -len / 2.),
4724  Point<3>(right, left, -len / 2.), Point<3>(right, rl2, -len / 2.),
4725  Point<3>(rl2, right, -len / 2.), Point<3>(left, right, -len / 2.),
4726  Point<3>(right, right, -len / 2.), Point<3>(rl2, left, -len / 2.),
4727  Point<3>(left, left, len / 2.), Point<3>(rl2, left, len / 2.),
4728  Point<3>(rl2, rl2, len / 2.), Point<3>(left, rl2, len / 2.),
4729  Point<3>(right, left, len / 2.), Point<3>(right, rl2, len / 2.),
4730  Point<3>(rl2, right, len / 2.), Point<3>(left, right, len / 2.),
4731  Point<3>(right, right, len / 2.), Point<3>(rl2, left, len / 2.)};
4732  const int cell_vertices[4][8] = {{0, 1, 3, 2, 10, 11, 13, 12},
4733  {9, 4, 2, 5, 19, 14, 12, 15},
4734  {3, 2, 7, 6, 13, 12, 17, 16},
4735  {2, 5, 6, 8, 12, 15, 16, 18}};
4736  std::vector<CellData<3>> cells(4, CellData<3>());
4737  for (unsigned int i = 0; i < 4; ++i)
4738  {
4739  for (unsigned int j = 0; j < 8; ++j)
4740  cells[i].vertices[j] = cell_vertices[i][j];
4741  cells[i].material_id = 0;
4742  }
4743  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
4744  std::end(vertices)),
4745  cells,
4746  SubCellData()); // no boundary information
4747 
4748  if (colorize)
4749  {
4750  Triangulation<3>::cell_iterator cell = tria.begin();
4751  cell->face(1)->set_boundary_id(1);
4752  ++cell;
4753  cell->face(0)->set_boundary_id(2);
4754  }
4755  }
4756 
4757 
4758 
4759  // Implementation for 3D only
4760  template <>
4762  const double left,
4763  const double right,
4764  const double thickness,
4765  const bool colorize)
4766  {
4767  Assert(left < right,
4768  ExcMessage("Invalid left-to-right bounds of enclosed hypercube"));
4769 
4770  std::vector<Point<3>> vertices(64);
4771  double coords[4];
4772  coords[0] = left - thickness;
4773  coords[1] = left;
4774  coords[2] = right;
4775  coords[3] = right + thickness;
4776 
4777  unsigned int k = 0;
4778  for (const double z : coords)
4779  for (const double y : coords)
4780  for (const double x : coords)
4781  vertices[k++] = Point<3>(x, y, z);
4782 
4783  const types::material_id materials[27] = {21, 20, 22, 17, 16, 18, 25,
4784  24, 26, 5, 4, 6, 1, 0,
4785  2, 9, 8, 10, 37, 36, 38,
4786  33, 32, 34, 41, 40, 42};
4787 
4788  std::vector<CellData<3>> cells(27);
4789  k = 0;
4790  for (unsigned int z = 0; z < 3; ++z)
4791  for (unsigned int y = 0; y < 3; ++y)
4792  for (unsigned int x = 0; x < 3; ++x)
4793  {
4794  cells[k].vertices[0] = x + 4 * y + 16 * z;
4795  cells[k].vertices[1] = x + 4 * y + 16 * z + 1;
4796  cells[k].vertices[2] = x + 4 * y + 16 * z + 4;
4797  cells[k].vertices[3] = x + 4 * y + 16 * z + 5;
4798  cells[k].vertices[4] = x + 4 * y + 16 * z + 16;
4799  cells[k].vertices[5] = x + 4 * y + 16 * z + 17;
4800  cells[k].vertices[6] = x + 4 * y + 16 * z + 20;
4801  cells[k].vertices[7] = x + 4 * y + 16 * z + 21;
4802  if (colorize)
4803  cells[k].material_id = materials[k];
4804  ++k;
4805  }
4806  tria.create_triangulation(vertices,
4807  cells,
4808  SubCellData()); // no boundary information
4809  }
4810 
4811 
4812 
4813  template <>
4814  void truncated_cone(Triangulation<3> &triangulation,
4815  const double radius_0,
4816  const double radius_1,
4817  const double half_length)
4818  {
4819  Assert(triangulation.n_cells() == 0,
4820  ExcMessage("The output triangulation object needs to be empty."));
4821  Assert(0 < radius_0, ExcMessage("The radii must be positive."));
4822  Assert(0 < radius_1, ExcMessage("The radii must be positive."));
4823  Assert(0 < half_length, ExcMessage("The half length must be positive."));
4824 
4825  const auto n_slices = 1 + static_cast<unsigned int>(std::ceil(
4826  half_length / std::max(radius_0, radius_1)));
4827 
4828  Triangulation<2> triangulation_2;
4829  GridGenerator::hyper_ball(triangulation_2, Point<2>(), radius_0);
4830  GridGenerator::extrude_triangulation(triangulation_2,
4831  n_slices,
4832  2 * half_length,
4833  triangulation);
4834  GridTools::rotate(numbers::PI / 2, 1, triangulation);
4835  GridTools::shift(Tensor<1, 3>({-half_length, 0.0, 0.0}), triangulation);
4836  // At this point we have a cylinder. Multiply the y and z coordinates by a
4837  // factor that scales (with x) linearly between radius_0 and radius_1 to fix
4838  // the circle radii and interior points:
4839  auto shift_radii = [=](const Point<3> &p) {
4840  const double slope = (radius_1 / radius_0 - 1.0) / (2.0 * half_length);
4841  const double factor = slope * (p[0] - -half_length) + 1.0;
4842  return Point<3>(p[0], factor * p[1], factor * p[2]);
4843  };
4844  GridTools::transform(shift_radii, triangulation);
4845 
4846  // Set boundary ids at -half_length to 1 and at half_length to 2. Set the
4847  // manifold id on hull faces (i.e., faces not on either end) to 0.
4848  for (const auto &face : triangulation.active_face_iterators())
4849  if (face->at_boundary())
4850  {
4851  if (std::abs(face->center()[0] - -half_length) < 1e-8 * half_length)
4852  face->set_boundary_id(1);
4853  else if (std::abs(face->center()[0] - half_length) <
4854  1e-8 * half_length)
4855  face->set_boundary_id(2);
4856  else
4857  face->set_all_manifold_ids(0);
4858  }
4859 
4860  triangulation.set_manifold(0, CylindricalManifold<3>());
4861  }
4862 
4863 
4864  // Implementation for 3D only
4865  template <>
4866  void hyper_L(Triangulation<3> &tria,
4867  const double a,
4868  const double b,
4869  const bool colorize)
4870  {
4871  // we slice out the top back right
4872  // part of the cube
4873  const Point<3> vertices[26] = {
4874  // front face of the big cube
4875  Point<3>(a, a, a),
4876  Point<3>((a + b) / 2, a, a),
4877  Point<3>(b, a, a),
4878  Point<3>(a, a, (a + b) / 2),
4879  Point<3>((a + b) / 2, a, (a + b) / 2),
4880  Point<3>(b, a, (a + b) / 2),
4881  Point<3>(a, a, b),
4882  Point<3>((a + b) / 2, a, b),
4883  Point<3>(b, a, b),
4884  // middle face of the big cube
4885  Point<3>(a, (a + b) / 2, a),
4886  Point<3>((a + b) / 2, (a + b) / 2, a),
4887  Point<3>(b, (a + b) / 2, a),
4888  Point<3>(a, (a + b) / 2, (a + b) / 2),
4889  Point<3>((a + b) / 2, (a + b) / 2, (a + b) / 2),
4890  Point<3>(b, (a + b) / 2, (a + b) / 2),
4891  Point<3>(a, (a + b) / 2, b),
4892  Point<3>((a + b) / 2, (a + b) / 2, b),
4893  Point<3>(b, (a + b) / 2, b),
4894  // back face of the big cube
4895  // last (top right) point is missing
4896  Point<3>(a, b, a),
4897  Point<3>((a + b) / 2, b, a),
4898  Point<3>(b, b, a),
4899  Point<3>(a, b, (a + b) / 2),
4900  Point<3>((a + b) / 2, b, (a + b) / 2),
4901  Point<3>(b, b, (a + b) / 2),
4902  Point<3>(a, b, b),
4903  Point<3>((a + b) / 2, b, b)};
4904  const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
4905  {1, 2, 10, 11, 4, 5, 13, 14},
4906  {3, 4, 12, 13, 6, 7, 15, 16},
4907  {4, 5, 13, 14, 7, 8, 16, 17},
4908  {9, 10, 18, 19, 12, 13, 21, 22},
4909  {10, 11, 19, 20, 13, 14, 22, 23},
4910  {12, 13, 21, 22, 15, 16, 24, 25}};
4911 
4912  std::vector<CellData<3>> cells(7, CellData<3>());
4913 
4914  for (unsigned int i = 0; i < 7; ++i)
4915  {
4916  for (unsigned int j = 0; j < 8; ++j)
4917  cells[i].vertices[j] = cell_vertices[i][j];
4918  cells[i].material_id = 0;
4919  }
4920 
4921  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
4922  std::end(vertices)),
4923  cells,
4924  SubCellData()); // no boundary information
4925 
4926  if (colorize)
4927  {
4928  Assert(false, ExcNotImplemented());
4929  }
4930  }
4931 
4932 
4933 
4934  // Implementation for 3D only
4935  template <>
4936  void hyper_ball(Triangulation<3> &tria,
4937  const Point<3> & p,
4938  const double radius,
4939  const bool internal_manifold)
4940  {
4941  const double a =
4942  1. / (1 + std::sqrt(3.0)); // equilibrate cell sizes at transition
4943  // from the inner part to the radial
4944  // cells
4945  const unsigned int n_vertices = 16;
4946  const Point<3> vertices[n_vertices] = {
4947  // first the vertices of the inner
4948  // cell
4949  p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0) * a),
4950  p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0) * a),
4951  p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0) * a),
4952  p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0) * a),
4953  p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0) * a),
4954  p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0) * a),
4955  p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0) * a),
4956  p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0) * a),
4957  // now the eight vertices at
4958  // the outer sphere
4959  p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0)),
4960  p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0)),
4961  p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0)),
4962  p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0)),
4963  p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0)),
4964  p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0)),
4965  p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0)),
4966  p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0)),
4967  };
4968 
4969  // one needs to draw the seven cubes to
4970  // understand what's going on here
4971  const unsigned int n_cells = 7;
4972  const int cell_vertices[n_cells][8] = {
4973  {0, 1, 4, 5, 3, 2, 7, 6}, // center
4974  {8, 9, 12, 13, 0, 1, 4, 5}, // bottom
4975  {9, 13, 1, 5, 10, 14, 2, 6}, // right
4976  {11, 10, 3, 2, 15, 14, 7, 6}, // top
4977  {8, 0, 12, 4, 11, 3, 15, 7}, // left
4978  {8, 9, 0, 1, 11, 10, 3, 2}, // front
4979  {12, 4, 13, 5, 15, 7, 14, 6}}; // back
4980 
4981  std::vector<CellData<3>> cells(n_cells, CellData<3>());
4982 
4983  for (unsigned int i = 0; i < n_cells; ++i)
4984  {
4985  for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j)
4986  cells[i].vertices[j] = cell_vertices[i][j];
4987  cells[i].material_id = 0;
4988  cells[i].manifold_id = i == 0 ? numbers::flat_manifold_id : 1;
4989  }
4990 
4991  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
4992  std::end(vertices)),
4993  cells,
4994  SubCellData()); // no boundary information
4996  tria.set_manifold(0, SphericalManifold<3>(p));
4997  if (internal_manifold)
4998  tria.set_manifold(1, SphericalManifold<3>(p));
4999  }
5000 
5001 
5002 
5003  template <int spacedim>
5005  const Point<spacedim> & p,
5006  const double radius)
5007  {
5008  Triangulation<spacedim> volume_mesh;
5009  GridGenerator::hyper_ball(volume_mesh, p, radius);
5010  std::set<types::boundary_id> boundary_ids;
5011  boundary_ids.insert(0);
5012  GridGenerator::extract_boundary_mesh(volume_mesh, tria, boundary_ids);
5013  tria.set_all_manifold_ids(0);
5015  }
5016 
5017 
5018 
5019  // Implementation for 3D only
5020  template <>
5021  void cylinder(Triangulation<3> &tria,
5022  const double radius,
5023  const double half_length)
5024  {
5025  // Copy the base from hyper_ball<3>
5026  // and transform it to yz
5027  const double d = radius / std::sqrt(2.0);
5028  const double a = d / (1 + std::sqrt(2.0));
5029  Point<3> vertices[24] = {
5030  Point<3>(-d, -half_length, -d),
5031  Point<3>(d, -half_length, -d),
5032  Point<3>(-a, -half_length, -a),
5033  Point<3>(a, -half_length, -a),
5034  Point<3>(-a, -half_length, a),
5035  Point<3>(a, -half_length, a),
5036  Point<3>(-d, -half_length, d),
5037  Point<3>(d, -half_length, d),
5038  Point<3>(-d, 0, -d),
5039  Point<3>(d, 0, -d),
5040  Point<3>(-a, 0, -a),
5041  Point<3>(a, 0, -a),
5042  Point<3>(-a, 0, a),
5043  Point<3>(a, 0, a),
5044  Point<3>(-d, 0, d),
5045  Point<3>(d, 0, d),
5046  Point<3>(-d, half_length, -d),
5047  Point<3>(d, half_length, -d),
5048  Point<3>(-a, half_length, -a),
5049  Point<3>(a, half_length, -a),
5050  Point<3>(-a, half_length, a),
5051  Point<3>(a, half_length, a),
5052  Point<3>(-d, half_length, d),
5053  Point<3>(d, half_length, d),
5054  };
5055  // Turn cylinder such that y->x
5056  for (auto &vertex : vertices)
5057  {
5058  const double h = vertex(1);
5059  vertex(1) = -vertex(0);
5060  vertex(0) = h;
5061  }
5062 
5063  int cell_vertices[10][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
5064  {0, 2, 8, 10, 6, 4, 14, 12},
5065  {2, 3, 10, 11, 4, 5, 12, 13},
5066  {1, 7, 9, 15, 3, 5, 11, 13},
5067  {6, 4, 14, 12, 7, 5, 15, 13}};
5068  for (unsigned int i = 0; i < 5; ++i)
5069  for (unsigned int j = 0; j < 8; ++j)
5070  cell_vertices[i + 5][j] = cell_vertices[i][j] + 8;
5071 
5072  std::vector<CellData<3>> cells(10, CellData<3>());
5073 
5074  for (unsigned int i = 0; i < 10; ++i)
5075  {
5076  for (unsigned int j = 0; j < 8; ++j)
5077  cells[i].vertices[j] = cell_vertices[i][j];
5078  cells[i].material_id = 0;
5079  }
5080 
5081  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
5082  std::end(vertices)),
5083  cells,
5084  SubCellData()); // no boundary information
5085 
5086  // set boundary indicators for the
5087  // faces at the ends to 1 and 2,
5088  // respectively. note that we also
5089  // have to deal with those lines
5090  // that are purely in the interior
5091  // of the ends. we determine whether
5092  // an edge is purely in the
5093  // interior if one of its vertices
5094  // is at coordinates '+-a' as set
5095  // above
5096  Triangulation<3>::cell_iterator cell = tria.begin();
5097  Triangulation<3>::cell_iterator end = tria.end();
5098 
5100 
5101  for (; cell != end; ++cell)
5102  for (unsigned int i : GeometryInfo<3>::face_indices())
5103  if (cell->at_boundary(i))
5104  {
5105  if (cell->face(i)->center()(0) > half_length - 1.e-5)
5106  {
5107  cell->face(i)->set_boundary_id(2);
5108  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
5109 
5110  for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5111  ++e)
5112  if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5113  (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5114  (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5115  (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5116  {
5117  cell->face(i)->line(e)->set_boundary_id(2);
5118  cell->face(i)->line(e)->set_manifold_id(
5120  }
5121  }
5122  else if (cell->face(i)->center()(0) < -half_length + 1.e-5)
5123  {
5124  cell->face(i)->set_boundary_id(1);
5125  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
5126 
5127  for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5128  ++e)
5129  if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5130  (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5131  (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5132  (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5133  {
5134  cell->face(i)->line(e)->set_boundary_id(1);
5135  cell->face(i)->line(e)->set_manifold_id(
5137  }
5138  }
5139  }
5141  }
5142 
5143 
5144  template <>
5146  const Point<3> & center,
5147  const double radius)
5148  {
5149  const unsigned int dim = 3;
5150 
5151  // equilibrate cell sizes at
5152  // transition from the inner part
5153  // to the radial cells
5154  const Point<dim> vertices[15] = {
5155  center + Point<dim>(0, 0, 0) * radius,
5156  center + Point<dim>(+1, 0, 0) * radius,
5157  center + Point<dim>(+1, 0, 0) * (radius / 2.),
5158  center + Point<dim>(0, +1, 0) * (radius / 2.),
5159  center + Point<dim>(+1, +1, 0) * (radius / (2 * std::sqrt(2.0))),
5160  center + Point<dim>(0, +1, 0) * radius,
5161  center + Point<dim>(+1, +1, 0) * (radius / std::sqrt(2.0)),
5162  center + Point<dim>(0, 0, 1) * radius / 2.,
5163  center + Point<dim>(+1, 0, 1) * radius / std::sqrt(2.0),
5164  center + Point<dim>(+1, 0, 1) * (radius / (2 * std::sqrt(2.0))),
5165  center + Point<dim>(0, +1, 1) * (radius / (2 * std::sqrt(2.0))),
5166  center + Point<dim>(+1, +1, 1) * (radius / (2 * std::sqrt(3.0))),
5167  center + Point<dim>(0, +1, 1) * radius / std::sqrt(2.0),
5168  center + Point<dim>(+1, +1, 1) * (radius / (std::sqrt(3.0))),
5169  center + Point<dim>(0, 0, 1) * radius};
5170  const int cell_vertices[4][8] = {{0, 2, 3, 4, 7, 9, 10, 11},
5171  {1, 6, 2, 4, 8, 13, 9, 11},
5172  {5, 3, 6, 4, 12, 10, 13, 11},
5173  {7, 9, 10, 11, 14, 8, 12, 13}};
5174 
5175  std::vector<CellData<dim>> cells(4, CellData<dim>());
5176 
5177  for (unsigned int i = 0; i < 4; ++i)
5178  {
5179  for (unsigned int j = 0; j < 8; ++j)
5180  cells[i].vertices[j] = cell_vertices[i][j];
5181  cells[i].material_id = 0;
5182  }
5183 
5184  tria.create_triangulation(std::vector<Point<dim>>(std::begin(vertices),
5185  std::end(vertices)),
5186  cells,
5187  SubCellData()); // no boundary information
5188 
5191 
5193  while (cell != end)
5194  {
5195  for (unsigned int i : GeometryInfo<dim>::face_indices())
5196  {
5197  if (cell->face(i)->boundary_id() ==
5199  continue;
5200 
5201  // If x,y or z is zero, then this is part of the plane
5202  if (cell->face(i)->center()(0) < center(0) + 1.e-5 * radius ||
5203  cell->face(i)->center()(1) < center(1) + 1.e-5 * radius ||
5204  cell->face(i)->center()(2) < center(2) + 1.e-5 * radius)
5205  {
5206  cell->face(i)->set_boundary_id(1);
5207  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
5208  // also set the boundary indicators of the bounding lines,
5209  // unless both vertices are on the perimeter
5210  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5211  ++j)
5212  {
5213  const Point<3> line_vertices[2] = {
5214  cell->face(i)->line(j)->vertex(0),
5215  cell->face(i)->line(j)->vertex(1)};
5216  if ((std::fabs(line_vertices[0].distance(center) - radius) >
5217  1e-5 * radius) ||
5218  (std::fabs(line_vertices[1].distance(center) - radius) >
5219  1e-5 * radius))
5220  {
5221  cell->face(i)->line(j)->set_boundary_id(1);
5222  cell->face(i)->line(j)->set_manifold_id(
5224  }
5225  }
5226  }
5227  }
5228  ++cell;
5229  }
5230  tria.set_manifold(0, SphericalManifold<3>(center));
5231  }
5232 
5233 
5234  // Implementation for 3D only
5235  template <>
5236  void half_hyper_ball(Triangulation<3> &tria,
5237  const Point<3> & center,
5238  const double radius)
5239  {
5240  // These are for the two lower squares
5241  const double d = radius / std::sqrt(2.0);
5242  const double a = d / (1 + std::sqrt(2.0));
5243  // These are for the two upper square
5244  const double b = a / 2.0;
5245  const double c = d / 2.0;
5246  // And so are these
5247  const double hb = radius * std::sqrt(3.0) / 4.0;
5248  const double hc = radius * std::sqrt(3.0) / 2.0;
5249 
5250  Point<3> vertices[16] = {
5251  center + Point<3>(0, d, -d),
5252  center + Point<3>(0, -d, -d),
5253  center + Point<3>(0, a, -a),
5254  center + Point<3>(0, -a, -a),
5255  center + Point<3>(0, a, a),
5256  center + Point<3>(0, -a, a),
5257  center + Point<3>(0, d, d),
5258  center + Point<3>(0, -d, d),
5259 
5260  center + Point<3>(hc, c, -c),
5261  center + Point<3>(hc, -c, -c),
5262  center + Point<3>(hb, b, -b),
5263  center + Point<3>(hb, -b, -b),
5264  center + Point<3>(hb, b, b),
5265  center + Point<3>(hb, -b, b),
5266  center + Point<3>(hc, c, c),
5267  center + Point<3>(hc, -c, c),
5268  };
5269 
5270  int cell_vertices[6][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
5271  {0, 2, 8, 10, 6, 4, 14, 12},
5272  {2, 3, 10, 11, 4, 5, 12, 13},
5273  {1, 7, 9, 15, 3, 5, 11, 13},
5274  {6, 4, 14, 12, 7, 5, 15, 13},
5275  {8, 10, 9, 11, 14, 12, 15, 13}};
5276 
5277  std::vector<CellData<3>> cells(6, CellData<3>());
5278 
5279  for (unsigned int i = 0; i < 6; ++i)
5280  {
5281  for (unsigned int j = 0; j < 8; ++j)
5282  cells[i].vertices[j] = cell_vertices[i][j];
5283  cells[i].material_id = 0;
5284  }
5285 
5286  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
5287  std::end(vertices)),
5288  cells,
5289  SubCellData()); // no boundary information
5290 
5291  Triangulation<3>::cell_iterator cell = tria.begin();
5292  Triangulation<3>::cell_iterator end = tria.end();
5293 
5295 
5296  // go over all faces. for the ones on the flat face, set boundary
5297  // indicator for face and edges to one; the rest will remain at
5298  // zero but we have to pay attention to those edges that are
5299  // at the perimeter of the flat face since they should not be
5300  // set to one
5301  while (cell != end)
5302  {
5303  for (unsigned int i : GeometryInfo<3>::face_indices())
5304  {
5305  if (!cell->at_boundary(i))
5306  continue;
5307 
5308  // If the center is on the plane x=0, this is a planar element. set
5309  // its boundary indicator. also set the boundary indicators of the
5310  // bounding faces unless both vertices are on the perimeter
5311  if (cell->face(i)->center()(0) < center(0) + 1.e-5 * radius)
5312  {
5313  cell->face(i)->set_boundary_id(1);
5314  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
5315  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5316  ++j)
5317  {
5318  const Point<3> line_vertices[2] = {
5319  cell->face(i)->line(j)->vertex(0),
5320  cell->face(i)->line(j)->vertex(1)};
5321  if ((std::fabs(line_vertices[0].distance(center) - radius) >
5322  1e-5 * radius) ||
5323  (std::fabs(line_vertices[1].distance(center) - radius) >
5324  1e-5 * radius))
5325  {
5326  cell->face(i)->line(j)->set_boundary_id(1);
5327  cell->face(i)->line(j)->set_manifold_id(
5329  }
5330  }
5331  }
5332  }
5333  ++cell;
5334  }
5335  tria.set_manifold(0, SphericalManifold<3>(center));
5336  }
5337 
5338 
5339  template <>
5340  void hyper_shell(Triangulation<3> & tria,
5341  const Point<3> & p,
5342  const double inner_radius,
5343  const double outer_radius,
5344  const unsigned int n_cells,
5345  const bool colorize)
5346  {
5347  Assert((inner_radius > 0) && (inner_radius < outer_radius),
5348  ExcInvalidRadii());
5349 
5350  const unsigned int n = (n_cells == 0) ? 6 : n_cells;
5351 
5352  const double irad = inner_radius / std::sqrt(3.0);
5353  const double orad = outer_radius / std::sqrt(3.0);
5354  std::vector<Point<3>> vertices;
5355  std::vector<CellData<3>> cells;
5356 
5357  // Corner points of the cube [-1,1]^3
5358  static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1}, //
5359  {+1, -1, -1}, //
5360  {-1, +1, -1}, //
5361  {+1, +1, -1}, //
5362  {-1, -1, +1}, //
5363  {+1, -1, +1}, //
5364  {-1, +1, +1}, //
5365  {+1, +1, +1}}};
5366 
5367  // Start with the shell bounded by
5368  // two nested cubes
5369  if (n == 6)
5370  {
5371  for (unsigned int i = 0; i < 8; ++i)
5372  vertices.push_back(p + hexahedron[i] * irad);
5373  for (unsigned int i = 0; i < 8; ++i)
5374  vertices.push_back(p + hexahedron[i] * orad);
5375 
5376  const unsigned int n_cells = 6;
5377  const int cell_vertices[n_cells][8] = {
5378  {8, 9, 10, 11, 0, 1, 2, 3}, // bottom
5379  {9, 11, 1, 3, 13, 15, 5, 7}, // right
5380  {12, 13, 4, 5, 14, 15, 6, 7}, // top
5381  {8, 0, 10, 2, 12, 4, 14, 6}, // left
5382  {8, 9, 0, 1, 12, 13, 4, 5}, // front
5383  {10, 2, 11, 3, 14, 6, 15, 7}}; // back
5384 
5385  cells.resize(n_cells, CellData<3>());
5386 
5387  for (unsigned int i = 0; i < n_cells; ++i)
5388  {
5389  for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell;
5390  ++j)
5391  cells[i].vertices[j] = cell_vertices[i][j];
5392  cells[i].material_id = 0;
5393  }
5394 
5395  tria.create_triangulation(vertices, cells, SubCellData());
5396  }
5397  // A more regular subdivision can
5398  // be obtained by two nested
5399  // rhombic dodecahedra
5400 
5401  else if (n == 12)
5402  {
5403  // Octahedron inscribed in the cube [-1,1]^3
5404  static const std::array<Point<3>, 6> octahedron = {{{-1, 0, 0}, //
5405  {1, 0, 0}, //
5406  {0, -1, 0}, //
5407  {0, 1, 0}, //
5408  {0, 0, -1}, //
5409  {0, 0, 1}}};
5410 
5411  for (unsigned int i = 0; i < 8; ++i)
5412  vertices.push_back(p + hexahedron[i] * irad);
5413  for (unsigned int i = 0; i < 6; ++i)
5414  vertices.push_back(p + octahedron[i] * inner_radius);
5415  for (unsigned int i = 0; i < 8; ++i)
5416  vertices.push_back(p + hexahedron[i] * orad);
5417  for (unsigned int i = 0; i < 6; ++i)
5418  vertices.push_back(p + octahedron[i] * outer_radius);
5419 
5420  const unsigned int n_cells = 12;
5421  const unsigned int rhombi[n_cells][4] = {{10, 4, 0, 8},
5422  {4, 13, 8, 6},
5423  {10, 5, 4, 13},
5424  {1, 9, 10, 5},
5425  {9, 7, 5, 13},
5426  {7, 11, 13, 6},
5427  {9, 3, 7, 11},
5428  {1, 12, 9, 3},
5429  {12, 2, 3, 11},
5430  {2, 8, 11, 6},
5431  {12, 0, 2, 8},
5432  {1, 10, 12, 0}};
5433 
5434  cells.resize(n_cells, CellData<3>());
5435 
5436  for (unsigned int i = 0; i < n_cells; ++i)
5437  {
5438  for (unsigned int j = 0; j < 4; ++j)
5439  {
5440  cells[i].vertices[j] = rhombi[i][j];
5441  cells[i].vertices[j + 4] = rhombi[i][j] + 14;
5442  }
5443  cells[i].material_id = 0;
5444  }
5445 
5446  tria.create_triangulation(vertices, cells, SubCellData());
5447  }
5448  else if (n == 96)
5449  {
5450  // create a triangulation based on the 12-cell version. This function
5451  // was needed before SphericalManifold was written: it manually
5452  // adjusted the interior vertices to lie along concentric
5453  // spheres. Nowadays we can just refine globally:
5454  Triangulation<3> tmp;
5455  hyper_shell(tmp, p, inner_radius, outer_radius, 12);
5456  tmp.refine_global(1);
5457 
5458  // now copy the resulting level 1 cells into the new triangulation,
5459  cells.resize(tmp.n_active_cells(), CellData<3>());
5460  for (const auto &cell : tmp.active_cell_iterators())
5461  {
5462  const unsigned int cell_index = cell->active_cell_index();
5463  for (unsigned int v = 0; v < GeometryInfo<3>::vertices_per_cell;
5464  ++v)
5465  cells[cell_index].vertices[v] = cell->vertex_index(v);
5466  cells[cell_index].material_id = 0;
5467  }
5468 
5469  tria.create_triangulation(tmp.get_vertices(), cells, SubCellData());
5470  }
5471  else
5472  {
5473  Assert(false, ExcMessage("Invalid number of coarse mesh cells."));
5474  }
5475 
5476  if (colorize)
5477  colorize_hyper_shell(tria, p, inner_radius, outer_radius);
5478  tria.set_all_manifold_ids(0);
5479  tria.set_manifold(0, SphericalManifold<3>(p));
5480  }
5481 
5482 
5483 
5484  // Implementation for 3D only
5485  template <>
5487  const Point<3> & center,
5488  const double inner_radius,
5489  const double outer_radius,
5490  const unsigned int /*n_cells*/,
5491  const bool colorize)
5492  {
5493  Assert((inner_radius > 0) && (inner_radius < outer_radius),
5494  ExcInvalidRadii());
5495 
5496  // These are for the two lower squares
5497  const double d = outer_radius / std::sqrt(2.0);
5498  const double a = inner_radius / std::sqrt(2.0);
5499  // These are for the two upper square
5500  const double b = a / 2.0;
5501  const double c = d / 2.0;
5502  // And so are these
5503  const double hb = inner_radius * std::sqrt(3.0) / 2.0;
5504  const double hc = outer_radius * std::sqrt(3.0) / 2.0;
5505 
5506  Point<3> vertices[16] = {
5507  center + Point<3>(0, d, -d),
5508  center + Point<3>(0, -d, -d),
5509  center + Point<3>(0, a, -a),
5510  center + Point<3>(0, -a, -a),
5511  center + Point<3>(0, a, a),
5512  center + Point<3>(0, -a, a),
5513  center + Point<3>(0, d, d),
5514  center + Point<3>(0, -d, d),
5515 
5516  center + Point<3>(hc, c, -c),
5517  center + Point<3>(hc, -c, -c),
5518  center + Point<3>(hb, b, -b),
5519  center + Point<3>(hb, -b, -b),
5520  center + Point<3>(hb, b, b),
5521  center + Point<3>(hb, -b, b),
5522  center + Point<3>(hc, c, c),
5523  center + Point<3>(hc, -c, c),
5524  };
5525 
5526  int cell_vertices[5][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
5527  {0, 2, 8, 10, 6, 4, 14, 12},
5528  {1, 7, 9, 15, 3, 5, 11, 13},
5529  {6, 4, 14, 12, 7, 5, 15, 13},
5530  {8, 10, 9, 11, 14, 12, 15, 13}};
5531 
5532  std::vector<CellData<3>> cells(5, CellData<3>());
5533 
5534  for (unsigned int i = 0; i < 5; ++i)
5535  {
5536  for (unsigned int j = 0; j < 8; ++j)
5537  cells[i].vertices[j] = cell_vertices[i][j];
5538  cells[i].material_id = 0;
5539  }
5540 
5541  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
5542  std::end(vertices)),
5543  cells,
5544  SubCellData()); // no boundary information
5545 
5546  if (colorize)
5547  {
5548  // We want to use a standard boundary description where
5549  // the boundary is not curved. Hence set boundary id 2 to
5550  // to all faces in a first step.
5551  Triangulation<3>::cell_iterator cell = tria.begin();
5552  for (; cell != tria.end(); ++cell)
5553  for (unsigned int i : GeometryInfo<3>::face_indices())
5554  if (cell->at_boundary(i))
5555  cell->face(i)->set_all_boundary_ids(2);
5556 
5557  // Next look for the curved boundaries. If the x value of the
5558  // center of the face is not equal to center(0), we're on a curved
5559  // boundary. Then decide whether the center is nearer to the inner
5560  // or outer boundary to set the correct boundary id.
5561  for (cell = tria.begin(); cell != tria.end(); ++cell)
5562  for (unsigned int i : GeometryInfo<3>::face_indices())
5563  if (cell->at_boundary(i))
5564  {
5565  const Triangulation<3>::face_iterator face = cell->face(i);
5566 
5567  const Point<3> face_center(face->center());
5568  if (std::abs(face_center(0) - center(0)) >
5569  1.e-6 * face_center.norm())
5570  {
5571  if (std::abs((face_center - center).norm() - inner_radius) <
5572  std::abs((face_center - center).norm() - outer_radius))
5573  face->set_all_boundary_ids(0);
5574  else
5575  face->set_all_boundary_ids(1);
5576  }
5577  }
5578  }
5579  tria.set_all_manifold_ids(0);
5580  tria.set_manifold(0, SphericalManifold<3>(center));
5581  }
5582 
5583 
5584  // Implementation for 3D only
5585  template <>
5587  const Point<3> & center,
5588  const double inner_radius,
5589  const double outer_radius,
5590  const unsigned int n,
5591  const bool colorize)
5592  {
5593  Assert((inner_radius > 0) && (inner_radius < outer_radius),
5594  ExcInvalidRadii());
5595  if (n == 0 || n == 3)
5596  {
5597  const double a = inner_radius * std::sqrt(2.0) / 2e0;
5598  const double b = outer_radius * std::sqrt(2.0) / 2e0;
5599  const double c = a * std::sqrt(3.0) / 2e0;
5600  const double d = b * std::sqrt(3.0) / 2e0;
5601  const double e = outer_radius / 2e0;
5602  const double h = inner_radius / 2e0;
5603 
5604  std::vector<Point<3>> vertices;
5605 
5606  vertices.push_back(center + Point<3>(0, inner_radius, 0)); // 0
5607  vertices.push_back(center + Point<3>(a, a, 0)); // 1
5608  vertices.push_back(center + Point<3>(b, b, 0)); // 2
5609  vertices.push_back(center + Point<3>(0, outer_radius, 0)); // 3
5610  vertices.push_back(center + Point<3>(0, a, a)); // 4
5611  vertices.push_back(center + Point<3>(c, c, h)); // 5
5612  vertices.push_back(center + Point<3>(d, d, e)); // 6
5613  vertices.push_back(center + Point<3>(0, b, b)); // 7
5614  vertices.push_back(center + Point<3>(inner_radius, 0, 0)); // 8
5615  vertices.push_back(center + Point<3>(outer_radius, 0, 0)); // 9
5616  vertices.push_back(center + Point<3>(a, 0, a)); // 10
5617  vertices.push_back(center + Point<3>(b, 0, b)); // 11
5618  vertices.push_back(center + Point<3>(0, 0, inner_radius)); // 12
5619  vertices.push_back(center + Point<3>(0, 0, outer_radius)); // 13
5620 
5621  const int cell_vertices[3][8] = {
5622  {0, 1, 3, 2, 4, 5, 7, 6},
5623  {1, 8, 2, 9, 5, 10, 6, 11},
5624  {4, 5, 7, 6, 12, 10, 13, 11},
5625  };
5626  std::vector<CellData<3>> cells(3);
5627 
5628  for (unsigned int i = 0; i < 3; ++i)
5629  {
5630  for (unsigned int j = 0; j < 8; ++j)
5631  cells[i].vertices[j] = cell_vertices[i][j];
5632  cells[i].material_id = 0;
5633  }
5634 
5635  tria.create_triangulation(vertices,
5636  cells,
5637  SubCellData()); // no boundary information
5638  }
5639  else
5640  {
5641  AssertThrow(false, ExcNotImplemented());
5642  }
5643 
5644  if (colorize)
5645  colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius);
5646 
5647  tria.set_all_manifold_ids(0);
5648  tria.set_manifold(0, SphericalManifold<3>(center));
5649  }
5650 
5651 
5652  // Implementation for 3D only
5653  template <>
5654  void cylinder_shell(Triangulation<3> & tria,
5655  const double length,
5656  const double inner_radius,
5657  const double outer_radius,
5658  const unsigned int n_radial_cells,
5659  const unsigned int n_axial_cells)
5660  {
5661  Assert((inner_radius > 0) && (inner_radius < outer_radius),
5662  ExcInvalidRadii());
5663 
5664  const double pi = numbers::PI;
5665 
5666  // determine the number of cells
5667  // for the grid. if not provided by
5668  // the user determine it such that
5669  // the length of each cell on the
5670  // median (in the middle between
5671  // the two circles) is equal to its
5672  // radial extent (which is the
5673  // difference between the two
5674  // radii)
5675  const unsigned int N_r =
5676  (n_radial_cells == 0 ? static_cast<unsigned int>(std::ceil(
5677  (2 * pi * (outer_radius + inner_radius) / 2) /
5678  (outer_radius - inner_radius))) :
5679  n_radial_cells);
5680  const unsigned int N_z =
5681  (n_axial_cells == 0 ?
5682  static_cast<unsigned int>(std::ceil(
5683  length / (2 * pi * (outer_radius + inner_radius) / 2 / N_r))) :
5684  n_axial_cells);
5685 
5686  // set up N vertices on the
5687  // outer and N vertices on
5688  // the inner circle. the
5689  // first N ones are on the
5690  // outer one, and all are
5691  // numbered counter-clockwise
5692  std::vector<Point<2>> vertices_2d(2 * N_r);
5693  for (unsigned int i = 0; i < N_r; ++i)
5694  {
5695  vertices_2d[i] =
5696  Point<2>(std::cos(2 * pi * i / N_r), std::sin(2 * pi * i / N_r)) *
5697  outer_radius;
5698  vertices_2d[i + N_r] = vertices_2d[i] * (inner_radius / outer_radius);
5699  }
5700 
5701  std::vector<Point<3>> vertices_3d;
5702  vertices_3d.reserve(2 * N_r * (N_z + 1));
5703  for (unsigned int j = 0; j <= N_z; ++j)
5704  for (unsigned int i = 0; i < 2 * N_r; ++i)
5705  {
5706  const Point<3> v(vertices_2d[i][0],
5707  vertices_2d[i][1],
5708  j * length / N_z);
5709  vertices_3d.push_back(v);
5710  }
5711 
5712  std::vector<CellData<3>> cells(N_r * N_z, CellData<3>());
5713 
5714  for (unsigned int j = 0; j < N_z; ++j)
5715  for (unsigned int i = 0; i < N_r; ++i)
5716  {
5717  cells[i + j * N_r].vertices[0] = i + (j + 1) * 2 * N_r;
5718  cells[i + j * N_r].vertices[1] = (i + 1) % N_r + (j + 1) * 2 * N_r;
5719  cells[i + j * N_r].vertices[2] = i + j * 2 * N_r;
5720  cells[i + j * N_r].vertices[3] = (i + 1) % N_r + j * 2 * N_r;
5721 
5722  cells[i + j * N_r].vertices[4] = N_r + i + (j + 1) * 2 * N_r;
5723  cells[i + j * N_r].vertices[5] =
5724  N_r + ((i + 1) % N_r) + (j + 1) * 2 * N_r;
5725  cells[i + j * N_r].vertices[6] = N_r + i + j * 2 * N_r;
5726  cells[i + j * N_r].vertices[7] = N_r + ((i + 1) % N_r) + j * 2 * N_r;
5727 
5728  cells[i + j * N_r].material_id = 0;
5729  }
5730 
5731  tria.create_triangulation(vertices_3d, cells, SubCellData());
5732  tria.set_all_manifold_ids(0);
5734  }
5735 
5736 
5737 
5738  template <int dim, int spacedim>
5739  void
5741  const std::initializer_list<const Triangulation<dim, spacedim> *const>
5742  & triangulations,
5744  const double duplicated_vertex_tolerance,
5745  const bool copy_manifold_ids)
5746  {
5747  std::vector<Point<spacedim>> vertices;
5748  std::vector<CellData<dim>> cells;
5749  SubCellData subcell_data;
5750 
5751  unsigned int n_accumulated_vertices = 0;
5752  for (const auto triangulation : triangulations)
5753  {
5754  Assert(triangulation->n_levels() == 1,
5755  ExcMessage("The input triangulations must be non-empty "
5756  "and must not be refined."));
5757 
5758  std::vector<Point<spacedim>> tria_vertices;
5759  std::vector<CellData<dim>> tria_cells;
5760  SubCellData tria_subcell_data;
5761  std::tie(tria_vertices, tria_cells, tria_subcell_data) =
5763 
5764  vertices.insert(vertices.end(),
5765  tria_vertices.begin(),
5766  tria_vertices.end());
5767  for (CellData<dim> &cell_data : tria_cells)
5768  {
5769  for (unsigned int &vertex_n : cell_data.vertices)
5770  vertex_n += n_accumulated_vertices;
5771  cells.push_back(cell_data);
5772  }
5773 
5774  // Skip copying lines with no manifold information.
5775  if (copy_manifold_ids)
5776  {
5777  for (CellData<1> &line_data : tria_subcell_data.boundary_lines)
5778  {
5779  if (line_data.manifold_id == numbers::flat_manifold_id)
5780  continue;
5781  for (unsigned int &vertex_n : line_data.vertices)
5782  vertex_n += n_accumulated_vertices;
5783  line_data.boundary_id =
5785  subcell_data.boundary_lines.push_back(line_data);
5786  }
5787 
5788  for (CellData<2> &quad_data : tria_subcell_data.boundary_quads)
5789  {
5790  if (quad_data.manifold_id == numbers::flat_manifold_id)
5791  continue;
5792  for (unsigned int &vertex_n : quad_data.vertices)
5793  vertex_n += n_accumulated_vertices;
5794  quad_data.boundary_id =
5796  subcell_data.boundary_quads.push_back(quad_data);
5797  }
5798  }
5799 
5800  n_accumulated_vertices += triangulation->n_vertices();
5801  }
5802 
5803  // throw out duplicated vertices
5804  std::vector<unsigned int> considered_vertices;
5806  cells,
5807  subcell_data,
5808  considered_vertices,
5809  duplicated_vertex_tolerance);
5810 
5811  // reorder the cells to ensure that they satisfy the convention for
5812  // edge and face directions
5814  result.clear();
5815  result.create_triangulation(vertices, cells, subcell_data);
5816  }
5817 
5818 
5819 
5820  template <int dim, int spacedim>
5821  void
5823  const Triangulation<dim, spacedim> &triangulation_2,
5825  const double duplicated_vertex_tolerance,
5826  const bool copy_manifold_ids)
5827  {
5828  // if either Triangulation is empty then merging is just a copy.
5829  if (triangulation_1.n_cells() == 0)
5830  {
5831  result.copy_triangulation(triangulation_2);
5832  return;
5833  }
5834  if (triangulation_2.n_cells() == 0)
5835  {
5836  result.copy_triangulation(triangulation_1);
5837  return;
5838  }
5839  merge_triangulations({&triangulation_1, &triangulation_2},
5840  result,
5841  duplicated_vertex_tolerance,
5842  copy_manifold_ids);
5843  }
5844 
5845 
5846 
5847  namespace
5848  {
5870  template <int structdim>
5871  void
5872  delete_duplicated_objects(std::vector<CellData<structdim>> &subcell_data)
5873  {
5874  static_assert(structdim == 1 || structdim == 2,
5875  "This function is only implemented for lines and "
5876  "quadrilaterals.");
5877  // start by making sure that all objects representing the same vertices
5878  // are numbered in the same way by canonicalizing the numberings. This
5879  // makes it possible to detect duplicates.
5880  for (CellData<structdim> &cell_data : subcell_data)
5881  {
5882  if (structdim == 1)
5883  std::sort(std::begin(cell_data.vertices),
5884  std::end(cell_data.vertices));
5885  else if (structdim == 2)
5886  {
5887  // rotate the vertex numbers so that the lowest one is first
5888  std::array<unsigned int, 4> renumbering;
5889  std::copy(std::begin(cell_data.vertices),
5890  std::end(cell_data.vertices),
5891  renumbering.begin());
5892 
5893  // convert to old style vertex numbering. This makes the
5894  // permutations easy since the valid configurations are
5895  //
5896  // 3 2 2 1 1 0 0 3
5897  // 0 1 3 0 2 3 1 2
5898  // (0123) (3012) (2310) (1230)
5899  //
5900  // rather than the lexical ordering which is harder to permute
5901  // by rotation.
5902  std::swap(renumbering[2], renumbering[3]);
5903  std::rotate(renumbering.begin(),
5904  std::min_element(renumbering.begin(),
5905  renumbering.end()),
5906  renumbering.end());
5907  // convert to new style
5908  std::swap(renumbering[2], renumbering[3]);
5909  // deal with cases where we might have
5910  //
5911  // 3 2 1 2
5912  // 0 1 0 3
5913  //
5914  // by forcing the second vertex (in lexical ordering) to be
5915  // smaller than the third
5916  if (renumbering[1] > renumbering[2])
5917  std::swap(renumbering[1], renumbering[2]);
5918  std::copy(renumbering.begin(),
5919  renumbering.end(),
5920  std::begin(cell_data.vertices));
5921  }
5922  }
5923 
5924  // Now that all cell objects have been canonicalized they can be sorted:
5925  auto compare = [](const CellData<structdim> &a,
5926  const CellData<structdim> &b) {
5927  return std::lexicographical_compare(std::begin(a.vertices),
5928  std::end(a.vertices),
5929  std::begin(b.vertices),
5930  std::end(b.vertices));
5931  };
5932  std::sort(subcell_data.begin(), subcell_data.end(), compare);
5933 
5934  // Finally, determine which objects are duplicates. Duplicates are
5935  // assumed to be interior objects, so delete all but one and change the
5936  // boundary id:
5937  auto left = subcell_data.begin();
5938  while (left != subcell_data.end())
5939  {
5940  const auto right =
5941  std::upper_bound(left, subcell_data.end(), *left, compare);
5942  // if the range has more than one item, then there are duplicates -
5943  // set all boundary ids in the range to the internal boundary id
5944  if (left + 1 != right)
5945  for (auto it = left; it != right; ++it)
5946  {
5948  Assert(it->manifold_id == left->manifold_id,
5949  ExcMessage(
5950  "In the process of grid generation a single "
5951  "line or quadrilateral has been assigned two "
5952  "different manifold ids. This can happen when "
5953  "a Triangulation is copied, e.g., via "
5954  "GridGenerator::replicate_triangulation() and "
5955  "not all external boundary faces have the same "
5956  "manifold id. Double check that all faces "
5957  "which you expect to be merged together have "
5958  "the same manifold id."));
5959  }
5960  left = right;
5961  }
5962 
5963  subcell_data.erase(std::unique(subcell_data.begin(), subcell_data.end()),
5964  subcell_data.end());
5965  }
5966  } // namespace
5967 
5968 
5969 
5970  template <int dim, int spacedim>
5971  void
5973  const std::vector<unsigned int> & extents,
5975  {
5976  AssertDimension(dim, extents.size());
5977 #ifdef DEBUG
5978  for (const auto &extent : extents)
5979  Assert(0 < extent,
5980  ExcMessage("The Triangulation must be copied at least one time in "
5981  "each coordinate dimension."));
5982 #endif
5983  const BoundingBox<spacedim> bbox(input.get_vertices());
5984  const auto & min = bbox.get_boundary_points().first;
5985  const auto & max = bbox.get_boundary_points().second;
5986 
5987  std::array<Tensor<1, spacedim>, dim> offsets;
5988  for (unsigned int d = 0; d < dim; ++d)
5989  offsets[d][d] = max[d] - min[d];
5990 
5991  Triangulation<dim, spacedim> tria_to_replicate;
5992  tria_to_replicate.copy_triangulation(input);
5993  for (unsigned int d = 0; d < dim; ++d)
5994  {
5995  std::vector<Point<spacedim>> input_vertices;
5996  std::vector<CellData<dim>> input_cell_data;
5997  SubCellData input_subcell_data;
5998  std::tie(input_vertices, input_cell_data, input_subcell_data) =
5999  GridTools::get_coarse_mesh_description(tria_to_replicate);
6000  std::vector<Point<spacedim>> output_vertices = input_vertices;
6001  std::vector<CellData<dim>> output_cell_data = input_cell_data;
6002  SubCellData output_subcell_data = input_subcell_data;
6003 
6004  for (unsigned int k = 1; k < extents[d]; ++k)
6005  {
6006  const std::size_t vertex_offset = k * input_vertices.size();
6007  // vertices
6008  for (const Point<spacedim> &point : input_vertices)
6009  output_vertices.push_back(point + double(k) * offsets[d]);
6010  // cell data
6011  for (const CellData<dim> &cell_data : input_cell_data)
6012  {
6013  output_cell_data.push_back(cell_data);
6014  for (unsigned int &vertex : output_cell_data.back().vertices)
6015  vertex += vertex_offset;
6016  }
6017  // subcell data
6018  for (const CellData<1> &boundary_line :
6019  input_subcell_data.boundary_lines)
6020  {
6021  output_subcell_data.boundary_lines.push_back(boundary_line);
6022  for (unsigned int &vertex :
6023  output_subcell_data.boundary_lines.back().vertices)
6024  vertex += vertex_offset;
6025  }
6026  for (const CellData<2> &boundary_quad :
6027  input_subcell_data.boundary_quads)
6028  {
6029  output_subcell_data.boundary_quads.push_back(boundary_quad);
6030  for (unsigned int &vertex :
6031  output_subcell_data.boundary_quads.back().vertices)
6032  vertex += vertex_offset;
6033  }
6034  }
6035  // check all vertices: since the grid is coarse, most will be on the
6036  // boundary anyway
6037  std::vector<unsigned int> boundary_vertices;
6039  output_vertices,
6040  output_cell_data,
6041  output_subcell_data,
6042  boundary_vertices,
6043  1e-6 * input.begin_active()->diameter());
6044  // delete_duplicated_vertices also deletes any unused vertices
6045  // deal with any reordering issues created by delete_duplicated_vertices
6046  GridReordering<dim>::reorder_cells(output_cell_data, true);
6047  // clean up the boundary ids of the boundary objects: note that we
6048  // have to do this after delete_duplicated_vertices so that boundary
6049  // objects are actually duplicated at this point
6050  if (dim == 2)
6051  delete_duplicated_objects(output_subcell_data.boundary_lines);
6052  else if (dim == 3)
6053  {
6054  delete_duplicated_objects(output_subcell_data.boundary_quads);
6055  for (CellData<1> &boundary_line :
6056  output_subcell_data.boundary_lines)
6057  // set boundary lines to the default value - let
6058  // create_triangulation figure out the rest.
6060  }
6061 
6062  tria_to_replicate.clear();
6063  tria_to_replicate.create_triangulation(output_vertices,
6064  output_cell_data,
6065  output_subcell_data);
6066  }
6067 
6068  result.copy_triangulation(tria_to_replicate);
6069  }
6070 
6071 
6072 
6073  template <int dim, int spacedim>
6074  void
6076  const Triangulation<dim, spacedim> &triangulation_1,
6077  const Triangulation<dim, spacedim> &triangulation_2,
6079  {
6080  Assert(GridTools::have_same_coarse_mesh(triangulation_1, triangulation_2),
6081  ExcMessage("The two input triangulations are not derived from "
6082  "the same coarse mesh as required."));
6083  Assert((dynamic_cast<
6085  &triangulation_1) == nullptr) &&
6086  (dynamic_cast<
6088  &triangulation_2) == nullptr),
6089  ExcMessage("The source triangulations for this function must both "
6090  "be available entirely locally, and not be distributed "
6091  "triangulations."));
6092 
6093  // first copy triangulation_1, and
6094  // then do as many iterations as
6095  // there are levels in
6096  // triangulation_2 to refine
6097  // additional cells. since this is
6098  // the maximum number of
6099  // refinements to get from the
6100  // coarse grid to triangulation_2,
6101  // it is clear that this is also
6102  // the maximum number of
6103  // refinements to get from any cell
6104  // on triangulation_1 to
6105  // triangulation_2
6106  result.clear();
6107  result.copy_triangulation(triangulation_1);
6108  for (unsigned int iteration = 0; iteration < triangulation_2.n_levels();
6109  ++iteration)
6110  {
6112  intergrid_map.make_mapping(result, triangulation_2);
6113 
6114  bool any_cell_flagged = false;
6116  result_cell = result.begin_active();
6117  result_cell != result.end();
6118  ++result_cell)
6119  if (intergrid_map[result_cell]->has_children())
6120  {
6121  any_cell_flagged = true;
6122  result_cell->set_refine_flag();
6123  }
6124 
6125  if (any_cell_flagged == false)
6126  break;
6127  else
6129  }
6130  }
6131 
6132 
6133 
6134  template <int dim, int spacedim>
6135  void
6137  const Triangulation<dim, spacedim> &input_triangulation,
6139  & cells_to_remove,
6141  {
6142  // simply copy the vertices; we will later strip those
6143  // that turn out to be unused
6144  std::vector<Point<spacedim>> vertices = input_triangulation.get_vertices();
6145 
6146  // the loop through the cells and copy stuff, excluding
6147  // the ones we are to remove
6148  std::vector<CellData<dim>> cells;
6150  input_triangulation.begin_active();
6151  cell != input_triangulation.end();
6152  ++cell)
6153  if (cells_to_remove.find(cell) == cells_to_remove.end())
6154  {
6155  Assert(static_cast<unsigned int>(cell->level()) ==
6156  input_triangulation.n_levels() - 1,
6157  ExcMessage(
6158  "Your input triangulation appears to have "
6159  "adaptively refined cells. This is not allowed. You can "
6160  "only call this function on a triangulation in which "
6161  "all cells are on the same refinement level."));
6162 
6163  CellData<dim> this_cell;
6164  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
6165  ++v)
6166  this_cell.vertices[v] = cell->vertex_index(v);
6167  this_cell.material_id = cell->material_id();
6168  cells.push_back(this_cell);
6169  }
6170 
6171  // throw out duplicated vertices from the two meshes, reorder vertices as
6172  // necessary and create the triangulation
6173  SubCellData subcell_data;
6174  std::vector<unsigned int> considered_vertices;
6176  cells,
6177  subcell_data,
6178  considered_vertices);
6179 
6180  // then clear the old triangulation and create the new one
6181  result.clear();
6182  result.create_triangulation(vertices, cells, subcell_data);
6183  }
6184 
6185 
6186 
6187  void
6189  const Triangulation<2, 2> & input,
6190  const unsigned int n_slices,
6191  const double height,
6192  Triangulation<3, 3> & result,
6193  const bool copy_manifold_ids,
6194  const std::vector<types::manifold_id> &manifold_priorities)
6195  {
6196  Assert(input.n_levels() == 1,
6197  ExcMessage(
6198  "The input triangulation must be a coarse mesh, i.e., it must "
6199  "not have been refined."));
6200  Assert(result.n_cells() == 0,
6201  ExcMessage("The output triangulation object needs to be empty."));
6202  Assert(height > 0,
6203  ExcMessage("The given height for extrusion must be positive."));
6204  Assert(n_slices >= 2,
6205  ExcMessage(
6206  "The number of slices for extrusion must be at least 2."));
6207 
6208  const double delta_h = height / (n_slices - 1);
6209  std::vector<double> slices_z_values;
6210  for (unsigned int i = 0; i < n_slices; ++i)
6211  slices_z_values.push_back(i * delta_h);
6213  input, slices_z_values, result, copy_manifold_ids, manifold_priorities);
6214  }
6215 
6216 
6217 
6218  void
6220  const Triangulation<2, 2> & input,
6221  const std::vector<double> & slice_coordinates,
6222  Triangulation<3, 3> & result,
6223  const bool copy_manifold_ids,
6224  const std::vector<types::manifold_id> &manifold_priorities)
6225  {
6226  Assert(input.n_levels() == 1,
6227  ExcMessage(
6228  "The input triangulation must be a coarse mesh, i.e., it must "
6229  "not have been refined."));
6230  Assert(result.n_cells() == 0,
6231  ExcMessage("The output triangulation object needs to be empty."));
6232  Assert(slice_coordinates.size() >= 2,
6233  ExcMessage(
6234  "The number of slices for extrusion must be at least 2."));
6235  Assert(std::is_sorted(slice_coordinates.begin(), slice_coordinates.end()),
6236  ExcMessage("Slice z-coordinates should be in ascending order"));
6237 
6238  const auto priorities = [&]() -> std::vector<types::manifold_id> {
6239  // if a non-empty (i.e., not the default) vector is given for
6240  // manifold_priorities then use it (but check its validity in debug
6241  // mode)
6242  if (0 < manifold_priorities.size())
6243  {
6244 #ifdef DEBUG
6245  // check that the provided manifold_priorities is valid
6246  std::vector<types::manifold_id> sorted_manifold_priorities =
6247  manifold_priorities;
6248  std::sort(sorted_manifold_priorities.begin(),
6249  sorted_manifold_priorities.end());
6250  Assert(std::unique(sorted_manifold_priorities.begin(),
6251  sorted_manifold_priorities.end()) ==
6252  sorted_manifold_priorities.end(),
6253  ExcMessage(
6254  "The given vector of manifold ids may not contain any "
6255  "duplicated entries."));
6256  std::vector<types::manifold_id> sorted_manifold_ids =
6257  input.get_manifold_ids();
6258  std::sort(sorted_manifold_ids.begin(), sorted_manifold_ids.end());
6259  if (sorted_manifold_priorities != sorted_manifold_ids)
6260  {
6261  std::ostringstream message;
6262  message << "The given triangulation has manifold ids {";
6263  for (const types::manifold_id manifold_id : sorted_manifold_ids)
6264  if (manifold_id != sorted_manifold_ids.back())
6265  message << manifold_id << ", ";
6266  message << sorted_manifold_ids.back() << "}, but \n"
6267  << " the given vector of manifold ids is {";
6268  for (const types::manifold_id manifold_id : manifold_priorities)
6269  if (manifold_id != manifold_priorities.back())
6270  message << manifold_id << ", ";
6271  message
6272  << manifold_priorities.back() << "}.\n"
6273  << " These vectors should contain the same elements.\n";
6274  const std::string m = message.str();
6275  Assert(false, ExcMessage(m));
6276  }
6277 #endif
6278  return manifold_priorities;
6279  }
6280  // otherwise use the default ranking: ascending order, but TFI manifolds
6281  // are at the end.
6282  std::vector<types::manifold_id> default_priorities =
6283  input.get_manifold_ids();
6284  const auto first_tfi_it = std::partition(
6285  default_priorities.begin(),
6286  default_priorities.end(),
6287  [&input](const types::manifold_id &id) {
6288  return dynamic_cast<const TransfiniteInterpolationManifold<2, 2> *>(
6289  &input.get_manifold(id)) == nullptr;
6290  });
6291  std::sort(default_priorities.begin(), first_tfi_it);
6292  std::sort(first_tfi_it, default_priorities.end());
6293 
6294  return default_priorities;
6295  }();
6296 
6297  const std::size_t n_slices = slice_coordinates.size();
6298  std::vector<Point<3>> points(n_slices * input.n_vertices());
6299  std::vector<CellData<3>> cells;
6300  cells.reserve((n_slices - 1) * input.n_active_cells());
6301 
6302  // copy the array of points as many times as there will be slices,
6303  // one slice at a time. The z-axis value are defined in slices_coordinates
6304  for (std::size_t slice_n = 0; slice_n < n_slices; ++slice_n)
6305  {
6306  for (std::size_t vertex_n = 0; vertex_n < input.n_vertices();
6307  ++vertex_n)
6308  {
6309  const Point<2> vertex = input.get_vertices()[vertex_n];
6310  points[slice_n * input.n_vertices() + vertex_n] =
6311  Point<3>(vertex[0], vertex[1], slice_coordinates[slice_n]);
6312  }
6313  }
6314 
6315  // then create the cells of each of the slices, one stack at a
6316  // time
6317  for (const auto &cell : input.active_cell_iterators())
6318  {
6319  for (std::size_t slice_n = 0; slice_n < n_slices - 1; ++slice_n)
6320  {
6321  CellData<3> this_cell;
6322  for (unsigned int vertex_n = 0;
6323  vertex_n < GeometryInfo<2>::vertices_per_cell;
6324  ++vertex_n)
6325  {
6326  this_cell.vertices[vertex_n] =
6327  cell->vertex_index(vertex_n) + slice_n * input.n_vertices();
6328  this_cell
6330  cell->vertex_index(vertex_n) +
6331  (slice_n + 1) * input.n_vertices();
6332  }
6333 
6334  this_cell.material_id = cell->material_id();
6335  if (copy_manifold_ids)
6336  this_cell.manifold_id = cell->manifold_id();
6337  cells.push_back(this_cell);
6338  }
6339  }
6340 
6341  // Next, create face data for all faces that are orthogonal to the x-y
6342  // plane
6343  SubCellData subcell_data;
6344  std::vector<CellData<2>> &quads = subcell_data.boundary_quads;
6345  types::boundary_id max_boundary_id = 0;
6346  quads.reserve(input.n_active_lines() * (n_slices - 1) +
6347  input.n_active_cells() * 2);
6348  for (const auto &face : input.active_face_iterators())
6349  {
6350  CellData<2> quad;
6351  quad.boundary_id = face->boundary_id();
6352  if (face->at_boundary())
6353  max_boundary_id = std::max(max_boundary_id, quad.boundary_id);
6354  if (copy_manifold_ids)
6355  quad.manifold_id = face->manifold_id();
6356  for (std::size_t slice_n = 0; slice_n < n_slices - 1; ++slice_n)
6357  {
6358  quad.vertices[0] =
6359  face->vertex_index(0) + slice_n * input.n_vertices();
6360  quad.vertices[1] =
6361  face->vertex_index(1) + slice_n * input.n_vertices();
6362  quad.vertices[2] =
6363  face->vertex_index(0) + (slice_n + 1) * input.n_vertices();
6364  quad.vertices[3] =
6365  face->vertex_index(1) + (slice_n + 1) * input.n_vertices();
6366  quads.push_back(quad);
6367  }
6368  }
6369 
6370  // if necessary, create face data for faces parallel to the x-y
6371  // plane. This is only necessary if we need to set manifolds.
6372  if (copy_manifold_ids)
6373  for (const auto &cell : input.active_cell_iterators())
6374  {
6375  CellData<2> quad;
6377  quad.manifold_id = cell->manifold_id(); // check is outside loop
6378  for (std::size_t slice_n = 1; slice_n < n_slices - 1; ++slice_n)
6379  {
6380  quad.vertices[0] =
6381  cell->vertex_index(0) + slice_n * input.n_vertices();
6382  quad.vertices[1] =
6383  cell->vertex_index(1) + slice_n * input.n_vertices();
6384  quad.vertices[2] =
6385  cell->vertex_index(2) + slice_n * input.n_vertices();
6386  quad.vertices[3] =
6387  cell->vertex_index(3) + slice_n * input.n_vertices();
6388  quads.push_back(quad);
6389  }
6390  }
6391 
6392  // then mark the bottom and top boundaries of the extruded mesh
6393  // with max_boundary_id+1 and max_boundary_id+2. check that this
6394  // remains valid
6395  Assert((max_boundary_id != numbers::invalid_boundary_id) &&
6396  (max_boundary_id + 1 != numbers::invalid_boundary_id) &&
6397  (max_boundary_id + 2 != numbers::invalid_boundary_id),
6398  ExcMessage(
6399  "The input triangulation to this function is using boundary "
6400  "indicators in a range that do not allow using "
6401  "max_boundary_id+1 and max_boundary_id+2 as boundary "
6402  "indicators for the bottom and top faces of the "
6403  "extruded triangulation."));
6404  const types::boundary_id bottom_boundary_id = max_boundary_id + 1;
6405  const types::boundary_id top_boundary_id = max_boundary_id + 2;
6406  for (const auto &cell : input.active_cell_iterators())
6407  {
6408  CellData<2> quad;
6409  quad.boundary_id = bottom_boundary_id;
6410  quad.vertices[0] = cell->vertex_index(0);
6411  quad.vertices[1] = cell->vertex_index(1);
6412  quad.vertices[2] = cell->vertex_index(2);
6413  quad.vertices[3] = cell->vertex_index(3);
6414  if (copy_manifold_ids)
6415  quad.manifold_id = cell->manifold_id();
6416  quads.push_back(quad);
6417 
6418  quad.boundary_id = top_boundary_id;
6419  for (unsigned int &vertex : quad.vertices)
6420  vertex += (n_slices - 1) * input.n_vertices();
6421  if (copy_manifold_ids)
6422  quad.manifold_id = cell->manifold_id();
6423  quads.push_back(quad);
6424  }
6425 
6426  // use all of this to finally create the extruded 3d
6427  // triangulation. it is not necessary to call
6428  // GridReordering<3,3>::reorder_cells because the cells we have
6429  // constructed above are automatically correctly oriented. this is
6430  // because the 2d base mesh is always correctly oriented, and
6431  // extruding it automatically yields a correctly oriented 3d mesh,
6432  // as discussed in the edge orientation paper mentioned in the
6433  // introduction to the GridReordering class.
6434  result.create_triangulation(points, cells, subcell_data);
6435 
6436  for (auto manifold_id_it = priorities.rbegin();
6437  manifold_id_it != priorities.rend();
6438  ++manifold_id_it)
6439  for (const auto &face : result.active_face_iterators())
6440  if (face->manifold_id() == *manifold_id_it)
6441  for (unsigned int line_n = 0;
6442  line_n < GeometryInfo<3>::lines_per_face;
6443  ++line_n)
6444  face->line(line_n)->set_manifold_id(*manifold_id_it);
6445  }
6446 
6447 
6448  template <>
6450  const double,
6451  const double,
6452  const double,
6453  const unsigned int,
6454  const bool)
6455  {
6456  Assert(false, ExcNotImplemented());
6457  }
6458 
6459 
6460 
6461  template <>
6463  const double inner_radius,
6464  const double outer_radius,
6465  const double, // width,
6466  const unsigned int, // width_repetition,
6467  const bool colorize)
6468  {
6469  const int dim = 2;
6470 
6471  Assert(inner_radius < outer_radius,
6472  ExcMessage("outer_radius has to be bigger than inner_radius."));
6473 
6474  Point<dim> center;
6475  // We create an hyper_shell in two dimensions, and then we modify it.
6476  hyper_shell(triangulation, center, inner_radius, outer_radius, 8);
6479  triangulation.begin_active(),
6480  endc = triangulation.end();
6481  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
6482  for (; cell != endc; ++cell)
6483  {
6484  for (auto f : GeometryInfo<dim>::face_indices())
6485  if (cell->face(f)->at_boundary())
6486  {
6487  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
6488  ++v)
6489  {
6490  unsigned int vv = cell->face(f)->vertex_index(v);
6491  if (treated_vertices[vv] == false)
6492  {
6493  treated_vertices[vv] = true;
6494  switch (vv)
6495  {
6496  case 1:
6497  cell->face(f)->vertex(v) =
6498  center + Point<dim>(outer_radius, outer_radius);
6499  break;
6500  case 3:
6501  cell->face(f)->vertex(v) =
6502  center + Point<dim>(-outer_radius, outer_radius);
6503  break;
6504  case 5:
6505  cell->face(f)->vertex(v) =
6506  center + Point<dim>(-outer_radius, -outer_radius);
6507  break;
6508  case 7:
6509  cell->face(f)->vertex(v) =
6510  center + Point<dim>(outer_radius, -outer_radius);
6511  break;
6512  default:
6513  break;
6514  }
6515  }
6516  }
6517  }
6518  }
6519  double eps = 1e-3 * outer_radius;
6520  cell = triangulation.begin_active();
6521  for (; cell != endc; ++cell)
6522  {
6523  for (auto f : GeometryInfo<dim>::face_indices())
6524  if (cell->face(f)->at_boundary())
6525  {
6526  double dx = cell->face(f)->center()(0) - center(0);
6527  double dy = cell->face(f)->center()(1) - center(1);
6528  if (colorize)
6529  {
6530  if (std::abs(dx + outer_radius) < eps)
6531  cell->face(f)->set_boundary_id(0);
6532  else if (std::abs(dx - outer_radius) < eps)
6533  cell->face(f)->set_boundary_id(1);
6534  else if (std::abs(dy + outer_radius) < eps)
6535  cell->face(f)->set_boundary_id(2);
6536  else if (std::abs(dy - outer_radius) < eps)
6537  cell->face(f)->set_boundary_id(3);
6538  else
6539  {
6540  cell->face(f)->set_boundary_id(4);
6541  cell->face(f)->set_manifold_id(0);
6542  }
6543  }
6544  else
6545  {
6546  double d = (cell->face(f)->center() - center).norm();
6547  if (d - inner_radius < 0)
6548  {
6549  cell->face(f)->set_boundary_id(1);
6550  cell->face(f)->set_manifold_id(0);
6551  }
6552  else
6553  cell->face(f)->set_boundary_id(0);
6554  }
6555  }
6556  }
6557  triangulation.set_manifold(0, PolarManifold<2>(center));
6558  }
6559 
6560 
6561 
6562  template <int dim>
6563  void
6565  const Point<dim> & center,
6566  const double inner_radius,
6567  const double outer_radius,
6568  const unsigned int n_shells,
6569  const double skewness,
6570  const unsigned int n_cells,
6571  const bool colorize)
6572  {
6573  Assert(dim == 2 || dim == 3, ExcNotImplemented());
6574  (void)colorize;
6575  (void)n_cells;
6576  Assert(inner_radius < outer_radius,
6577  ExcMessage("outer_radius has to be bigger than inner_radius."));
6578  if (n_shells == 0)
6579  return; // empty Triangulation
6580 
6581  std::vector<double> radii;
6582  radii.push_back(inner_radius);
6583  for (unsigned int shell_n = 1; shell_n < n_shells; ++shell_n)
6584  if (skewness == 0.0)
6585  // same as below, but works in the limiting case of zero skewness
6586  radii.push_back(inner_radius +
6587  (outer_radius - inner_radius) *
6588  (1.0 - (1.0 - double(shell_n) / n_shells)));
6589  else
6590  radii.push_back(
6591  inner_radius +
6592  (outer_radius - inner_radius) *
6593  (1.0 - std::tanh(skewness * (1.0 - double(shell_n) / n_shells)) /
6594  std::tanh(skewness)));
6595  radii.push_back(outer_radius);
6596 
6597  double grid_vertex_tolerance = 0.0;
6598  for (unsigned int shell_n = 0; shell_n < radii.size() - 1; ++shell_n)
6599  {
6600  Triangulation<dim> current_shell;
6601  GridGenerator::hyper_shell(current_shell,
6602  center,
6603  radii[shell_n],
6604  radii[shell_n + 1],
6605  n_cells == 0 ? (dim == 2 ? 8 : 12) :
6606  n_cells);
6607 
6608  // The innermost shell has the smallest cells: use that to set the
6609  // vertex merging tolerance
6610  if (grid_vertex_tolerance == 0.0)
6611  grid_vertex_tolerance =
6612  0.5 * internal::minimal_vertex_distance(current_shell);
6613 
6614  Triangulation<dim> temp(std::move(triangulation));
6615  triangulation.clear();
6617  temp,
6618  triangulation,
6619  grid_vertex_tolerance);
6620  }
6621 
6622  const types::manifold_id manifold_id = 0;
6623  triangulation.set_all_manifold_ids(manifold_id);
6624  if (dim == 2)
6625  triangulation.set_manifold(manifold_id, PolarManifold<dim>(center));
6626  else if (dim == 3)
6627  triangulation.set_manifold(manifold_id, SphericalManifold<dim>(center));
6628 
6629  // We use boundary vertex positions to see if things are on the inner or
6630  // outer boundary.
6631  constexpr double radial_vertex_tolerance =
6632  100.0 * std::numeric_limits<double>::epsilon();
6633  auto assert_vertex_distance_within_tolerance =
6634  [center, radial_vertex_tolerance](
6635  const TriaIterator<TriaAccessor<dim - 1, dim, dim>> face,
6636  const double radius) {
6637  (void)center;
6638  (void)radial_vertex_tolerance;
6639  (void)face;
6640  (void)radius;
6641  for (unsigned int vertex_n = 0;
6642  vertex_n < GeometryInfo<dim>::vertices_per_face;
6643  ++vertex_n)
6644  {
6645  Assert(std::abs((face->vertex(vertex_n) - center).norm() - radius) <
6646  (center.norm() + radius) * radial_vertex_tolerance,
6647  ExcInternalError());
6648  }
6649  };
6650  if (colorize)
6651  for (const auto &cell : triangulation.active_cell_iterators())
6652  for (unsigned int face_n = 0;
6653  face_n < GeometryInfo<dim>::faces_per_cell;
6654  ++face_n)
6655  {
6656  auto face = cell->face(face_n);
6657  if (face->at_boundary())
6658  {
6659  if (((face->vertex(0) - center).norm() - inner_radius) <
6660  (center.norm() + inner_radius) * radial_vertex_tolerance)
6661  {
6662  // we must be at an inner face, but check
6663  assert_vertex_distance_within_tolerance(face, inner_radius);
6664  face->set_all_boundary_ids(0);
6665  }
6666  else
6667  {
6668  // we must be at an outer face, but check
6669  assert_vertex_distance_within_tolerance(face, outer_radius);
6670  face->set_all_boundary_ids(1);
6671  }
6672  }
6673  }
6674  }
6675 
6676 
6677 
6678  template <>
6680  const double inner_radius,
6681  const double outer_radius,
6682  const double L,
6683  const unsigned int Nz,
6684  const bool colorize)
6685  {
6686  const int dim = 3;
6687 
6688  Assert(inner_radius < outer_radius,
6689  ExcMessage("outer_radius has to be bigger than inner_radius."));
6690  Assert(L > 0, ExcMessage("Must give positive extension L"));
6691  Assert(Nz >= 1, ExcLowerRange(1, Nz));
6692 
6693  cylinder_shell(triangulation, L, inner_radius, outer_radius, 8, Nz);
6694  triangulation.set_all_manifold_ids