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