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