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