Reference documentation for deal.II version 9.1.0-pre
grid_generator.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/distributed/shared_tria.h>
17 #include <deal.II/distributed/tria.h>
18 
19 #include <deal.II/grid/grid_generator.h>
20 #include <deal.II/grid/grid_reordering.h>
21 #include <deal.II/grid/grid_tools.h>
22 #include <deal.II/grid/intergrid_map.h>
23 #include <deal.II/grid/manifold_lib.h>
24 #include <deal.II/grid/tria.h>
25 #include <deal.II/grid/tria_accessor.h>
26 #include <deal.II/grid/tria_iterator.h>
27 
28 #include <cmath>
29 #include <iostream>
30 #include <limits>
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 
35 namespace GridGenerator
36 {
37  namespace
38  {
39  // Corner points of the cube [-1,1]^3
40  const Point<3> hexahedron[8] = {Point<3>(-1, -1, -1),
41  Point<3>(+1, -1, -1),
42  Point<3>(-1, +1, -1),
43  Point<3>(+1, +1, -1),
44  Point<3>(-1, -1, +1),
45  Point<3>(+1, -1, +1),
46  Point<3>(-1, +1, +1),
47  Point<3>(+1, +1, +1)};
48 
49  // Octahedron inscribed in the cube
50  // [-1,1]^3
51  const Point<3> octahedron[6] = {Point<3>(-1, 0, 0),
52  Point<3>(1, 0, 0),
53  Point<3>(0, -1, 0),
54  Point<3>(0, 1, 0),
55  Point<3>(0, 0, -1),
56  Point<3>(0, 0, 1)};
57 
58 
63  template <int dim, int spacedim>
64  void
65  colorize_hyper_rectangle(Triangulation<dim, spacedim> &tria)
66  {
67  // there is nothing to do in 1d
68  if (dim > 1)
69  {
70  // there is only one cell, so
71  // simple task
72  const typename Triangulation<dim, spacedim>::cell_iterator cell =
73  tria.begin();
74  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
75  cell->face(f)->set_boundary_id(f);
76  }
77  }
78 
79 
80 
81  template <int spacedim>
82  void colorize_subdivided_hyper_rectangle(Triangulation<1, spacedim> &tria,
83  const Point<spacedim> &,
84  const Point<spacedim> &,
85  const double)
86  {
87  for (typename Triangulation<1, spacedim>::cell_iterator cell =
88  tria.begin();
89  cell != tria.end();
90  ++cell)
91  if (cell->center()(0) > 0)
92  cell->set_material_id(1);
93  // boundary indicators are set to
94  // 0 (left) and 1 (right) by default.
95  }
96 
97 
98 
99  template <int dim, int spacedim>
100  void
101  colorize_subdivided_hyper_rectangle(Triangulation<dim, spacedim> &tria,
102  const Point<spacedim> & p1,
103  const Point<spacedim> & p2,
104  const double epsilon)
105  {
106  // run through all faces and check
107  // if one of their center coordinates matches
108  // one of the corner points. Comparisons
109  // are made using an epsilon which
110  // should be smaller than the smallest cell
111  // diameter.
112 
113  typename Triangulation<dim, spacedim>::face_iterator face =
114  tria.begin_face(),
115  endface =
116  tria.end_face();
117  for (; face != endface; ++face)
118  if (face->at_boundary())
119  if (face->boundary_id() == 0)
120  {
121  const Point<spacedim> center(face->center());
122 
123  if (std::abs(center(0) - p1[0]) < epsilon)
124  face->set_boundary_id(0);
125  else if (std::abs(center(0) - p2[0]) < epsilon)
126  face->set_boundary_id(1);
127  else if (dim > 1 && std::abs(center(1) - p1[1]) < epsilon)
128  face->set_boundary_id(2);
129  else if (dim > 1 && std::abs(center(1) - p2[1]) < epsilon)
130  face->set_boundary_id(3);
131  else if (dim > 2 && std::abs(center(2) - p1[2]) < epsilon)
132  face->set_boundary_id(4);
133  else if (dim > 2 && std::abs(center(2) - p2[2]) < epsilon)
134  face->set_boundary_id(5);
135  else
136  // triangulation says it
137  // is on the boundary,
138  // but we could not find
139  // on which boundary.
140  Assert(false, ExcInternalError());
141  }
142 
143  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
144  tria.begin();
145  cell != tria.end();
146  ++cell)
147  {
148  char id = 0;
149  for (unsigned int d = 0; d < dim; ++d)
150  if (cell->center()(d) > 0)
151  id += (1 << d);
152  cell->set_material_id(id);
153  }
154  }
155 
156 
161  void colorize_hyper_shell(Triangulation<2> &tria,
162  const Point<2> &,
163  const double,
164  const double)
165  {
166  // In spite of receiving geometrical
167  // data, we do this only based on
168  // topology.
169 
170  // For the mesh based on cube,
171  // this is highly irregular
172  for (Triangulation<2>::cell_iterator cell = tria.begin();
173  cell != tria.end();
174  ++cell)
175  {
176  Assert(cell->face(2)->at_boundary(), ExcInternalError());
177  cell->face(2)->set_all_boundary_ids(1);
178  }
179  }
180 
181 
186  void colorize_hyper_shell(Triangulation<3> &tria,
187  const Point<3> &,
188  const double,
189  const double)
190  {
191  // the following uses a good amount
192  // of knowledge about the
193  // orientation of cells. this is
194  // probably not good style...
195  if (tria.n_cells() == 6)
196  {
197  Triangulation<3>::cell_iterator cell = tria.begin();
198 
199  Assert(cell->face(4)->at_boundary(), ExcInternalError());
200  cell->face(4)->set_all_boundary_ids(1);
201 
202  ++cell;
203  Assert(cell->face(2)->at_boundary(), ExcInternalError());
204  cell->face(2)->set_all_boundary_ids(1);
205 
206  ++cell;
207  Assert(cell->face(2)->at_boundary(), ExcInternalError());
208  cell->face(2)->set_all_boundary_ids(1);
209 
210  ++cell;
211  Assert(cell->face(0)->at_boundary(), ExcInternalError());
212  cell->face(0)->set_all_boundary_ids(1);
213 
214  ++cell;
215  Assert(cell->face(2)->at_boundary(), ExcInternalError());
216  cell->face(2)->set_all_boundary_ids(1);
217 
218  ++cell;
219  Assert(cell->face(0)->at_boundary(), ExcInternalError());
220  cell->face(0)->set_all_boundary_ids(1);
221  }
222  else if (tria.n_cells() == 12)
223  {
224  // again use some internal
225  // knowledge
226  for (Triangulation<3>::cell_iterator cell = tria.begin();
227  cell != tria.end();
228  ++cell)
229  {
230  Assert(cell->face(5)->at_boundary(), ExcInternalError());
231  cell->face(5)->set_all_boundary_ids(1);
232  }
233  }
234  else if (tria.n_cells() == 96)
235  {
236  // the 96-cell hypershell is
237  // based on a once refined
238  // 12-cell mesh. consequently,
239  // since the outer faces all
240  // are face_no==5 above, so
241  // they are here (unless they
242  // are in the interior). Use
243  // this to assign boundary
244  // indicators, but also make
245  // sure that we encounter
246  // exactly 48 such faces
247  unsigned int count = 0;
248  for (Triangulation<3>::cell_iterator cell = tria.begin();
249  cell != tria.end();
250  ++cell)
251  if (cell->face(5)->at_boundary())
252  {
253  cell->face(5)->set_all_boundary_ids(1);
254  ++count;
255  }
256  Assert(count == 48, ExcInternalError());
257  }
258  else
259  Assert(false, ExcNotImplemented());
260  }
261 
262 
263 
269  void colorize_quarter_hyper_shell(Triangulation<3> &tria,
270  const Point<3> & center,
271  const double inner_radius,
272  const double outer_radius)
273  {
274  if (tria.n_cells() != 3)
275  AssertThrow(false, ExcNotImplemented());
276 
277  double middle = (outer_radius - inner_radius) / 2e0 + inner_radius;
278  double eps = 1e-3 * middle;
279  Triangulation<3>::cell_iterator cell = tria.begin();
280 
281  for (; cell != tria.end(); ++cell)
282  for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; ++f)
283  {
284  if (!cell->face(f)->at_boundary())
285  continue;
286 
287  double radius = cell->face(f)->center().norm() - center.norm();
288  if (std::fabs(cell->face(f)->center()(0)) <
289  eps) // x = 0 set boundary 2
290  {
291  cell->face(f)->set_boundary_id(2);
292  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
293  ++j)
294  if (cell->face(f)->line(j)->at_boundary())
295  if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
296  cell->face(f)->line(j)->vertex(1).norm()) >
297  eps)
298  cell->face(f)->line(j)->set_boundary_id(2);
299  }
300  else if (std::fabs(cell->face(f)->center()(1)) <
301  eps) // y = 0 set boundary 3
302  {
303  cell->face(f)->set_boundary_id(3);
304  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
305  ++j)
306  if (cell->face(f)->line(j)->at_boundary())
307  if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
308  cell->face(f)->line(j)->vertex(1).norm()) >
309  eps)
310  cell->face(f)->line(j)->set_boundary_id(3);
311  }
312  else if (std::fabs(cell->face(f)->center()(2)) <
313  eps) // z = 0 set boundary 4
314  {
315  cell->face(f)->set_boundary_id(4);
316  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
317  ++j)
318  if (cell->face(f)->line(j)->at_boundary())
319  if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
320  cell->face(f)->line(j)->vertex(1).norm()) >
321  eps)
322  cell->face(f)->line(j)->set_boundary_id(4);
323  }
324  else if (radius < middle) // inner radius set boundary 0
325  {
326  cell->face(f)->set_boundary_id(0);
327  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
328  ++j)
329  if (cell->face(f)->line(j)->at_boundary())
330  if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
331  cell->face(f)->line(j)->vertex(1).norm()) <
332  eps)
333  cell->face(f)->line(j)->set_boundary_id(0);
334  }
335  else if (radius > middle) // outer radius set boundary 1
336  {
337  cell->face(f)->set_boundary_id(1);
338  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
339  ++j)
340  if (cell->face(f)->line(j)->at_boundary())
341  if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
342  cell->face(f)->line(j)->vertex(1).norm()) <
343  eps)
344  cell->face(f)->line(j)->set_boundary_id(1);
345  }
346  else
347  Assert(false, ExcInternalError());
348  }
349  }
350 
351  } // namespace
352 
353 
354  template <int dim, int spacedim>
355  void
356  hyper_rectangle(Triangulation<dim, spacedim> &tria,
357  const Point<dim> & p_1,
358  const Point<dim> & p_2,
359  const bool colorize)
360  {
361  // First, extend dimensions from dim to spacedim and
362  // normalize such that p1 is lower in all coordinate
363  // directions. Additional entries will be 0.
364  Point<spacedim> p1, p2;
365  for (unsigned int i = 0; i < dim; ++i)
366  {
367  p1(i) = std::min(p_1(i), p_2(i));
368  p2(i) = std::max(p_1(i), p_2(i));
369  }
370 
371  std::vector<Point<spacedim>> vertices(GeometryInfo<dim>::vertices_per_cell);
372  switch (dim)
373  {
374  case 1:
375  vertices[0] = p1;
376  vertices[1] = p2;
377  break;
378  case 2:
379  vertices[0] = vertices[1] = p1;
380  vertices[2] = vertices[3] = p2;
381 
382  vertices[1](0) = p2(0);
383  vertices[2](0) = p1(0);
384  break;
385  case 3:
386  vertices[0] = vertices[1] = vertices[2] = vertices[3] = p1;
387  vertices[4] = vertices[5] = vertices[6] = vertices[7] = p2;
388 
389  vertices[1](0) = p2(0);
390  vertices[2](1) = p2(1);
391  vertices[3](0) = p2(0);
392  vertices[3](1) = p2(1);
393 
394  vertices[4](0) = p1(0);
395  vertices[4](1) = p1(1);
396  vertices[5](1) = p1(1);
397  vertices[6](0) = p1(0);
398 
399  break;
400  default:
401  Assert(false, ExcNotImplemented());
402  }
403 
404  // Prepare cell data
405  std::vector<CellData<dim>> cells(1);
406  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
407  cells[0].vertices[i] = i;
408  cells[0].material_id = 0;
409 
410  tria.create_triangulation(vertices, cells, SubCellData());
411 
412  // Assign boundary indicators
413  if (colorize)
414  colorize_hyper_rectangle(tria);
415  }
416 
417 
418  template <int dim, int spacedim>
419  void
420  hyper_cube(Triangulation<dim, spacedim> &tria,
421  const double left,
422  const double right,
423  const bool colorize)
424  {
425  Assert(left < right,
426  ExcMessage("Invalid left-to-right bounds of hypercube"));
427 
428  Point<dim> p1, p2;
429  for (unsigned int i = 0; i < dim; ++i)
430  {
431  p1(i) = left;
432  p2(i) = right;
433  }
434  hyper_rectangle(tria, p1, p2, colorize);
435  }
436 
437  template <int dim>
438  void
439  simplex(Triangulation<dim> &tria, const std::vector<Point<dim>> &vertices)
440  {
441  AssertDimension(vertices.size(), dim + 1);
442  Assert(dim > 1, ExcNotImplemented());
443  Assert(dim < 4, ExcNotImplemented());
444 
445 #ifdef DEBUG
446  Tensor<2, dim> vector_matrix;
447  for (unsigned int d = 0; d < dim; ++d)
448  for (unsigned int c = 1; c <= dim; ++c)
449  vector_matrix[c - 1][d] = vertices[c](d) - vertices[0](d);
450  Assert(determinant(vector_matrix) > 0.,
451  ExcMessage("Vertices of simplex must form a right handed system"));
452 #endif
453 
454  // Set up the vertices by first copying into points.
455  std::vector<Point<dim>> points = vertices;
456  Point<dim> center;
457  // Compute the edge midpoints and add up everything to compute the
458  // center point.
459  for (unsigned int i = 0; i <= dim; ++i)
460  {
461  points.push_back(0.5 * (points[i] + points[(i + 1) % (dim + 1)]));
462  center += points[i];
463  }
464  if (dim > 2)
465  {
466  // In 3D, we have some more edges to deal with
467  for (unsigned int i = 1; i < dim; ++i)
468  points.push_back(0.5 * (points[i - 1] + points[i + 1]));
469  // And we need face midpoints
470  for (unsigned int i = 0; i <= dim; ++i)
471  points.push_back(1. / 3. *
472  (points[i] + points[(i + 1) % (dim + 1)] +
473  points[(i + 2) % (dim + 1)]));
474  }
475  points.push_back((1. / (dim + 1)) * center);
476 
477  std::vector<CellData<dim>> cells(dim + 1);
478  switch (dim)
479  {
480  case 2:
481  AssertDimension(points.size(), 7);
482  cells[0].vertices[0] = 0;
483  cells[0].vertices[1] = 3;
484  cells[0].vertices[2] = 5;
485  cells[0].vertices[3] = 6;
486  cells[0].material_id = 0;
487 
488  cells[1].vertices[0] = 3;
489  cells[1].vertices[1] = 1;
490  cells[1].vertices[2] = 6;
491  cells[1].vertices[3] = 4;
492  cells[1].material_id = 0;
493 
494  cells[2].vertices[0] = 5;
495  cells[2].vertices[1] = 6;
496  cells[2].vertices[2] = 2;
497  cells[2].vertices[3] = 4;
498  cells[2].material_id = 0;
499  break;
500  case 3:
501  AssertDimension(points.size(), 15);
502  cells[0].vertices[0] = 0;
503  cells[0].vertices[1] = 4;
504  cells[0].vertices[2] = 8;
505  cells[0].vertices[3] = 10;
506  cells[0].vertices[4] = 7;
507  cells[0].vertices[5] = 13;
508  cells[0].vertices[6] = 12;
509  cells[0].vertices[7] = 14;
510  cells[0].material_id = 0;
511 
512  cells[1].vertices[0] = 4;
513  cells[1].vertices[1] = 1;
514  cells[1].vertices[2] = 10;
515  cells[1].vertices[3] = 5;
516  cells[1].vertices[4] = 13;
517  cells[1].vertices[5] = 9;
518  cells[1].vertices[6] = 14;
519  cells[1].vertices[7] = 11;
520  cells[1].material_id = 0;
521 
522  cells[2].vertices[0] = 8;
523  cells[2].vertices[1] = 10;
524  cells[2].vertices[2] = 2;
525  cells[2].vertices[3] = 5;
526  cells[2].vertices[4] = 12;
527  cells[2].vertices[5] = 14;
528  cells[2].vertices[6] = 6;
529  cells[2].vertices[7] = 11;
530  cells[2].material_id = 0;
531 
532  cells[3].vertices[0] = 7;
533  cells[3].vertices[1] = 13;
534  cells[3].vertices[2] = 12;
535  cells[3].vertices[3] = 14;
536  cells[3].vertices[4] = 3;
537  cells[3].vertices[5] = 9;
538  cells[3].vertices[6] = 6;
539  cells[3].vertices[7] = 11;
540  cells[3].material_id = 0;
541  break;
542  default:
543  Assert(false, ExcNotImplemented());
544  }
545  tria.create_triangulation(points, cells, SubCellData());
546  }
547 
548 
549  void moebius(Triangulation<3> & tria,
550  const unsigned int n_cells,
551  const unsigned int n_rotations,
552  const double R,
553  const double r)
554  {
555  const unsigned int dim = 3;
556  Assert(n_cells > 4,
557  ExcMessage(
558  "More than 4 cells are needed to create a moebius grid."));
559  Assert(r > 0 && R > 0,
560  ExcMessage("Outer and inner radius must be positive."));
561  Assert(R > r,
562  ExcMessage("Outer radius must be greater than inner radius."));
563 
564 
565  std::vector<Point<dim>> vertices(4 * n_cells);
566  double beta_step = n_rotations * numbers::PI / 2.0 / n_cells;
567  double alpha_step = 2.0 * numbers::PI / n_cells;
568 
569  for (unsigned int i = 0; i < n_cells; ++i)
570  for (unsigned int j = 0; j < 4; ++j)
571  {
572  vertices[4 * i + j][0] =
573  R * std::cos(i * alpha_step) +
574  r * std::cos(i * beta_step + j * numbers::PI / 2.0) *
575  std::cos(i * alpha_step);
576  vertices[4 * i + j][1] =
577  R * std::sin(i * alpha_step) +
578  r * std::cos(i * beta_step + j * numbers::PI / 2.0) *
579  std::sin(i * alpha_step);
580  vertices[4 * i + j][2] =
581  r * std::sin(i * beta_step + j * numbers::PI / 2.0);
582  }
583 
584  unsigned int offset = 0;
585 
586  std::vector<CellData<dim>> cells(n_cells);
587  for (unsigned int i = 0; i < n_cells; ++i)
588  {
589  for (unsigned int j = 0; j < 2; ++j)
590  {
591  cells[i].vertices[0 + 4 * j] = offset + 0 + 4 * j;
592  cells[i].vertices[1 + 4 * j] = offset + 3 + 4 * j;
593  cells[i].vertices[2 + 4 * j] = offset + 2 + 4 * j;
594  cells[i].vertices[3 + 4 * j] = offset + 1 + 4 * j;
595  }
596  offset += 4;
597  cells[i].material_id = 0;
598  }
599 
600  // now correct the last four vertices
601  cells[n_cells - 1].vertices[4] = (0 + n_rotations) % 4;
602  cells[n_cells - 1].vertices[5] = (3 + n_rotations) % 4;
603  cells[n_cells - 1].vertices[6] = (2 + n_rotations) % 4;
604  cells[n_cells - 1].vertices[7] = (1 + n_rotations) % 4;
605 
606  GridReordering<dim>::invert_all_cells_of_negative_grid(vertices, cells);
607  tria.create_triangulation_compatibility(vertices, cells, SubCellData());
608  }
609 
610 
611 
612  template <>
613  void torus<2, 3>(Triangulation<2, 3> &tria, const double R, const double r)
614  {
615  Assert(R > r,
616  ExcMessage("Outer radius R must be greater than the inner "
617  "radius r."));
618  Assert(r > 0.0, ExcMessage("The inner radius r must be positive."));
619 
620  const unsigned int dim = 2;
621  const unsigned int spacedim = 3;
622  std::vector<Point<spacedim>> vertices(16);
623 
624  vertices[0] = Point<spacedim>(R - r, 0, 0);
625  vertices[1] = Point<spacedim>(R, -r, 0);
626  vertices[2] = Point<spacedim>(R + r, 0, 0);
627  vertices[3] = Point<spacedim>(R, r, 0);
628  vertices[4] = Point<spacedim>(0, 0, R - r);
629  vertices[5] = Point<spacedim>(0, -r, R);
630  vertices[6] = Point<spacedim>(0, 0, R + r);
631  vertices[7] = Point<spacedim>(0, r, R);
632  vertices[8] = Point<spacedim>(-(R - r), 0, 0);
633  vertices[9] = Point<spacedim>(-R, -r, 0);
634  vertices[10] = Point<spacedim>(-(R + r), 0, 0);
635  vertices[11] = Point<spacedim>(-R, r, 0);
636  vertices[12] = Point<spacedim>(0, 0, -(R - r));
637  vertices[13] = Point<spacedim>(0, -r, -R);
638  vertices[14] = Point<spacedim>(0, 0, -(R + r));
639  vertices[15] = Point<spacedim>(0, r, -R);
640 
641  std::vector<CellData<dim>> cells(16);
642  // Right Hand Orientation
643  cells[0].vertices[0] = 0;
644  cells[0].vertices[1] = 4;
645  cells[0].vertices[2] = 7;
646  cells[0].vertices[3] = 3;
647  cells[0].material_id = 0;
648 
649  cells[1].vertices[0] = 1;
650  cells[1].vertices[1] = 5;
651  cells[1].vertices[2] = 4;
652  cells[1].vertices[3] = 0;
653  cells[1].material_id = 0;
654 
655  cells[2].vertices[0] = 2;
656  cells[2].vertices[1] = 6;
657  cells[2].vertices[2] = 5;
658  cells[2].vertices[3] = 1;
659  cells[2].material_id = 0;
660 
661  cells[3].vertices[0] = 3;
662  cells[3].vertices[1] = 7;
663  cells[3].vertices[2] = 6;
664  cells[3].vertices[3] = 2;
665  cells[3].material_id = 0;
666 
667  cells[4].vertices[0] = 4;
668  cells[4].vertices[1] = 8;
669  cells[4].vertices[2] = 11;
670  cells[4].vertices[3] = 7;
671  cells[4].material_id = 0;
672 
673  cells[5].vertices[0] = 5;
674  cells[5].vertices[1] = 9;
675  cells[5].vertices[2] = 8;
676  cells[5].vertices[3] = 4;
677  cells[5].material_id = 0;
678 
679  cells[6].vertices[0] = 6;
680  cells[6].vertices[1] = 10;
681  cells[6].vertices[2] = 9;
682  cells[6].vertices[3] = 5;
683  cells[6].material_id = 0;
684 
685  cells[7].vertices[0] = 7;
686  cells[7].vertices[1] = 11;
687  cells[7].vertices[2] = 10;
688  cells[7].vertices[3] = 6;
689  cells[7].material_id = 0;
690 
691  cells[8].vertices[0] = 8;
692  cells[8].vertices[1] = 12;
693  cells[8].vertices[2] = 15;
694  cells[8].vertices[3] = 11;
695  cells[8].material_id = 0;
696 
697  cells[9].vertices[0] = 9;
698  cells[9].vertices[1] = 13;
699  cells[9].vertices[2] = 12;
700  cells[9].vertices[3] = 8;
701  cells[9].material_id = 0;
702 
703  cells[10].vertices[0] = 10;
704  cells[10].vertices[1] = 14;
705  cells[10].vertices[2] = 13;
706  cells[10].vertices[3] = 9;
707  cells[10].material_id = 0;
708 
709  cells[11].vertices[0] = 11;
710  cells[11].vertices[1] = 15;
711  cells[11].vertices[2] = 14;
712  cells[11].vertices[3] = 10;
713  cells[11].material_id = 0;
714 
715  cells[12].vertices[0] = 12;
716  cells[12].vertices[1] = 0;
717  cells[12].vertices[2] = 3;
718  cells[12].vertices[3] = 15;
719  cells[12].material_id = 0;
720 
721  cells[13].vertices[0] = 13;
722  cells[13].vertices[1] = 1;
723  cells[13].vertices[2] = 0;
724  cells[13].vertices[3] = 12;
725  cells[13].material_id = 0;
726 
727  cells[14].vertices[0] = 14;
728  cells[14].vertices[1] = 2;
729  cells[14].vertices[2] = 1;
730  cells[14].vertices[3] = 13;
731  cells[14].material_id = 0;
732 
733  cells[15].vertices[0] = 15;
734  cells[15].vertices[1] = 3;
735  cells[15].vertices[2] = 2;
736  cells[15].vertices[3] = 14;
737  cells[15].material_id = 0;
738 
739  // Must call this to be able to create a
740  // correct triangulation in dealii, read
741  // GridReordering<> doc
742  GridReordering<dim, spacedim>::reorder_cells(cells);
743  tria.create_triangulation_compatibility(vertices, cells, SubCellData());
744 
745  tria.set_all_manifold_ids(0);
746  tria.set_manifold(0, TorusManifold<2>(R, r));
747  }
748 
749  template <>
750  void torus<3, 3>(Triangulation<3, 3> &tria, const double R, const double r)
751  {
752  Assert(R > r,
753  ExcMessage("Outer radius R must be greater than the inner "
754  "radius r."));
755  Assert(r > 0.0, ExcMessage("The inner radius r must be positive."));
756 
757  // abuse the moebius function to generate a torus for us
758  GridGenerator::moebius(tria, 6 /*n_cells*/, 0 /*n_rotations*/, R, r);
759 
760  // rotate by 90 degrees around the x axis to make the torus sit in the
761  // x-z plane instead of the x-y plane to be consistent with the other
762  // torus() function.
763  GridTools::rotate(numbers::PI / 2.0, 0, tria);
764 
765  // set manifolds as documented
766  tria.set_all_manifold_ids(1);
767  tria.set_all_manifold_ids_on_boundary(0);
768  tria.set_manifold(0, TorusManifold<3>(R, r));
769  }
770 
771 
772 
773  template <int dim>
774  void
775  general_cell(Triangulation<dim> & tria,
776  const std::vector<Point<dim>> &vertices,
777  const bool colorize)
778  {
779  Assert(vertices.size() == ::GeometryInfo<dim>::vertices_per_cell,
780  ExcMessage("Wrong number of vertices."));
781 
782  // First create a hyper_rectangle and then deform it.
783  hyper_cube(tria, 0, 1, colorize);
784 
785  typename Triangulation<dim>::active_cell_iterator cell =
786  tria.begin_active();
787  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
788  cell->vertex(i) = vertices[i];
789 
790  // Check that the order of the vertices makes sense, i.e., the volume of the
791  // cell is positive.
792  Assert(GridTools::volume(tria) > 0.,
793  ExcMessage(
794  "The volume of the cell is not greater than zero. "
795  "This could be due to the wrong ordering of the vertices."));
796  }
797 
798 
799 
800  template <>
801  void parallelogram(Triangulation<3> &,
802  const Point<3> (&/*corners*/)[3],
803  const bool /*colorize*/)
804  {
805  Assert(false, ExcNotImplemented());
806  }
807 
808  template <>
809  void parallelogram(Triangulation<1> &,
810  const Point<1> (&/*corners*/)[1],
811  const bool /*colorize*/)
812  {
813  Assert(false, ExcNotImplemented());
814  }
815 
816  // Implementation for 2D only
817  template <>
818  void parallelogram(Triangulation<2> &tria,
819  const Point<2> (&corners)[2],
820  const bool colorize)
821  {
822  Point<2> origin;
823  std::array<Tensor<1, 2>, 2> edges;
824  edges[0] = corners[0];
825  edges[1] = corners[1];
826  std::vector<unsigned int> subdivisions;
827  subdivided_parallelepiped<2, 2>(
828  tria, origin, edges, subdivisions, colorize);
829  }
830 
831 
832 
833  template <int dim>
834  void
835  parallelepiped(Triangulation<dim> &tria,
836  const Point<dim> (&corners)[dim],
837  const bool colorize)
838  {
839  unsigned int n_subdivisions[dim];
840  for (unsigned int i = 0; i < dim; ++i)
841  n_subdivisions[i] = 1;
842 
843  // and call the function below
844  subdivided_parallelepiped(tria, n_subdivisions, corners, colorize);
845  }
846 
847  template <int dim>
848  void
849  subdivided_parallelepiped(Triangulation<dim> &tria,
850  const unsigned int n_subdivisions,
851  const Point<dim> (&corners)[dim],
852  const bool colorize)
853  {
854  // Equalize number of subdivisions in each dim-direction, their
855  // validity will be checked later
856  unsigned int n_subdivisions_[dim];
857  for (unsigned int i = 0; i < dim; ++i)
858  n_subdivisions_[i] = n_subdivisions;
859 
860  // and call the function below
861  subdivided_parallelepiped(tria, n_subdivisions_, corners, colorize);
862  }
863 
864  template <int dim>
865  void
866  subdivided_parallelepiped(Triangulation<dim> &tria,
867 #ifndef _MSC_VER
868  const unsigned int (&n_subdivisions)[dim],
869 #else
870  const unsigned int *n_subdivisions,
871 #endif
872  const Point<dim> (&corners)[dim],
873  const bool colorize)
874  {
875  Point<dim> origin;
876  std::vector<unsigned int> subdivisions;
877  std::array<Tensor<1, dim>, dim> edges;
878  for (unsigned int i = 0; i < dim; ++i)
879  {
880  subdivisions.push_back(n_subdivisions[i]);
881  edges[i] = corners[i];
882  }
883 
884  subdivided_parallelepiped<dim, dim>(
885  tria, origin, edges, subdivisions, colorize);
886  }
887 
888  // Parallelepiped implementation in 1d, 2d, and 3d. @note The
889  // implementation in 1d is similar to hyper_rectangle(), and in 2d is
890  // similar to parallelogram().
891  //
892  // The GridReordering::reorder_grid is made use of towards the end of
893  // this function. Thus the triangulation is explicitly constructed for
894  // all dim here since it is slightly different in that respect
895  // (cf. hyper_rectangle(), parallelogram()).
896  template <int dim, int spacedim>
897  void
898  subdivided_parallelepiped(Triangulation<dim, spacedim> & tria,
899  const Point<spacedim> & origin,
900  const std::array<Tensor<1, spacedim>, dim> &edges,
901  const std::vector<unsigned int> &subdivisions,
902  const bool colorize)
903  {
904  std::vector<unsigned int> compute_subdivisions = subdivisions;
905  if (compute_subdivisions.size() == 0)
906  {
907  compute_subdivisions.resize(dim, 1);
908  }
909 
910  Assert(compute_subdivisions.size() == dim,
911  ExcMessage("One subdivision must be provided for each dimension."));
912  // check subdivisions
913  for (unsigned int i = 0; i < dim; ++i)
914  {
915  Assert(compute_subdivisions[i] > 0,
916  ExcInvalidRepetitions(subdivisions[i]));
917  Assert(
918  edges[i].norm() > 0,
919  ExcMessage(
920  "Edges in subdivided_parallelepiped() must not be degenerate."));
921  }
922 
923  /*
924  * Verify that the edge points to the right in 1D, vectors are oriented in
925  * a counter clockwise direction in 2D, or form a right handed system in
926  * 3D.
927  */
928  bool twisted_data = false;
929  switch (dim)
930  {
931  case 1:
932  {
933  twisted_data = (edges[0][0] < 0);
934  break;
935  }
936  case 2:
937  {
938  if (spacedim == 2) // this check does not make sense otherwise
939  {
940  const double plane_normal =
941  edges[0][0] * edges[1][1] - edges[0][1] * edges[1][0];
942  twisted_data = (plane_normal < 0.0);
943  }
944  break;
945  }
946  case 3:
947  {
948  // Check that the first two vectors are not linear combinations to
949  // avoid zero division later on.
950  Assert(std::abs(edges[0] * edges[1] /
951  (edges[0].norm() * edges[1].norm()) -
952  1.0) > 1.0e-15,
953  ExcMessage(
954  "Edges in subdivided_parallelepiped() must point in"
955  " different directions."));
956  const Tensor<1, spacedim> plane_normal =
957  cross_product_3d(edges[0], edges[1]);
958 
959  /*
960  * Ensure that edges 1, 2, and 3 form a right-handed set of
961  * vectors. This works by applying the definition of the dot product
962  *
963  * cos(theta) = dot(x, y)/(norm(x)*norm(y))
964  *
965  * and then, since the normal vector and third edge should both
966  * point away from the plane formed by the first two edges, the
967  * angle between them must be between 0 and pi/2; hence we just need
968  *
969  * 0 < dot(x, y).
970  */
971  twisted_data = (plane_normal * edges[2] < 0.0);
972  break;
973  }
974  default:
975  Assert(false, ExcInternalError());
976  }
977  (void)twisted_data; // make the static analyzer happy
978  Assert(
979  !twisted_data,
980  ExcInvalidInputOrientation(
981  "The triangulation you are trying to create will consist of cells"
982  " with negative measures. This is usually the result of input data"
983  " that does not define a right-handed coordinate system. The usual"
984  " fix for this is to ensure that in 1D the given point is to the"
985  " right of the origin (or the given edge tensor is positive), in 2D"
986  " that the two edges (and their cross product) obey the right-hand"
987  " rule (which may usually be done by switching the order of the"
988  " points or edge tensors), or in 3D that the edges form a"
989  " right-handed coordinate system (which may also be accomplished by"
990  " switching the order of the first two points or edge tensors)."));
991 
992  // Check corners do not overlap (unique)
993  for (unsigned int i = 0; i < dim; ++i)
994  for (unsigned int j = i + 1; j < dim; ++j)
995  Assert((edges[i] != edges[j]),
996  ExcMessage(
997  "Degenerate edges of subdivided_parallelepiped encountered."));
998 
999  // Create a list of points
1000  std::vector<Point<spacedim>> points;
1001 
1002  switch (dim)
1003  {
1004  case 1:
1005  for (unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
1006  points.push_back(origin + edges[0] / compute_subdivisions[0] * x);
1007  break;
1008 
1009  case 2:
1010  for (unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
1011  for (unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
1012  points.push_back(origin + edges[0] / compute_subdivisions[0] * x +
1013  edges[1] / compute_subdivisions[1] * y);
1014  break;
1015 
1016  case 3:
1017  for (unsigned int z = 0; z <= compute_subdivisions[2]; ++z)
1018  for (unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
1019  for (unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
1020  points.push_back(origin +
1021  edges[0] / compute_subdivisions[0] * x +
1022  edges[1] / compute_subdivisions[1] * y +
1023  edges[2] / compute_subdivisions[2] * z);
1024  break;
1025 
1026  default:
1027  Assert(false, ExcNotImplemented());
1028  }
1029 
1030  // Prepare cell data
1031  unsigned int n_cells = 1;
1032  for (unsigned int i = 0; i < dim; ++i)
1033  n_cells *= compute_subdivisions[i];
1034  std::vector<CellData<dim>> cells(n_cells);
1035 
1036  // Create fixed ordering of
1037  switch (dim)
1038  {
1039  case 1:
1040  for (unsigned int x = 0; x < compute_subdivisions[0]; ++x)
1041  {
1042  cells[x].vertices[0] = x;
1043  cells[x].vertices[1] = x + 1;
1044 
1045  // wipe material id
1046  cells[x].material_id = 0;
1047  }
1048  break;
1049 
1050  case 2:
1051  {
1052  // Shorthand
1053  const unsigned int n_dy = compute_subdivisions[1];
1054  const unsigned int n_dx = compute_subdivisions[0];
1055 
1056  for (unsigned int y = 0; y < n_dy; ++y)
1057  for (unsigned int x = 0; x < n_dx; ++x)
1058  {
1059  const unsigned int c = y * n_dx + x;
1060  cells[c].vertices[0] = y * (n_dx + 1) + x;
1061  cells[c].vertices[1] = y * (n_dx + 1) + x + 1;
1062  cells[c].vertices[2] = (y + 1) * (n_dx + 1) + x;
1063  cells[c].vertices[3] = (y + 1) * (n_dx + 1) + x + 1;
1064 
1065  // wipe material id
1066  cells[c].material_id = 0;
1067  }
1068  }
1069  break;
1070 
1071  case 3:
1072  {
1073  // Shorthand
1074  const unsigned int n_dz = compute_subdivisions[2];
1075  const unsigned int n_dy = compute_subdivisions[1];
1076  const unsigned int n_dx = compute_subdivisions[0];
1077 
1078  for (unsigned int z = 0; z < n_dz; ++z)
1079  for (unsigned int y = 0; y < n_dy; ++y)
1080  for (unsigned int x = 0; x < n_dx; ++x)
1081  {
1082  const unsigned int c = z * n_dy * n_dx + y * n_dx + x;
1083 
1084  cells[c].vertices[0] =
1085  z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
1086  cells[c].vertices[1] =
1087  z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x + 1;
1088  cells[c].vertices[2] =
1089  z * (n_dy + 1) * (n_dx + 1) + (y + 1) * (n_dx + 1) + x;
1090  cells[c].vertices[3] = z * (n_dy + 1) * (n_dx + 1) +
1091  (y + 1) * (n_dx + 1) + x + 1;
1092  cells[c].vertices[4] =
1093  (z + 1) * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
1094  cells[c].vertices[5] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
1095  y * (n_dx + 1) + x + 1;
1096  cells[c].vertices[6] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
1097  (y + 1) * (n_dx + 1) + x;
1098  cells[c].vertices[7] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
1099  (y + 1) * (n_dx + 1) + x + 1;
1100 
1101  // wipe material id
1102  cells[c].material_id = 0;
1103  }
1104  break;
1105  }
1106 
1107  default:
1108  Assert(false, ExcNotImplemented());
1109  }
1110 
1111  // Create triangulation
1112  // reorder the cells to ensure that they satisfy the convention for
1113  // edge and face directions
1114  GridReordering<dim>::reorder_cells(cells, true);
1115  tria.create_triangulation(points, cells, SubCellData());
1116 
1117  // Finally assign boundary indicators according to hyper_rectangle
1118  if (colorize)
1119  {
1120  typename Triangulation<dim>::active_cell_iterator cell =
1121  tria.begin_active(),
1122  endc = tria.end();
1123  for (; cell != endc; ++cell)
1124  {
1125  for (unsigned int face = 0;
1126  face < GeometryInfo<dim>::faces_per_cell;
1127  ++face)
1128  {
1129  if (cell->face(face)->at_boundary())
1130  cell->face(face)->set_boundary_id(face);
1131  }
1132  }
1133  }
1134  }
1135 
1136 
1137  template <int dim, int spacedim>
1138  void
1139  subdivided_hyper_cube(Triangulation<dim, spacedim> &tria,
1140  const unsigned int repetitions,
1141  const double left,
1142  const double right)
1143  {
1144  Assert(repetitions >= 1, ExcInvalidRepetitions(repetitions));
1145  Assert(left < right,
1146  ExcMessage("Invalid left-to-right bounds of hypercube"));
1147 
1148  Point<dim> p0, p1;
1149  for (unsigned int i = 0; i < dim; ++i)
1150  {
1151  p0[i] = left;
1152  p1[i] = right;
1153  }
1154 
1155  std::vector<unsigned int> reps(dim, repetitions);
1156  subdivided_hyper_rectangle(tria, reps, p0, p1);
1157  }
1158 
1159 
1160 
1161  template <int dim, int spacedim>
1162  void
1163  subdivided_hyper_rectangle(Triangulation<dim, spacedim> & tria,
1164  const std::vector<unsigned int> &repetitions,
1165  const Point<dim> & p_1,
1166  const Point<dim> & p_2,
1167  const bool colorize)
1168  {
1169  Assert(repetitions.size() == dim, ExcInvalidRepetitionsDimension(dim));
1170 
1171  // First, extend dimensions from dim to spacedim and
1172  // normalize such that p1 is lower in all coordinate
1173  // directions. Additional entries will be 0.
1174  Point<spacedim> p1, p2;
1175  for (unsigned int i = 0; i < dim; ++i)
1176  {
1177  p1(i) = std::min(p_1(i), p_2(i));
1178  p2(i) = std::max(p_1(i), p_2(i));
1179  }
1180 
1181  // calculate deltas and validate input
1182  std::vector<Point<spacedim>> delta(dim);
1183  for (unsigned int i = 0; i < dim; ++i)
1184  {
1185  Assert(repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
1186 
1187  delta[i][i] = (p2[i] - p1[i]) / repetitions[i];
1188  Assert(
1189  delta[i][i] > 0.0,
1190  ExcMessage(
1191  "The first dim entries of coordinates of p1 and p2 need to be different."));
1192  }
1193 
1194  // then generate the points
1195  std::vector<Point<spacedim>> points;
1196  switch (dim)
1197  {
1198  case 1:
1199  for (unsigned int x = 0; x <= repetitions[0]; ++x)
1200  points.push_back(p1 + (double)x * delta[0]);
1201  break;
1202 
1203  case 2:
1204  for (unsigned int y = 0; y <= repetitions[1]; ++y)
1205  for (unsigned int x = 0; x <= repetitions[0]; ++x)
1206  points.push_back(p1 + (double)x * delta[0] +
1207  (double)y * delta[1]);
1208  break;
1209 
1210  case 3:
1211  for (unsigned int z = 0; z <= repetitions[2]; ++z)
1212  for (unsigned int y = 0; y <= repetitions[1]; ++y)
1213  for (unsigned int x = 0; x <= repetitions[0]; ++x)
1214  points.push_back(p1 + (double)x * delta[0] +
1215  (double)y * delta[1] + (double)z * delta[2]);
1216  break;
1217 
1218  default:
1219  Assert(false, ExcNotImplemented());
1220  }
1221 
1222  // next create the cells
1223  std::vector<CellData<dim>> cells;
1224  switch (dim)
1225  {
1226  case 1:
1227  {
1228  cells.resize(repetitions[0]);
1229  for (unsigned int x = 0; x < repetitions[0]; ++x)
1230  {
1231  cells[x].vertices[0] = x;
1232  cells[x].vertices[1] = x + 1;
1233  cells[x].material_id = 0;
1234  }
1235  break;
1236  }
1237 
1238  case 2:
1239  {
1240  cells.resize(repetitions[1] * repetitions[0]);
1241  for (unsigned int y = 0; y < repetitions[1]; ++y)
1242  for (unsigned int x = 0; x < repetitions[0]; ++x)
1243  {
1244  const unsigned int c = x + y * repetitions[0];
1245  cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
1246  cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
1247  cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
1248  cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
1249  cells[c].material_id = 0;
1250  }
1251  break;
1252  }
1253 
1254  case 3:
1255  {
1256  const unsigned int n_x = (repetitions[0] + 1);
1257  const unsigned int n_xy =
1258  (repetitions[0] + 1) * (repetitions[1] + 1);
1259 
1260  cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
1261  for (unsigned int z = 0; z < repetitions[2]; ++z)
1262  for (unsigned int y = 0; y < repetitions[1]; ++y)
1263  for (unsigned int x = 0; x < repetitions[0]; ++x)
1264  {
1265  const unsigned int c = x + y * repetitions[0] +
1266  z * repetitions[0] * repetitions[1];
1267  cells[c].vertices[0] = z * n_xy + y * n_x + x;
1268  cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
1269  cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
1270  cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
1271  cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
1272  cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
1273  cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
1274  cells[c].vertices[7] =
1275  (z + 1) * n_xy + (y + 1) * n_x + x + 1;
1276  cells[c].material_id = 0;
1277  }
1278  break;
1279  }
1280 
1281  default:
1282  Assert(false, ExcNotImplemented());
1283  }
1284 
1285  tria.create_triangulation(points, cells, SubCellData());
1286 
1287  if (colorize)
1288  {
1289  // to colorize, run through all
1290  // faces of all cells and set
1291  // boundary indicator to the
1292  // correct value if it was 0.
1293 
1294  // use a large epsilon to
1295  // compare numbers to avoid
1296  // roundoff problems.
1297  double epsilon = 10;
1298  for (unsigned int i = 0; i < dim; ++i)
1299  epsilon = std::min(epsilon, 0.01 * delta[i][i]);
1300  Assert(epsilon > 0,
1301  ExcMessage(
1302  "The distance between corner points must be positive."))
1303 
1304  // actual code is external since
1305  // 1-D is different from 2/3D.
1306  colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
1307  }
1308  }
1309 
1310 
1311 
1312  template <int dim>
1313  void
1314  subdivided_hyper_rectangle(Triangulation<dim> & tria,
1315  const std::vector<std::vector<double>> &step_sz,
1316  const Point<dim> & p_1,
1317  const Point<dim> & p_2,
1318  const bool colorize)
1319  {
1320  Assert(step_sz.size() == dim, ExcInvalidRepetitionsDimension(dim));
1321 
1322  // First, normalize input such that
1323  // p1 is lower in all coordinate
1324  // directions and check the consistency of
1325  // step sizes, i.e. that they all
1326  // add up to the sizes specified by
1327  // p_1 and p_2
1328  Point<dim> p1(p_1);
1329  Point<dim> p2(p_2);
1330  std::vector<std::vector<double>> step_sizes(step_sz);
1331 
1332  for (unsigned int i = 0; i < dim; ++i)
1333  {
1334  if (p1(i) > p2(i))
1335  {
1336  std::swap(p1(i), p2(i));
1337  std::reverse(step_sizes[i].begin(), step_sizes[i].end());
1338  }
1339 
1340  double x = 0;
1341  for (unsigned int j = 0; j < step_sizes.at(i).size(); j++)
1342  x += step_sizes[i][j];
1343  Assert(std::fabs(x - (p2(i) - p1(i))) <= 1e-12 * std::fabs(x),
1344  ExcMessage(
1345  "The sequence of step sizes in coordinate direction " +
1346  Utilities::int_to_string(i) +
1347  " must be equal to the distance of the two given "
1348  "points in this coordinate direction."));
1349  }
1350 
1351 
1352  // then generate the necessary
1353  // points
1354  std::vector<Point<dim>> points;
1355  switch (dim)
1356  {
1357  case 1:
1358  {
1359  double x = 0;
1360  for (unsigned int i = 0;; ++i)
1361  {
1362  points.push_back(Point<dim>(p1[0] + x));
1363 
1364  // form partial sums. in
1365  // the last run through
1366  // avoid accessing
1367  // non-existent values
1368  // and exit early instead
1369  if (i == step_sizes[0].size())
1370  break;
1371 
1372  x += step_sizes[0][i];
1373  }
1374  break;
1375  }
1376 
1377  case 2:
1378  {
1379  double y = 0;
1380  for (unsigned int j = 0;; ++j)
1381  {
1382  double x = 0;
1383  for (unsigned int i = 0;; ++i)
1384  {
1385  points.push_back(Point<dim>(p1[0] + x, p1[1] + y));
1386  if (i == step_sizes[0].size())
1387  break;
1388 
1389  x += step_sizes[0][i];
1390  }
1391 
1392  if (j == step_sizes[1].size())
1393  break;
1394 
1395  y += step_sizes[1][j];
1396  }
1397  break;
1398  }
1399  case 3:
1400  {
1401  double z = 0;
1402  for (unsigned int k = 0;; ++k)
1403  {
1404  double y = 0;
1405  for (unsigned int j = 0;; ++j)
1406  {
1407  double x = 0;
1408  for (unsigned int i = 0;; ++i)
1409  {
1410  points.push_back(
1411  Point<dim>(p1[0] + x, p1[1] + y, p1[2] + z));
1412  if (i == step_sizes[0].size())
1413  break;
1414 
1415  x += step_sizes[0][i];
1416  }
1417 
1418  if (j == step_sizes[1].size())
1419  break;
1420 
1421  y += step_sizes[1][j];
1422  }
1423 
1424  if (k == step_sizes[2].size())
1425  break;
1426 
1427  z += step_sizes[2][k];
1428  }
1429  break;
1430  }
1431 
1432  default:
1433  Assert(false, ExcNotImplemented());
1434  }
1435 
1436  // next create the cells
1437  // Prepare cell data
1438  std::vector<CellData<dim>> cells;
1439  switch (dim)
1440  {
1441  case 1:
1442  {
1443  cells.resize(step_sizes[0].size());
1444  for (unsigned int x = 0; x < step_sizes[0].size(); ++x)
1445  {
1446  cells[x].vertices[0] = x;
1447  cells[x].vertices[1] = x + 1;
1448  cells[x].material_id = 0;
1449  }
1450  break;
1451  }
1452 
1453  case 2:
1454  {
1455  cells.resize(step_sizes[1].size() * step_sizes[0].size());
1456  for (unsigned int y = 0; y < step_sizes[1].size(); ++y)
1457  for (unsigned int x = 0; x < step_sizes[0].size(); ++x)
1458  {
1459  const unsigned int c = x + y * step_sizes[0].size();
1460  cells[c].vertices[0] = y * (step_sizes[0].size() + 1) + x;
1461  cells[c].vertices[1] = y * (step_sizes[0].size() + 1) + x + 1;
1462  cells[c].vertices[2] =
1463  (y + 1) * (step_sizes[0].size() + 1) + x;
1464  cells[c].vertices[3] =
1465  (y + 1) * (step_sizes[0].size() + 1) + x + 1;
1466  cells[c].material_id = 0;
1467  }
1468  break;
1469  }
1470 
1471  case 3:
1472  {
1473  const unsigned int n_x = (step_sizes[0].size() + 1);
1474  const unsigned int n_xy =
1475  (step_sizes[0].size() + 1) * (step_sizes[1].size() + 1);
1476 
1477  cells.resize(step_sizes[2].size() * step_sizes[1].size() *
1478  step_sizes[0].size());
1479  for (unsigned int z = 0; z < step_sizes[2].size(); ++z)
1480  for (unsigned int y = 0; y < step_sizes[1].size(); ++y)
1481  for (unsigned int x = 0; x < step_sizes[0].size(); ++x)
1482  {
1483  const unsigned int c =
1484  x + y * step_sizes[0].size() +
1485  z * step_sizes[0].size() * step_sizes[1].size();
1486  cells[c].vertices[0] = z * n_xy + y * n_x + x;
1487  cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
1488  cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
1489  cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
1490  cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
1491  cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
1492  cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
1493  cells[c].vertices[7] =
1494  (z + 1) * n_xy + (y + 1) * n_x + x + 1;
1495  cells[c].material_id = 0;
1496  }
1497  break;
1498  }
1499 
1500  default:
1501  Assert(false, ExcNotImplemented());
1502  }
1503 
1504  tria.create_triangulation(points, cells, SubCellData());
1505 
1506  if (colorize)
1507  {
1508  // to colorize, run through all
1509  // faces of all cells and set
1510  // boundary indicator to the
1511  // correct value if it was 0.
1512 
1513  // use a large epsilon to
1514  // compare numbers to avoid
1515  // roundoff problems.
1516  double min_size =
1517  *std::min_element(step_sizes[0].begin(), step_sizes[0].end());
1518  for (unsigned int i = 1; i < dim; ++i)
1519  min_size = std::min(min_size,
1520  *std::min_element(step_sizes[i].begin(),
1521  step_sizes[i].end()));
1522  const double epsilon = 0.01 * min_size;
1523 
1524  // actual code is external since
1525  // 1-D is different from 2/3D.
1526  colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
1527  }
1528  }
1529 
1530 
1531 
1532  template <>
1533  void
1534  subdivided_hyper_rectangle(Triangulation<1> & tria,
1535  const std::vector<std::vector<double>> &spacing,
1536  const Point<1> & p,
1537  const Table<1, types::material_id> &material_id,
1538  const bool colorize)
1539  {
1540  Assert(spacing.size() == 1, ExcInvalidRepetitionsDimension(1));
1541 
1542  const unsigned int n_cells = material_id.size(0);
1543 
1544  Assert(spacing[0].size() == n_cells, ExcInvalidRepetitionsDimension(1));
1545 
1546  double delta = std::numeric_limits<double>::max();
1547  for (unsigned int i = 0; i < n_cells; i++)
1548  {
1549  Assert(spacing[0][i] >= 0, ExcInvalidRepetitions(-1));
1550  delta = std::min(delta, spacing[0][i]);
1551  }
1552 
1553  // generate the necessary points
1554  std::vector<Point<1>> points;
1555  double ax = p[0];
1556  for (unsigned int x = 0; x <= n_cells; ++x)
1557  {
1558  points.emplace_back(ax);
1559  if (x < n_cells)
1560  ax += spacing[0][x];
1561  }
1562  // create the cells
1563  unsigned int n_val_cells = 0;
1564  for (unsigned int i = 0; i < n_cells; i++)
1565  if (material_id[i] != numbers::invalid_material_id)
1566  n_val_cells++;
1567 
1568  std::vector<CellData<1>> cells(n_val_cells);
1569  unsigned int id = 0;
1570  for (unsigned int x = 0; x < n_cells; ++x)
1571  if (material_id[x] != numbers::invalid_material_id)
1572  {
1573  cells[id].vertices[0] = x;
1574  cells[id].vertices[1] = x + 1;
1575  cells[id].material_id = material_id[x];
1576  id++;
1577  }
1578  // create triangulation
1579  SubCellData t;
1580  GridTools::delete_unused_vertices(points, cells, t);
1581 
1582  tria.create_triangulation(points, cells, t);
1583 
1584  // set boundary indicator
1585  if (colorize)
1586  Assert(false, ExcNotImplemented());
1587  }
1588 
1589 
1590  template <>
1591  void
1592  subdivided_hyper_rectangle(Triangulation<2> & tria,
1593  const std::vector<std::vector<double>> &spacing,
1594  const Point<2> & p,
1595  const Table<2, types::material_id> &material_id,
1596  const bool colorize)
1597  {
1598  Assert(spacing.size() == 2, ExcInvalidRepetitionsDimension(2));
1599 
1600  std::vector<unsigned int> repetitions(2);
1601  unsigned int n_cells = 1;
1602  double delta = std::numeric_limits<double>::max();
1603  for (unsigned int i = 0; i < 2; i++)
1604  {
1605  repetitions[i] = spacing[i].size();
1606  n_cells *= repetitions[i];
1607  for (unsigned int j = 0; j < repetitions[i]; j++)
1608  {
1609  Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
1610  delta = std::min(delta, spacing[i][j]);
1611  }
1612  Assert(material_id.size(i) == repetitions[i],
1613  ExcInvalidRepetitionsDimension(i));
1614  }
1615 
1616  // generate the necessary points
1617  std::vector<Point<2>> points;
1618  double ay = p[1];
1619  for (unsigned int y = 0; y <= repetitions[1]; ++y)
1620  {
1621  double ax = p[0];
1622  for (unsigned int x = 0; x <= repetitions[0]; ++x)
1623  {
1624  points.emplace_back(ax, ay);
1625  if (x < repetitions[0])
1626  ax += spacing[0][x];
1627  }
1628  if (y < repetitions[1])
1629  ay += spacing[1][y];
1630  }
1631 
1632  // create the cells
1633  unsigned int n_val_cells = 0;
1634  for (unsigned int i = 0; i < material_id.size(0); i++)
1635  for (unsigned int j = 0; j < material_id.size(1); j++)
1636  if (material_id[i][j] != numbers::invalid_material_id)
1637  n_val_cells++;
1638 
1639  std::vector<CellData<2>> cells(n_val_cells);
1640  unsigned int id = 0;
1641  for (unsigned int y = 0; y < repetitions[1]; ++y)
1642  for (unsigned int x = 0; x < repetitions[0]; ++x)
1643  if (material_id[x][y] != numbers::invalid_material_id)
1644  {
1645  cells[id].vertices[0] = y * (repetitions[0] + 1) + x;
1646  cells[id].vertices[1] = y * (repetitions[0] + 1) + x + 1;
1647  cells[id].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
1648  cells[id].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
1649  cells[id].material_id = material_id[x][y];
1650  id++;
1651  }
1652 
1653  // create triangulation
1654  SubCellData t;
1655  GridTools::delete_unused_vertices(points, cells, t);
1656 
1657  tria.create_triangulation(points, cells, t);
1658 
1659  // set boundary indicator
1660  if (colorize)
1661  {
1662  double eps = 0.01 * delta;
1663  Triangulation<2>::cell_iterator cell = tria.begin(), endc = tria.end();
1664  for (; cell != endc; ++cell)
1665  {
1666  Point<2> cell_center = cell->center();
1667  for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
1668  if (cell->face(f)->boundary_id() == 0)
1669  {
1670  Point<2> face_center = cell->face(f)->center();
1671  for (unsigned int i = 0; i < 2; ++i)
1672  {
1673  if (face_center[i] < cell_center[i] - eps)
1674  cell->face(f)->set_boundary_id(i * 2);
1675  if (face_center[i] > cell_center[i] + eps)
1676  cell->face(f)->set_boundary_id(i * 2 + 1);
1677  }
1678  }
1679  }
1680  }
1681  }
1682 
1683 
1684  template <>
1685  void
1686  subdivided_hyper_rectangle(Triangulation<3> & tria,
1687  const std::vector<std::vector<double>> &spacing,
1688  const Point<3> & p,
1689  const Table<3, types::material_id> &material_id,
1690  const bool colorize)
1691  {
1692  const unsigned int dim = 3;
1693 
1694  Assert(spacing.size() == dim, ExcInvalidRepetitionsDimension(dim));
1695 
1696  std::vector<unsigned int> repetitions(dim);
1697  unsigned int n_cells = 1;
1698  double delta = std::numeric_limits<double>::max();
1699  for (unsigned int i = 0; i < dim; i++)
1700  {
1701  repetitions[i] = spacing[i].size();
1702  n_cells *= repetitions[i];
1703  for (unsigned int j = 0; j < repetitions[i]; j++)
1704  {
1705  Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
1706  delta = std::min(delta, spacing[i][j]);
1707  }
1708  Assert(material_id.size(i) == repetitions[i],
1709  ExcInvalidRepetitionsDimension(i));
1710  }
1711 
1712  // generate the necessary points
1713  std::vector<Point<dim>> points;
1714  double az = p[2];
1715  for (unsigned int z = 0; z <= repetitions[2]; ++z)
1716  {
1717  double ay = p[1];
1718  for (unsigned int y = 0; y <= repetitions[1]; ++y)
1719  {
1720  double ax = p[0];
1721  for (unsigned int x = 0; x <= repetitions[0]; ++x)
1722  {
1723  points.emplace_back(ax, ay, az);
1724  if (x < repetitions[0])
1725  ax += spacing[0][x];
1726  }
1727  if (y < repetitions[1])
1728  ay += spacing[1][y];
1729  }
1730  if (z < repetitions[2])
1731  az += spacing[2][z];
1732  }
1733 
1734  // create the cells
1735  unsigned int n_val_cells = 0;
1736  for (unsigned int i = 0; i < material_id.size(0); i++)
1737  for (unsigned int j = 0; j < material_id.size(1); j++)
1738  for (unsigned int k = 0; k < material_id.size(2); k++)
1739  if (material_id[i][j][k] != numbers::invalid_material_id)
1740  n_val_cells++;
1741 
1742  std::vector<CellData<dim>> cells(n_val_cells);
1743  unsigned int id = 0;
1744  const unsigned int n_x = (repetitions[0] + 1);
1745  const unsigned int n_xy = (repetitions[0] + 1) * (repetitions[1] + 1);
1746  for (unsigned int z = 0; z < repetitions[2]; ++z)
1747  for (unsigned int y = 0; y < repetitions[1]; ++y)
1748  for (unsigned int x = 0; x < repetitions[0]; ++x)
1749  if (material_id[x][y][z] != numbers::invalid_material_id)
1750  {
1751  cells[id].vertices[0] = z * n_xy + y * n_x + x;
1752  cells[id].vertices[1] = z * n_xy + y * n_x + x + 1;
1753  cells[id].vertices[2] = z * n_xy + (y + 1) * n_x + x;
1754  cells[id].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
1755  cells[id].vertices[4] = (z + 1) * n_xy + y * n_x + x;
1756  cells[id].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
1757  cells[id].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
1758  cells[id].vertices[7] = (z + 1) * n_xy + (y + 1) * n_x + x + 1;
1759  cells[id].material_id = material_id[x][y][z];
1760  id++;
1761  }
1762 
1763  // create triangulation
1764  SubCellData t;
1765  GridTools::delete_unused_vertices(points, cells, t);
1766 
1767  tria.create_triangulation(points, cells, t);
1768 
1769  // set boundary indicator
1770  if (colorize)
1771  {
1772  double eps = 0.01 * delta;
1773  Triangulation<dim>::cell_iterator cell = tria.begin(),
1774  endc = tria.end();
1775  for (; cell != endc; ++cell)
1776  {
1777  Point<dim> cell_center = cell->center();
1778  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
1779  if (cell->face(f)->boundary_id() == 0)
1780  {
1781  Point<dim> face_center = cell->face(f)->center();
1782  for (unsigned int i = 0; i < dim; ++i)
1783  {
1784  if (face_center[i] < cell_center[i] - eps)
1785  cell->face(f)->set_boundary_id(i * 2);
1786  if (face_center[i] > cell_center[i] + eps)
1787  cell->face(f)->set_boundary_id(i * 2 + 1);
1788  }
1789  }
1790  }
1791  }
1792  }
1793 
1794  template <int dim, int spacedim>
1795  void
1796  cheese(Triangulation<dim, spacedim> & tria,
1797  const std::vector<unsigned int> &holes)
1798  {
1799  AssertDimension(holes.size(), dim);
1800  // The corner points of the first cell. If there is a desire at
1801  // some point to change the geometry of the cells, they can be
1802  // made an argument to the function.
1803 
1804  Point<spacedim> p1;
1805  Point<spacedim> p2;
1806  for (unsigned int d = 0; d < dim; ++d)
1807  p2(d) = 1.;
1808 
1809  // then check that all repetitions
1810  // are >= 1, and calculate deltas
1811  // convert repetitions from double
1812  // to int by taking the ceiling.
1813  std::vector<Point<spacedim>> delta(dim);
1814  unsigned int repetitions[dim];
1815  for (unsigned int i = 0; i < dim; ++i)
1816  {
1817  Assert(holes[i] >= 1,
1818  ExcMessage("At least one hole needed in each direction"));
1819  repetitions[i] = 2 * holes[i] + 1;
1820  delta[i][i] = (p2[i] - p1[i]);
1821  }
1822 
1823  // then generate the necessary
1824  // points
1825  std::vector<Point<spacedim>> points;
1826  switch (dim)
1827  {
1828  case 1:
1829  for (unsigned int x = 0; x <= repetitions[0]; ++x)
1830  points.push_back(p1 + (double)x * delta[0]);
1831  break;
1832 
1833  case 2:
1834  for (unsigned int y = 0; y <= repetitions[1]; ++y)
1835  for (unsigned int x = 0; x <= repetitions[0]; ++x)
1836  points.push_back(p1 + (double)x * delta[0] +
1837  (double)y * delta[1]);
1838  break;
1839 
1840  case 3:
1841  for (unsigned int z = 0; z <= repetitions[2]; ++z)
1842  for (unsigned int y = 0; y <= repetitions[1]; ++y)
1843  for (unsigned int x = 0; x <= repetitions[0]; ++x)
1844  points.push_back(p1 + (double)x * delta[0] +
1845  (double)y * delta[1] + (double)z * delta[2]);
1846  break;
1847 
1848  default:
1849  Assert(false, ExcNotImplemented());
1850  }
1851 
1852  // next create the cells
1853  // Prepare cell data
1854  std::vector<CellData<dim>> cells;
1855  switch (dim)
1856  {
1857  case 2:
1858  {
1859  cells.resize(repetitions[1] * repetitions[0] - holes[1] * holes[0]);
1860  unsigned int c = 0;
1861  for (unsigned int y = 0; y < repetitions[1]; ++y)
1862  for (unsigned int x = 0; x < repetitions[0]; ++x)
1863  {
1864  if ((x % 2 == 1) && (y % 2 == 1))
1865  continue;
1866  Assert(c < cells.size(), ExcInternalError());
1867  cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
1868  cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
1869  cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
1870  cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
1871  cells[c].material_id = 0;
1872  ++c;
1873  }
1874  break;
1875  }
1876 
1877  case 3:
1878  {
1879  const unsigned int n_x = (repetitions[0] + 1);
1880  const unsigned int n_xy =
1881  (repetitions[0] + 1) * (repetitions[1] + 1);
1882 
1883  cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
1884 
1885  unsigned int c = 0;
1886  for (unsigned int z = 0; z < repetitions[2]; ++z)
1887  for (unsigned int y = 0; y < repetitions[1]; ++y)
1888  for (unsigned int x = 0; x < repetitions[0]; ++x)
1889  {
1890  Assert(c < cells.size(), ExcInternalError());
1891  cells[c].vertices[0] = z * n_xy + y * n_x + x;
1892  cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
1893  cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
1894  cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
1895  cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
1896  cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
1897  cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
1898  cells[c].vertices[7] =
1899  (z + 1) * n_xy + (y + 1) * n_x + x + 1;
1900  cells[c].material_id = 0;
1901  ++c;
1902  }
1903  break;
1904  }
1905 
1906  default:
1907  Assert(false, ExcNotImplemented());
1908  }
1909 
1910  tria.create_triangulation(points, cells, SubCellData());
1911  }
1912 
1913 
1914 
1915  template <>
1916  void plate_with_a_hole(Triangulation<1> & /*tria*/,
1917  const double /*inner_radius*/,
1918  const double /*outer_radius*/,
1919  const double /*pad_bottom*/,
1920  const double /*pad_top*/,
1921  const double /*pad_left*/,
1922  const double /*pad_right*/,
1923  const Point<1> /*center*/,
1924  const types::manifold_id /*polar_manifold_id*/,
1925  const types::manifold_id /*tfi_manifold_id*/,
1926  const double /*L*/,
1927  const unsigned int /*n_slices*/,
1928  const bool /*colorize*/)
1929  {
1930  Assert(false, ExcNotImplemented());
1931  }
1932 
1933 
1934  namespace internal
1935  {
1936  // helper function to check if point is in 2D box
1937  bool inline point_in_2d_box(const Point<2> &p,
1938  const Point<2> &c,
1939  const double radius)
1940  {
1941  return (std::abs(p[0] - c[0]) < radius) &&
1942  (std::abs(p[1] - c[1]) < radius);
1943  }
1944 
1945  // same as above but will ingore the third component
1946  // of both points
1947  bool inline point_in_2d_box(const Point<3> &p,
1948  const Point<3> &center,
1949  const double radius)
1950  {
1951  return point_in_2d_box(Point<2>(p[0], p[1]),
1952  Point<2>(center[0], center[1]),
1953  radius);
1954  }
1955 
1961  void plate_with_a_hole(Triangulation<2> & tria,
1962  const double inner_radius,
1963  const double outer_radius,
1964  const double pad_bottom,
1965  const double pad_top,
1966  const double pad_left,
1967  const double pad_right,
1968  const Point<2> new_center,
1969  const types::manifold_id polar_manifold_id,
1970  const types::manifold_id tfi_manifold_id,
1971  const double L,
1972  const unsigned int /*n_slices*/,
1973  const bool colorize)
1974  {
1975  Assert((pad_bottom > 0 || pad_top > 0 || pad_left > 0 || pad_right > 0),
1976  ExcMessage("At least one padding parameter has to be non-zero."));
1977  Assert(pad_bottom >= 0., ExcMessage("Negative bottom padding."));
1978  Assert(pad_top >= 0., ExcMessage("Negative top padding."));
1979  Assert(pad_left >= 0., ExcMessage("Negative left padding."));
1980  Assert(pad_right >= 0., ExcMessage("Negative right padding."));
1981 
1982  const Point<2> center;
1983 
1984  auto min_line_length = [](const Triangulation<2> &tria) -> double {
1985  double length = std::numeric_limits<double>::max();
1986  for (const auto &cell : tria.active_cell_iterators())
1987  for (unsigned int n = 0; n < GeometryInfo<2>::lines_per_cell; ++n)
1988  length = std::min(length, cell->line(n)->diameter());
1989  return length;
1990  };
1991 
1992  // start by setting up the cylinder triangulation
1993  Triangulation<2> cylinder_tria;
1994  GridGenerator::hyper_cube_with_cylindrical_hole(cylinder_tria,
1995  inner_radius,
1996  outer_radius,
1997  L,
1998  /*repetitions*/ 1,
1999  colorize);
2000 
2001  // Mark the cylinder material ids as 2 to distinguish them from the bulk
2002  // cells (these are read again in the loop that sets manifold ids)
2003  for (const auto &cell : cylinder_tria.active_cell_iterators())
2004  cell->set_material_id(2);
2005 
2006  // hyper_cube_with_cylindrical_hole will have 2 cells along
2007  // each face, so he element size is outer_radius
2008 
2009  auto add_sizes = [](std::vector<double> &step_sizes,
2010  const double padding,
2011  const double h) -> void {
2012  // use std::round instead of std::ceil to improve aspect ratio
2013  // in case padding is only slightly larger than h.
2014  const unsigned int rounded = std::round(padding / h);
2015  // in case padding is much smaller than h, make sure we
2016  // have at least 1 element
2017  const unsigned int num = (padding > 0. && rounded == 0) ? 1 : rounded;
2018  for (unsigned int i = 0; i < num; ++i)
2019  step_sizes.push_back(padding / num);
2020  };
2021 
2022  std::vector<std::vector<double>> step_sizes(2);
2023  // x-coord
2024  // left:
2025  add_sizes(step_sizes[0], pad_left, outer_radius);
2026  // center
2027  step_sizes[0].push_back(outer_radius);
2028  step_sizes[0].push_back(outer_radius);
2029  // right
2030  add_sizes(step_sizes[0], pad_right, outer_radius);
2031  // y-coord
2032  // bottom
2033  add_sizes(step_sizes[1], pad_bottom, outer_radius);
2034  // center
2035  step_sizes[1].push_back(outer_radius);
2036  step_sizes[1].push_back(outer_radius);
2037  // top
2038  add_sizes(step_sizes[1], pad_top, outer_radius);
2039 
2040  // now create bulk
2041  Triangulation<2> bulk_tria;
2042  const Point<2> bl(-outer_radius - pad_left, -outer_radius - pad_bottom);
2043  const Point<2> tr(outer_radius + pad_right, outer_radius + pad_top);
2044  GridGenerator::subdivided_hyper_rectangle(
2045  bulk_tria, step_sizes, bl, tr, colorize);
2046 
2047  // now remove cells reserved from the cylindrical hole
2048  std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
2049  for (const auto &cell : bulk_tria.active_cell_iterators())
2050  if (point_in_2d_box(cell->center(), center, outer_radius))
2051  cells_to_remove.insert(cell);
2052 
2053  Triangulation<2> tria_without_cylinder;
2054  GridGenerator::create_triangulation_with_removed_cells(
2055  bulk_tria, cells_to_remove, tria_without_cylinder);
2056 
2057  const double tolerance = std::min(min_line_length(tria_without_cylinder),
2058  min_line_length(cylinder_tria)) /
2059  2.0;
2060 
2061  // set material ID in bulk part to something other than 2
2062  for (const auto &cell : tria_without_cylinder.active_cell_iterators())
2063  cell->set_material_id(1);
2064 
2065  GridGenerator::merge_triangulations(tria_without_cylinder,
2066  cylinder_tria,
2067  tria,
2068  tolerance);
2069 
2070  // now set manifold ids:
2071  for (const auto &cell : tria.active_cell_iterators())
2072  {
2073  // set all non-boundary manifold ids on the cells that came from the
2074  // grid around the cylinder to the new TFI manifold id.
2075  if (cell->material_id() == 2)
2076  {
2077  cell->set_manifold_id(tfi_manifold_id);
2078  for (unsigned int face_n = 0;
2079  face_n < GeometryInfo<2>::faces_per_cell;
2080  ++face_n)
2081  {
2082  const auto &face = cell->face(face_n);
2083  if (face->at_boundary() &&
2084  point_in_2d_box(face->center(), center, outer_radius))
2085  face->set_manifold_id(polar_manifold_id);
2086  else
2087  face->set_manifold_id(tfi_manifold_id);
2088  }
2089  }
2090  else
2091  {
2092  // ensure that all other manifold ids (including the faces
2093  // opposite the cylinder) are set to the flat id
2094  cell->set_all_manifold_ids(numbers::flat_manifold_id);
2095  }
2096  }
2097 
2098  static constexpr double tol =
2099  std::numeric_limits<double>::epsilon() * 10000;
2100  if (colorize)
2101  for (auto &cell : tria.active_cell_iterators())
2102  for (unsigned int face_n = 0;
2103  face_n < GeometryInfo<2>::faces_per_cell;
2104  ++face_n)
2105  {
2106  const auto face = cell->face(face_n);
2107  if (face->at_boundary())
2108  {
2109  const Point<2> center = face->center();
2110  // left side
2111  if (std::abs(center[0] - bl[0]) < tol * std::abs(bl[0]))
2112  face->set_boundary_id(0);
2113  // right side
2114  else if (std::abs(center[0] - tr[0]) < tol * std::abs(tr[0]))
2115  face->set_boundary_id(1);
2116  // bottom
2117  else if (std::abs(center[1] - bl[1]) < tol * std::abs(bl[1]))
2118  face->set_boundary_id(2);
2119  // top
2120  else if (std::abs(center[1] - tr[1]) < tol * std::abs(tr[1]))
2121  face->set_boundary_id(3);
2122  // cylinder boundary
2123  else
2124  {
2125  Assert(cell->material_id() == 2, ExcInternalError());
2126  face->set_boundary_id(4);
2127  }
2128  }
2129  }
2130 
2131  // move to the new center
2132  GridTools::shift(new_center, tria);
2133 
2134  PolarManifold<2> polar_manifold(new_center);
2135  tria.set_manifold(polar_manifold_id, polar_manifold);
2136  TransfiniteInterpolationManifold<2> inner_manifold;
2137  inner_manifold.initialize(tria);
2138  tria.set_manifold(tfi_manifold_id, inner_manifold);
2139  }
2140 
2141  } // namespace internal
2142 
2143  template <>
2144  void plate_with_a_hole(Triangulation<2> & tria,
2145  const double inner_radius,
2146  const double outer_radius,
2147  const double pad_bottom,
2148  const double pad_top,
2149  const double pad_left,
2150  const double pad_right,
2151  const Point<2> new_center,
2152  const types::manifold_id polar_manifold_id,
2153  const types::manifold_id tfi_manifold_id,
2154  const double L,
2155  const unsigned int n_slices,
2156  const bool colorize)
2157  {
2158  internal::plate_with_a_hole(tria,
2159  inner_radius,
2160  outer_radius,
2161  pad_bottom,
2162  pad_top,
2163  pad_left,
2164  pad_right,
2165  new_center,
2166  polar_manifold_id,
2167  tfi_manifold_id,
2168  L,
2169  n_slices,
2170  colorize);
2171 
2172  // reset material id
2173  for (auto &cell : tria.active_cell_iterators())
2174  cell->set_material_id(0);
2175  }
2176 
2177 
2178 
2179  template <>
2180  void plate_with_a_hole(Triangulation<3> & tria,
2181  const double inner_radius,
2182  const double outer_radius,
2183  const double pad_bottom,
2184  const double pad_top,
2185  const double pad_left,
2186  const double pad_right,
2187  const Point<3> new_center,
2188  const types::manifold_id polar_manifold_id,
2189  const types::manifold_id tfi_manifold_id,
2190  const double L,
2191  const unsigned int n_slices,
2192  const bool colorize)
2193  {
2194  Triangulation<2> tria_2;
2195  internal::plate_with_a_hole(tria_2,
2196  inner_radius,
2197  outer_radius,
2198  pad_bottom,
2199  pad_top,
2200  pad_left,
2201  pad_right,
2202  Point<2>(new_center[0], new_center[1]),
2203  polar_manifold_id,
2204  tfi_manifold_id,
2205  L,
2206  n_slices,
2207  colorize);
2208 
2209  // extrude to 3D
2210  extrude_triangulation(tria_2, n_slices, L, tria);
2211 
2212  // shift in Z direction to match specified center
2213  GridTools::shift(Point<3>(0, 0, new_center[2] - L / 2.), tria);
2214 
2215  // extrude will keep boundary IDs but can not keep manifolds, set them
2216  // below:
2217  const FlatManifold<3> flat_manifold;
2218  for (const auto &cell : tria.active_cell_iterators())
2219  if (cell->material_id() == 2)
2220  {
2221  cell->set_all_manifold_ids(tfi_manifold_id);
2222  for (unsigned int face_n = 0;
2223  face_n < GeometryInfo<3>::faces_per_cell;
2224  ++face_n)
2225  {
2226  // similar check to 2D version
2227  const auto face = cell->face(face_n);
2228  if (face->at_boundary() &&
2229  internal::point_in_2d_box(face->center(),
2230  new_center,
2231  outer_radius) &&
2232  std::abs(
2233  flat_manifold.normal_vector(face, face->center())[2]) <
2234  1.0e-10)
2235  face->set_manifold_id(polar_manifold_id);
2236  else
2237  face->set_manifold_id(tfi_manifold_id);
2238  }
2239  }
2240  else
2241  {
2242  // otherwise set everything to flat manifold
2243  cell->set_all_manifold_ids(numbers::flat_manifold_id);
2244  }
2245 
2246  // set up the new manifolds
2247  const Tensor<1, 3> direction{{0.0, 0.0, 1.0}};
2248  const CylindricalManifold<3> cylindrical_manifold(
2249  direction, /*axial_point*/ new_center);
2250  TransfiniteInterpolationManifold<3> inner_manifold;
2251  inner_manifold.initialize(tria);
2252  tria.set_manifold(polar_manifold_id, cylindrical_manifold);
2253  tria.set_manifold(tfi_manifold_id, inner_manifold);
2254 
2255  // finally reset material ids
2256  for (auto &cell : tria.active_cell_iterators())
2257  cell->set_material_id(0);
2258  }
2259 
2260 
2261 
2262  template <int dim, int spacedim>
2263  void
2264  hyper_cross(Triangulation<dim, spacedim> & tria,
2265  const std::vector<unsigned int> &sizes,
2266  const bool colorize)
2267  {
2268  AssertDimension(sizes.size(), GeometryInfo<dim>::faces_per_cell);
2269  Assert(dim > 1, ExcNotImplemented());
2270  Assert(dim < 4, ExcNotImplemented());
2271 
2272  // If there is a desire at some point to change the geometry of
2273  // the cells, this tensor can be made an argument to the function.
2274  Tensor<1, dim> dimensions;
2275  for (unsigned int d = 0; d < dim; ++d)
2276  dimensions[d] = 1.;
2277 
2278  std::vector<Point<spacedim>> points;
2279  unsigned int n_cells = 1;
2280  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
2281  n_cells += sizes[i];
2282 
2283  std::vector<CellData<dim>> cells(n_cells);
2284  // Vertices of the center cell
2285  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
2286  {
2287  Point<spacedim> p;
2288  for (unsigned int d = 0; d < dim; ++d)
2289  p(d) = 0.5 * dimensions[d] *
2290  GeometryInfo<dim>::unit_normal_orientation
2291  [GeometryInfo<dim>::vertex_to_face[i][d]];
2292  points.push_back(p);
2293  cells[0].vertices[i] = i;
2294  }
2295  cells[0].material_id = 0;
2296 
2297  // The index of the first cell of the leg.
2298  unsigned int cell_index = 1;
2299  // The legs of the cross
2300  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2301  ++face)
2302  {
2303  const unsigned int oface = GeometryInfo<dim>::opposite_face[face];
2304  const unsigned int dir = GeometryInfo<dim>::unit_normal_direction[face];
2305 
2306  // We are moving in the direction of face
2307  for (unsigned int j = 0; j < sizes[face]; ++j, ++cell_index)
2308  {
2309  const unsigned int last_cell = (j == 0) ? 0U : (cell_index - 1);
2310 
2311  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
2312  ++v)
2313  {
2314  const unsigned int cellv =
2315  GeometryInfo<dim>::face_to_cell_vertices(face, v);
2316  const unsigned int ocellv =
2317  GeometryInfo<dim>::face_to_cell_vertices(oface, v);
2318  // First the vertices which already exist
2319  cells[cell_index].vertices[ocellv] =
2320  cells[last_cell].vertices[cellv];
2321 
2322  // Now the new vertices
2323  cells[cell_index].vertices[cellv] = points.size();
2324 
2325  Point<spacedim> p = points[cells[cell_index].vertices[ocellv]];
2326  p(dir) += GeometryInfo<dim>::unit_normal_orientation[face] *
2327  dimensions[dir];
2328  points.push_back(p);
2329  }
2330  cells[cell_index].material_id = (colorize) ? (face + 1U) : 0U;
2331  }
2332  }
2333  tria.create_triangulation(points, cells, SubCellData());
2334  }
2335 
2336 
2337  template <>
2338  void
2339  hyper_cube_slit(Triangulation<1> &, const double, const double, const bool)
2340  {
2341  Assert(false, ExcNotImplemented());
2342  }
2343 
2344 
2345 
2346  template <>
2347  void enclosed_hyper_cube(Triangulation<1> &,
2348  const double,
2349  const double,
2350  const double,
2351  const bool)
2352  {
2353  Assert(false, ExcNotImplemented());
2354  }
2355 
2356 
2357 
2358  template <>
2359  void hyper_L(Triangulation<1> &, const double, const double, const bool)
2360  {
2361  Assert(false, ExcNotImplemented());
2362  }
2363 
2364 
2365 
2366  template <>
2367  void
2368  hyper_ball(Triangulation<1> &, const Point<1> &, const double, const bool)
2369  {
2370  Assert(false, ExcNotImplemented());
2371  }
2372 
2373 
2374 
2375  template <>
2376  void cylinder(Triangulation<1> &, const double, const double)
2377  {
2378  Assert(false, ExcNotImplemented());
2379  }
2380 
2381 
2382 
2383  template <>
2384  void
2385  truncated_cone(Triangulation<1> &, const double, const double, const double)
2386  {
2387  Assert(false, ExcNotImplemented());
2388  }
2389 
2390 
2391 
2392  template <>
2393  void hyper_shell(Triangulation<1> &,
2394  const Point<1> &,
2395  const double,
2396  const double,
2397  const unsigned int,
2398  const bool)
2399  {
2400  Assert(false, ExcNotImplemented());
2401  }
2402 
2403 
2404  template <>
2405  void cylinder_shell(Triangulation<1> &,
2406  const double,
2407  const double,
2408  const double,
2409  const unsigned int,
2410  const unsigned int)
2411  {
2412  Assert(false, ExcNotImplemented());
2413  }
2414 
2415 
2416  template <>
2417  void quarter_hyper_ball(Triangulation<1> &, const Point<1> &, const double)
2418  {
2419  Assert(false, ExcNotImplemented());
2420  }
2421 
2422 
2423  template <>
2424  void half_hyper_ball(Triangulation<1> &, const Point<1> &, const double)
2425  {
2426  Assert(false, ExcNotImplemented());
2427  }
2428 
2429 
2430  template <>
2431  void half_hyper_shell(Triangulation<1> &,
2432  const Point<1> &,
2433  const double,
2434  const double,
2435  const unsigned int,
2436  const bool)
2437  {
2438  Assert(false, ExcNotImplemented());
2439  }
2440 
2441  template <>
2442  void quarter_hyper_shell(Triangulation<1> &,
2443  const Point<1> &,
2444  const double,
2445  const double,
2446  const unsigned int,
2447  const bool)
2448  {
2449  Assert(false, ExcNotImplemented());
2450  }
2451 
2452  template <>
2453  void enclosed_hyper_cube(Triangulation<2> &tria,
2454  const double left,
2455  const double right,
2456  const double thickness,
2457  const bool colorize)
2458  {
2459  Assert(left < right,
2460  ExcMessage("Invalid left-to-right bounds of enclosed hypercube"));
2461 
2462  std::vector<Point<2>> vertices(16);
2463  double coords[4];
2464  coords[0] = left - thickness;
2465  coords[1] = left;
2466  coords[2] = right;
2467  coords[3] = right + thickness;
2468 
2469  unsigned int k = 0;
2470  for (unsigned int i0 = 0; i0 < 4; ++i0)
2471  for (unsigned int i1 = 0; i1 < 4; ++i1)
2472  vertices[k++] = Point<2>(coords[i1], coords[i0]);
2473 
2474  const types::material_id materials[9] = {5, 4, 6, 1, 0, 2, 9, 8, 10};
2475 
2476  std::vector<CellData<2>> cells(9);
2477  k = 0;
2478  for (unsigned int i0 = 0; i0 < 3; ++i0)
2479  for (unsigned int i1 = 0; i1 < 3; ++i1)
2480  {
2481  cells[k].vertices[0] = i1 + 4 * i0;
2482  cells[k].vertices[1] = i1 + 4 * i0 + 1;
2483  cells[k].vertices[2] = i1 + 4 * i0 + 4;
2484  cells[k].vertices[3] = i1 + 4 * i0 + 5;
2485  if (colorize)
2486  cells[k].material_id = materials[k];
2487  ++k;
2488  }
2489  tria.create_triangulation(vertices,
2490  cells,
2491  SubCellData()); // no boundary information
2492  }
2493 
2494 
2495 
2496  // Implementation for 2D only
2497  template <>
2498  void hyper_cube_slit(Triangulation<2> &tria,
2499  const double left,
2500  const double right,
2501  const bool colorize)
2502  {
2503  const double rl2 = (right + left) / 2;
2504  const Point<2> vertices[10] = {Point<2>(left, left),
2505  Point<2>(rl2, left),
2506  Point<2>(rl2, rl2),
2507  Point<2>(left, rl2),
2508  Point<2>(right, left),
2509  Point<2>(right, rl2),
2510  Point<2>(rl2, right),
2511  Point<2>(left, right),
2512  Point<2>(right, right),
2513  Point<2>(rl2, left)};
2514  const int cell_vertices[4][4] = {{0, 1, 3, 2},
2515  {9, 4, 2, 5},
2516  {3, 2, 7, 6},
2517  {2, 5, 6, 8}};
2518  std::vector<CellData<2>> cells(4, CellData<2>());
2519  for (unsigned int i = 0; i < 4; ++i)
2520  {
2521  for (unsigned int j = 0; j < 4; ++j)
2522  cells[i].vertices[j] = cell_vertices[i][j];
2523  cells[i].material_id = 0;
2524  }
2525  tria.create_triangulation(std::vector<Point<2>>(std::begin(vertices),
2526  std::end(vertices)),
2527  cells,
2528  SubCellData()); // no boundary information
2529 
2530  if (colorize)
2531  {
2532  Triangulation<2>::cell_iterator cell = tria.begin();
2533  cell->face(1)->set_boundary_id(1);
2534  ++cell;
2535  cell->face(0)->set_boundary_id(2);
2536  }
2537  }
2538 
2539 
2540 
2541  template <>
2542  void truncated_cone(Triangulation<2> &triangulation,
2543  const double radius_0,
2544  const double radius_1,
2545  const double half_length)
2546  {
2547  Point<2> vertices_tmp[4];
2548 
2549  vertices_tmp[0] = Point<2>(-half_length, -radius_0);
2550  vertices_tmp[1] = Point<2>(half_length, -radius_1);
2551  vertices_tmp[2] = Point<2>(-half_length, radius_0);
2552  vertices_tmp[3] = Point<2>(half_length, radius_1);
2553 
2554  const std::vector<Point<2>> vertices(std::begin(vertices_tmp),
2555  std::end(vertices_tmp));
2556  unsigned int cell_vertices[1][GeometryInfo<2>::vertices_per_cell];
2557 
2558  for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
2559  cell_vertices[0][i] = i;
2560 
2561  std::vector<CellData<2>> cells(1, CellData<2>());
2562 
2563  for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
2564  cells[0].vertices[i] = cell_vertices[0][i];
2565 
2566  cells[0].material_id = 0;
2567  triangulation.create_triangulation(vertices, cells, SubCellData());
2568 
2569  Triangulation<2>::cell_iterator cell = triangulation.begin();
2570 
2571  cell->face(0)->set_boundary_id(1);
2572  cell->face(1)->set_boundary_id(2);
2573 
2574  for (unsigned int i = 2; i < 4; ++i)
2575  cell->face(i)->set_boundary_id(0);
2576  }
2577 
2578 
2579 
2580  // Implementation for 2D only
2581  template <>
2582  void hyper_L(Triangulation<2> &tria,
2583  const double a,
2584  const double b,
2585  const bool colorize)
2586  {
2587  const Point<2> vertices[8] = {Point<2>(a, a),
2588  Point<2>((a + b) / 2, a),
2589  Point<2>(b, a),
2590  Point<2>(a, (a + b) / 2),
2591  Point<2>((a + b) / 2, (a + b) / 2),
2592  Point<2>(b, (a + b) / 2),
2593  Point<2>(a, b),
2594  Point<2>((a + b) / 2, b)};
2595  const int cell_vertices[3][4] = {{0, 1, 3, 4}, {1, 2, 4, 5}, {3, 4, 6, 7}};
2596 
2597  std::vector<CellData<2>> cells(3, CellData<2>());
2598 
2599  for (unsigned int i = 0; i < 3; ++i)
2600  {
2601  for (unsigned int j = 0; j < 4; ++j)
2602  cells[i].vertices[j] = cell_vertices[i][j];
2603  cells[i].material_id = 0;
2604  }
2605 
2606  tria.create_triangulation(std::vector<Point<2>>(std::begin(vertices),
2607  std::end(vertices)),
2608  cells,
2609  SubCellData());
2610 
2611  if (colorize)
2612  {
2613  Triangulation<2>::cell_iterator cell = tria.begin();
2614 
2615  cell->face(0)->set_boundary_id(0);
2616  cell->face(2)->set_boundary_id(1);
2617  cell++;
2618 
2619  cell->face(1)->set_boundary_id(2);
2620  cell->face(2)->set_boundary_id(1);
2621  cell->face(3)->set_boundary_id(3);
2622  cell++;
2623 
2624  cell->face(0)->set_boundary_id(0);
2625  cell->face(1)->set_boundary_id(4);
2626  cell->face(3)->set_boundary_id(5);
2627  }
2628  }
2629 
2630 
2631 
2632  // Implementation for 2D only
2633  template <>
2634  void hyper_ball(Triangulation<2> &tria,
2635  const Point<2> & p,
2636  const double radius,
2637  const bool internal_manifolds)
2638  {
2639  // equilibrate cell sizes at
2640  // transition from the inner part
2641  // to the radial cells
2642  const double a = 1. / (1 + std::sqrt(2.0));
2643  const Point<2> vertices[8] = {
2644  p + Point<2>(-1, -1) * (radius / std::sqrt(2.0)),
2645  p + Point<2>(+1, -1) * (radius / std::sqrt(2.0)),
2646  p + Point<2>(-1, -1) * (radius / std::sqrt(2.0) * a),
2647  p + Point<2>(+1, -1) * (radius / std::sqrt(2.0) * a),
2648  p + Point<2>(-1, +1) * (radius / std::sqrt(2.0) * a),
2649  p + Point<2>(+1, +1) * (radius / std::sqrt(2.0) * a),
2650  p + Point<2>(-1, +1) * (radius / std::sqrt(2.0)),
2651  p + Point<2>(+1, +1) * (radius / std::sqrt(2.0))};
2652 
2653  const int cell_vertices[5][4] = {
2654  {0, 1, 2, 3}, {0, 2, 6, 4}, {2, 3, 4, 5}, {1, 7, 3, 5}, {6, 4, 7, 5}};
2655 
2656  std::vector<CellData<2>> cells(5, CellData<2>());
2657 
2658  for (unsigned int i = 0; i < 5; ++i)
2659  {
2660  for (unsigned int j = 0; j < 4; ++j)
2661  cells[i].vertices[j] = cell_vertices[i][j];
2662  cells[i].material_id = 0;
2663  cells[i].manifold_id = i == 2 ? 1 : numbers::flat_manifold_id;
2664  }
2665 
2666  tria.create_triangulation(std::vector<Point<2>>(std::begin(vertices),
2667  std::end(vertices)),
2668  cells,
2669  SubCellData()); // no boundary information
2670  tria.set_all_manifold_ids_on_boundary(0);
2671  tria.set_manifold(0, SphericalManifold<2>(p));
2672  if (internal_manifolds)
2673  tria.set_manifold(1, SphericalManifold<2>(p));
2674  }
2675 
2676 
2677 
2678  template <>
2679  void hyper_shell(Triangulation<2> & tria,
2680  const Point<2> & center,
2681  const double inner_radius,
2682  const double outer_radius,
2683  const unsigned int n_cells,
2684  const bool colorize)
2685  {
2686  Assert((inner_radius > 0) && (inner_radius < outer_radius),
2687  ExcInvalidRadii());
2688 
2689  const double pi = numbers::PI;
2690 
2691  // determine the number of cells
2692  // for the grid. if not provided by
2693  // the user determine it such that
2694  // the length of each cell on the
2695  // median (in the middle between
2696  // the two circles) is equal to its
2697  // radial extent (which is the
2698  // difference between the two
2699  // radii)
2700  const unsigned int N =
2701  (n_cells == 0 ? static_cast<unsigned int>(
2702  std::ceil((2 * pi * (outer_radius + inner_radius) / 2) /
2703  (outer_radius - inner_radius))) :
2704  n_cells);
2705 
2706  // set up N vertices on the
2707  // outer and N vertices on
2708  // the inner circle. the
2709  // first N ones are on the
2710  // outer one, and all are
2711  // numbered counter-clockwise
2712  std::vector<Point<2>> vertices(2 * N);
2713  for (unsigned int i = 0; i < N; ++i)
2714  {
2715  vertices[i] =
2716  Point<2>(std::cos(2 * pi * i / N), std::sin(2 * pi * i / N)) *
2717  outer_radius;
2718  vertices[i + N] = vertices[i] * (inner_radius / outer_radius);
2719 
2720  vertices[i] += center;
2721  vertices[i + N] += center;
2722  }
2723 
2724  std::vector<CellData<2>> cells(N, CellData<2>());
2725 
2726  for (unsigned int i = 0; i < N; ++i)
2727  {
2728  cells[i].vertices[0] = i;
2729  cells[i].vertices[1] = (i + 1) % N;
2730  cells[i].vertices[2] = N + i;
2731  cells[i].vertices[3] = N + ((i + 1) % N);
2732 
2733  cells[i].material_id = 0;
2734  }
2735 
2736  tria.create_triangulation(vertices, cells, SubCellData());
2737 
2738  if (colorize)
2739  colorize_hyper_shell(tria, center, inner_radius, outer_radius);
2740 
2741  tria.set_all_manifold_ids(0);
2742  tria.set_manifold(0, SphericalManifold<2>(center));
2743  }
2744 
2745 
2746  // Implementation for 2D only
2747  template <>
2748  void cylinder(Triangulation<2> &tria,
2749  const double radius,
2750  const double half_length)
2751  {
2752  Point<2> p1(-half_length, -radius);
2753  Point<2> p2(half_length, radius);
2754 
2755  hyper_rectangle(tria, p1, p2, true);
2756 
2757  Triangulation<2>::face_iterator f = tria.begin_face();
2758  Triangulation<2>::face_iterator end = tria.end_face();
2759  while (f != end)
2760  {
2761  switch (f->boundary_id())
2762  {
2763  case 0:
2764  f->set_boundary_id(1);
2765  break;
2766  case 1:
2767  f->set_boundary_id(2);
2768  break;
2769  default:
2770  f->set_boundary_id(0);
2771  break;
2772  }
2773  ++f;
2774  }
2775  }
2776 
2777 
2778 
2779  // Implementation for 2D only
2780  template <>
2781  void cylinder_shell(Triangulation<2> &,
2782  const double,
2783  const double,
2784  const double,
2785  const unsigned int,
2786  const unsigned int)
2787  {
2788  Assert(false, ExcNotImplemented());
2789  }
2790 
2791 
2792  template <>
2793  void quarter_hyper_ball(Triangulation<2> &tria,
2794  const Point<2> & p,
2795  const double radius)
2796  {
2797  const unsigned int dim = 2;
2798 
2799  // equilibrate cell sizes at
2800  // transition from the inner part
2801  // to the radial cells
2802  const Point<dim> vertices[7] = {p + Point<dim>(0, 0) * radius,
2803  p + Point<dim>(+1, 0) * radius,
2804  p + Point<dim>(+1, 0) * (radius / 2),
2805  p + Point<dim>(0, +1) * (radius / 2),
2806  p + Point<dim>(+1, +1) *
2807  (radius / (2 * sqrt(2.0))),
2808  p + Point<dim>(0, +1) * radius,
2809  p + Point<dim>(+1, +1) *
2810  (radius / std::sqrt(2.0))};
2811 
2812  const int cell_vertices[3][4] = {{0, 2, 3, 4}, {1, 6, 2, 4}, {5, 3, 6, 4}};
2813 
2814  std::vector<CellData<dim>> cells(3, CellData<dim>());
2815 
2816  for (unsigned int i = 0; i < 3; ++i)
2817  {
2818  for (unsigned int j = 0; j < 4; ++j)
2819  cells[i].vertices[j] = cell_vertices[i][j];
2820  cells[i].material_id = 0;
2821  }
2822 
2823  tria.create_triangulation(std::vector<Point<dim>>(std::begin(vertices),
2824  std::end(vertices)),
2825  cells,
2826  SubCellData()); // no boundary information
2827 
2828  Triangulation<dim>::cell_iterator cell = tria.begin();
2829  Triangulation<dim>::cell_iterator end = tria.end();
2830 
2831  tria.set_all_manifold_ids_on_boundary(0);
2832 
2833  while (cell != end)
2834  {
2835  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
2836  {
2837  if (cell->face(i)->boundary_id() ==
2838  numbers::internal_face_boundary_id)
2839  continue;
2840 
2841  // If one the components is the same as the respective
2842  // component of the center, then this is part of the plane
2843  if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius ||
2844  cell->face(i)->center()(1) < p(1) + 1.e-5 * radius)
2845  {
2846  cell->face(i)->set_boundary_id(1);
2847  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
2848  }
2849  }
2850  ++cell;
2851  }
2852  tria.set_manifold(0, SphericalManifold<2>(p));
2853  }
2854 
2855 
2856  template <>
2857  void half_hyper_ball(Triangulation<2> &tria,
2858  const Point<2> & p,
2859  const double radius)
2860  {
2861  // equilibrate cell sizes at
2862  // transition from the inner part
2863  // to the radial cells
2864  const double a = 1. / (1 + std::sqrt(2.0));
2865  const Point<2> vertices[8] = {
2866  p + Point<2>(0, -1) * radius,
2867  p + Point<2>(+1, -1) * (radius / std::sqrt(2.0)),
2868  p + Point<2>(0, -1) * (radius / std::sqrt(2.0) * a),
2869  p + Point<2>(+1, -1) * (radius / std::sqrt(2.0) * a),
2870  p + Point<2>(0, +1) * (radius / std::sqrt(2.0) * a),
2871  p + Point<2>(+1, +1) * (radius / std::sqrt(2.0) * a),
2872  p + Point<2>(0, +1) * radius,
2873  p + Point<2>(+1, +1) * (radius / std::sqrt(2.0))};
2874 
2875  const int cell_vertices[5][4] = {{0, 1, 2, 3},
2876  {2, 3, 4, 5},
2877  {1, 7, 3, 5},
2878  {6, 4, 7, 5}};
2879 
2880  std::vector<CellData<2>> cells(4, CellData<2>());
2881 
2882  for (unsigned int i = 0; i < 4; ++i)
2883  {
2884  for (unsigned int j = 0; j < 4; ++j)
2885  cells[i].vertices[j] = cell_vertices[i][j];
2886  cells[i].material_id = 0;
2887  }
2888 
2889  tria.create_triangulation(std::vector<Point<2>>(std::begin(vertices),
2890  std::end(vertices)),
2891  cells,
2892  SubCellData()); // no boundary information
2893 
2894  Triangulation<2>::cell_iterator cell = tria.begin();
2895  Triangulation<2>::cell_iterator end = tria.end();
2896 
2897  tria.set_all_manifold_ids_on_boundary(0);
2898 
2899  while (cell != end)
2900  {
2901  for (unsigned int i = 0; i < GeometryInfo<2>::faces_per_cell; ++i)
2902  {
2903  if (cell->face(i)->boundary_id() ==
2904  numbers::internal_face_boundary_id)
2905  continue;
2906 
2907  // If x is zero, then this is part of the plane
2908  if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius)
2909  {
2910  cell->face(i)->set_boundary_id(1);
2911  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
2912  }
2913  }
2914  ++cell;
2915  }
2916  tria.set_manifold(0, SphericalManifold<2>(p));
2917  }
2918 
2919 
2920 
2921  // Implementation for 2D only
2922  template <>
2923  void half_hyper_shell(Triangulation<2> & tria,
2924  const Point<2> & center,
2925  const double inner_radius,
2926  const double outer_radius,
2927  const unsigned int n_cells,
2928  const bool colorize)
2929  {
2930  Assert((inner_radius > 0) && (inner_radius < outer_radius),
2931  ExcInvalidRadii());
2932 
2933  const double pi = numbers::PI;
2934  // determine the number of cells
2935  // for the grid. if not provided by
2936  // the user determine it such that
2937  // the length of each cell on the
2938  // median (in the middle between
2939  // the two circles) is equal to its
2940  // radial extent (which is the
2941  // difference between the two
2942  // radii)
2943  const unsigned int N =
2944  (n_cells == 0 ? static_cast<unsigned int>(
2945  std::ceil((pi * (outer_radius + inner_radius) / 2) /
2946  (outer_radius - inner_radius))) :
2947  n_cells);
2948 
2949  // set up N+1 vertices on the
2950  // outer and N+1 vertices on
2951  // the inner circle. the
2952  // first N+1 ones are on the
2953  // outer one, and all are
2954  // numbered counter-clockwise
2955  std::vector<Point<2>> vertices(2 * (N + 1));
2956  for (unsigned int i = 0; i <= N; ++i)
2957  {
2958  // enforce that the x-coordinates
2959  // of the first and last point of
2960  // each half-circle are exactly
2961  // zero (contrary to what we may
2962  // compute using the imprecise
2963  // value of pi)
2964  vertices[i] =
2965  Point<2>(((i == 0) || (i == N) ? 0 : std::cos(pi * i / N - pi / 2)),
2966  std::sin(pi * i / N - pi / 2)) *
2967  outer_radius;
2968  vertices[i + N + 1] = vertices[i] * (inner_radius / outer_radius);
2969 
2970  vertices[i] += center;
2971  vertices[i + N + 1] += center;
2972  }
2973 
2974 
2975  std::vector<CellData<2>> cells(N, CellData<2>());
2976 
2977  for (unsigned int i = 0; i < N; ++i)
2978  {
2979  cells[i].vertices[0] = i;
2980  cells[i].vertices[1] = (i + 1) % (N + 1);
2981  cells[i].vertices[2] = N + 1 + i;
2982  cells[i].vertices[3] = N + 1 + ((i + 1) % (N + 1));
2983 
2984  cells[i].material_id = 0;
2985  }
2986 
2987  tria.create_triangulation(vertices, cells, SubCellData());
2988 
2989  if (colorize)
2990  {
2991  Triangulation<2>::cell_iterator cell = tria.begin();
2992  for (; cell != tria.end(); ++cell)
2993  {
2994  cell->face(2)->set_boundary_id(1);
2995  }
2996  tria.begin()->face(0)->set_boundary_id(3);
2997 
2998  tria.last()->face(1)->set_boundary_id(2);
2999  }
3000  tria.set_all_manifold_ids(0);
3001  tria.set_manifold(0, SphericalManifold<2>(center));
3002  }
3003 
3004 
3005  template <>
3006  void quarter_hyper_shell(Triangulation<2> & tria,
3007  const Point<2> & center,
3008  const double inner_radius,
3009  const double outer_radius,
3010  const unsigned int n_cells,
3011  const bool colorize)
3012  {
3013  Assert((inner_radius > 0) && (inner_radius < outer_radius),
3014  ExcInvalidRadii());
3015 
3016  const double pi = numbers::PI;
3017  // determine the number of cells
3018  // for the grid. if not provided by
3019  // the user determine it such that
3020  // the length of each cell on the
3021  // median (in the middle between
3022  // the two circles) is equal to its
3023  // radial extent (which is the
3024  // difference between the two
3025  // radii)
3026  const unsigned int N =
3027  (n_cells == 0 ? static_cast<unsigned int>(
3028  std::ceil((pi * (outer_radius + inner_radius) / 4) /
3029  (outer_radius - inner_radius))) :
3030  n_cells);
3031 
3032  // set up N+1 vertices on the
3033  // outer and N+1 vertices on
3034  // the inner circle. the
3035  // first N+1 ones are on the
3036  // outer one, and all are
3037  // numbered counter-clockwise
3038  std::vector<Point<2>> vertices(2 * (N + 1));
3039  for (unsigned int i = 0; i <= N; ++i)
3040  {
3041  // enforce that the x-coordinates
3042  // of the last point is exactly
3043  // zero (contrary to what we may
3044  // compute using the imprecise
3045  // value of pi)
3046  vertices[i] = Point<2>(((i == N) ? 0 : std::cos(pi * i / N / 2)),
3047  std::sin(pi * i / N / 2)) *
3048  outer_radius;
3049  vertices[i + N + 1] = vertices[i] * (inner_radius / outer_radius);
3050 
3051  vertices[i] += center;
3052  vertices[i + N + 1] += center;
3053  }
3054 
3055 
3056  std::vector<CellData<2>> cells(N, CellData<2>());
3057 
3058  for (unsigned int i = 0; i < N; ++i)
3059  {
3060  cells[i].vertices[0] = i;
3061  cells[i].vertices[1] = (i + 1) % (N + 1);
3062  cells[i].vertices[2] = N + 1 + i;
3063  cells[i].vertices[3] = N + 1 + ((i + 1) % (N + 1));
3064 
3065  cells[i].material_id = 0;
3066  }
3067 
3068  tria.create_triangulation(vertices, cells, SubCellData());
3069 
3070  if (colorize)
3071  {
3072  Triangulation<2>::cell_iterator cell = tria.begin();
3073  for (; cell != tria.end(); ++cell)
3074  {
3075  cell->face(2)->set_boundary_id(1);
3076  }
3077  tria.begin()->face(0)->set_boundary_id(3);
3078 
3079  tria.last()->face(1)->set_boundary_id(2);
3080  }
3081 
3082  tria.set_all_manifold_ids(0);
3083  tria.set_manifold(0, SphericalManifold<2>(center));
3084  }
3085 
3086 
3087 
3088  // Implementation for 3D only
3089  template <>
3090  void hyper_cube_slit(Triangulation<3> &tria,
3091  const double left,
3092  const double right,
3093  const bool colorize)
3094  {
3095  const double rl2 = (right + left) / 2;
3096  const double len = (right - left) / 2.;
3097 
3098  const Point<3> vertices[20] = {
3099  Point<3>(left, left, -len / 2.), Point<3>(rl2, left, -len / 2.),
3100  Point<3>(rl2, rl2, -len / 2.), Point<3>(left, rl2, -len / 2.),
3101  Point<3>(right, left, -len / 2.), Point<3>(right, rl2, -len / 2.),
3102  Point<3>(rl2, right, -len / 2.), Point<3>(left, right, -len / 2.),
3103  Point<3>(right, right, -len / 2.), Point<3>(rl2, left, -len / 2.),
3104  Point<3>(left, left, len / 2.), Point<3>(rl2, left, len / 2.),
3105  Point<3>(rl2, rl2, len / 2.), Point<3>(left, rl2, len / 2.),
3106  Point<3>(right, left, len / 2.), Point<3>(right, rl2, len / 2.),
3107  Point<3>(rl2, right, len / 2.), Point<3>(left, right, len / 2.),
3108  Point<3>(right, right, len / 2.), Point<3>(rl2, left, len / 2.)};
3109  const int cell_vertices[4][8] = {{0, 1, 3, 2, 10, 11, 13, 12},
3110  {9, 4, 2, 5, 19, 14, 12, 15},
3111  {3, 2, 7, 6, 13, 12, 17, 16},
3112  {2, 5, 6, 8, 12, 15, 16, 18}};
3113  std::vector<CellData<3>> cells(4, CellData<3>());
3114  for (unsigned int i = 0; i < 4; ++i)
3115  {
3116  for (unsigned int j = 0; j < 8; ++j)
3117  cells[i].vertices[j] = cell_vertices[i][j];
3118  cells[i].material_id = 0;
3119  }
3120  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
3121  std::end(vertices)),
3122  cells,
3123  SubCellData()); // no boundary information
3124 
3125  if (colorize)
3126  {
3127  Triangulation<3>::cell_iterator cell = tria.begin();
3128  cell->face(1)->set_boundary_id(1);
3129  ++cell;
3130  cell->face(0)->set_boundary_id(2);
3131  }
3132  }
3133 
3134 
3135 
3136  // Implementation for 3D only
3137  template <>
3138  void enclosed_hyper_cube(Triangulation<3> &tria,
3139  const double left,
3140  const double right,
3141  const double thickness,
3142  const bool colorize)
3143  {
3144  Assert(left < right,
3145  ExcMessage("Invalid left-to-right bounds of enclosed hypercube"));
3146 
3147  std::vector<Point<3>> vertices(64);
3148  double coords[4];
3149  coords[0] = left - thickness;
3150  coords[1] = left;
3151  coords[2] = right;
3152  coords[3] = right + thickness;
3153 
3154  unsigned int k = 0;
3155  for (unsigned int z = 0; z < 4; ++z)
3156  for (unsigned int y = 0; y < 4; ++y)
3157  for (unsigned int x = 0; x < 4; ++x)
3158  vertices[k++] = Point<3>(coords[x], coords[y], coords[z]);
3159 
3160  const types::material_id materials[27] = {21, 20, 22, 17, 16, 18, 25,
3161  24, 26, 5, 4, 6, 1, 0,
3162  2, 9, 8, 10, 37, 36, 38,
3163  33, 32, 34, 41, 40, 42};
3164 
3165  std::vector<CellData<3>> cells(27);
3166  k = 0;
3167  for (unsigned int z = 0; z < 3; ++z)
3168  for (unsigned int y = 0; y < 3; ++y)
3169  for (unsigned int x = 0; x < 3; ++x)
3170  {
3171  cells[k].vertices[0] = x + 4 * y + 16 * z;
3172  cells[k].vertices[1] = x + 4 * y + 16 * z + 1;
3173  cells[k].vertices[2] = x + 4 * y + 16 * z + 4;
3174  cells[k].vertices[3] = x + 4 * y + 16 * z + 5;
3175  cells[k].vertices[4] = x + 4 * y + 16 * z + 16;
3176  cells[k].vertices[5] = x + 4 * y + 16 * z + 17;
3177  cells[k].vertices[6] = x + 4 * y + 16 * z + 20;
3178  cells[k].vertices[7] = x + 4 * y + 16 * z + 21;
3179  if (colorize)
3180  cells[k].material_id = materials[k];
3181  ++k;
3182  }
3183  tria.create_triangulation(vertices,
3184  cells,
3185  SubCellData()); // no boundary information
3186  }
3187 
3188 
3189 
3190  template <>
3191  void truncated_cone(Triangulation<3> &triangulation,
3192  const double radius_0,
3193  const double radius_1,
3194  const double half_length)
3195  {
3196  // Determine number of cells and vertices
3197  const unsigned int n_cells = static_cast<unsigned int>(
3198  std::ceil(half_length / std::max(radius_0, radius_1)));
3199  const unsigned int n_vertices = 4 * (n_cells + 1);
3200  std::vector<Point<3>> vertices_tmp(n_vertices);
3201 
3202  vertices_tmp[0] = Point<3>(-half_length, 0, -radius_0);
3203  vertices_tmp[1] = Point<3>(-half_length, radius_0, 0);
3204  vertices_tmp[2] = Point<3>(-half_length, -radius_0, 0);
3205  vertices_tmp[3] = Point<3>(-half_length, 0, radius_0);
3206 
3207  const double dx = 2 * half_length / n_cells;
3208 
3209  for (unsigned int i = 0; i < n_cells; ++i)
3210  {
3211  vertices_tmp[4 * (i + 1)] =
3212  vertices_tmp[4 * i] +
3213  Point<3>(dx, 0, 0.5 * (radius_0 - radius_1) * dx / half_length);
3214  vertices_tmp[4 * i + 5] =
3215  vertices_tmp[4 * i + 1] +
3216  Point<3>(dx, 0.5 * (radius_1 - radius_0) * dx / half_length, 0);
3217  vertices_tmp[4 * i + 6] =
3218  vertices_tmp[4 * i + 2] +
3219  Point<3>(dx, 0.5 * (radius_0 - radius_1) * dx / half_length, 0);
3220  vertices_tmp[4 * i + 7] =
3221  vertices_tmp[4 * i + 3] +
3222  Point<3>(dx, 0, 0.5 * (radius_1 - radius_0) * dx / half_length);
3223  }
3224 
3225  const std::vector<Point<3>> vertices(vertices_tmp.begin(),
3226  vertices_tmp.end());
3227  Table<2, unsigned int> cell_vertices(n_cells,
3228  GeometryInfo<3>::vertices_per_cell);
3229 
3230  for (unsigned int i = 0; i < n_cells; ++i)
3231  for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j)
3232  cell_vertices[i][j] = 4 * i + j;
3233 
3234  std::vector<CellData<3>> cells(n_cells, CellData<3>());
3235 
3236  for (unsigned int i = 0; i < n_cells; ++i)
3237  {
3238  for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j)
3239  cells[i].vertices[j] = cell_vertices[i][j];
3240 
3241  cells[i].material_id = 0;
3242  }
3243 
3244  triangulation.create_triangulation(vertices, cells, SubCellData());
3245  triangulation.set_all_manifold_ids_on_boundary(0);
3246 
3247  for (Triangulation<3>::cell_iterator cell = triangulation.begin();
3248  cell != triangulation.end();
3249  ++cell)
3250  {
3251  if (cell->vertex(0)(0) == -half_length)
3252  {
3253  cell->face(4)->set_boundary_id(1);
3254  cell->face(4)->set_manifold_id(numbers::flat_manifold_id);
3255 
3256  for (unsigned int i = 0; i < 4; ++i)
3257  cell->line(i)->set_boundary_id(0);
3258  }
3259 
3260  if (cell->vertex(4)(0) == half_length)
3261  {
3262  cell->face(5)->set_boundary_id(2);
3263  cell->face(5)->set_manifold_id(numbers::flat_manifold_id);
3264 
3265  for (unsigned int i = 4; i < 8; ++i)
3266  cell->line(i)->set_boundary_id(0);
3267  }
3268 
3269  for (unsigned int i = 0; i < 4; ++i)
3270  cell->face(i)->set_boundary_id(0);
3271  }
3272 
3273  triangulation.set_manifold(0, CylindricalManifold<3>());
3274  }
3275 
3276 
3277  // Implementation for 3D only
3278  template <>
3279  void hyper_L(Triangulation<3> &tria,
3280  const double a,
3281  const double b,
3282  const bool colorize)
3283  {
3284  // we slice out the top back right
3285  // part of the cube
3286  const Point<3> vertices[26] = {
3287  // front face of the big cube
3288  Point<3>(a, a, a),
3289  Point<3>((a + b) / 2, a, a),
3290  Point<3>(b, a, a),
3291  Point<3>(a, a, (a + b) / 2),
3292  Point<3>((a + b) / 2, a, (a + b) / 2),
3293  Point<3>(b, a, (a + b) / 2),
3294  Point<3>(a, a, b),
3295  Point<3>((a + b) / 2, a, b),
3296  Point<3>(b, a, b),
3297  // middle face of the big cube
3298  Point<3>(a, (a + b) / 2, a),
3299  Point<3>((a + b) / 2, (a + b) / 2, a),
3300  Point<3>(b, (a + b) / 2, a),
3301  Point<3>(a, (a + b) / 2, (a + b) / 2),
3302  Point<3>((a + b) / 2, (a + b) / 2, (a + b) / 2),
3303  Point<3>(b, (a + b) / 2, (a + b) / 2),
3304  Point<3>(a, (a + b) / 2, b),
3305  Point<3>((a + b) / 2, (a + b) / 2, b),
3306  Point<3>(b, (a + b) / 2, b),
3307  // back face of the big cube
3308  // last (top right) point is missing
3309  Point<3>(a, b, a),
3310  Point<3>((a + b) / 2, b, a),
3311  Point<3>(b, b, a),
3312  Point<3>(a, b, (a + b) / 2),
3313  Point<3>((a + b) / 2, b, (a + b) / 2),
3314  Point<3>(b, b, (a + b) / 2),
3315  Point<3>(a, b, b),
3316  Point<3>((a + b) / 2, b, b)};
3317  const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
3318  {1, 2, 10, 11, 4, 5, 13, 14},
3319  {3, 4, 12, 13, 6, 7, 15, 16},
3320  {4, 5, 13, 14, 7, 8, 16, 17},
3321  {9, 10, 18, 19, 12, 13, 21, 22},
3322  {10, 11, 19, 20, 13, 14, 22, 23},
3323  {12, 13, 21, 22, 15, 16, 24, 25}};
3324 
3325  std::vector<CellData<3>> cells(7, CellData<3>());
3326 
3327  for (unsigned int i = 0; i < 7; ++i)
3328  {
3329  for (unsigned int j = 0; j < 8; ++j)
3330  cells[i].vertices[j] = cell_vertices[i][j];
3331  cells[i].material_id = 0;
3332  }
3333 
3334  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
3335  std::end(vertices)),
3336  cells,
3337  SubCellData()); // no boundary information
3338 
3339  if (colorize)
3340  {
3341  Assert(false, ExcNotImplemented());
3342  }
3343  }
3344 
3345 
3346 
3347  // Implementation for 3D only
3348  template <>
3349  void hyper_ball(Triangulation<3> &tria,
3350  const Point<3> & p,
3351  const double radius,
3352  const bool internal_manifold)
3353  {
3354  const double a =
3355  1. / (1 + std::sqrt(3.0)); // equilibrate cell sizes at transition
3356  // from the inner part to the radial
3357  // cells
3358  const unsigned int n_vertices = 16;
3359  const Point<3> vertices[n_vertices] = {
3360  // first the vertices of the inner
3361  // cell
3362  p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0) * a),
3363  p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0) * a),
3364  p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0) * a),
3365  p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0) * a),
3366  p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0) * a),
3367  p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0) * a),
3368  p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0) * a),
3369  p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0) * a),
3370  // now the eight vertices at
3371  // the outer sphere
3372  p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0)),
3373  p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0)),
3374  p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0)),
3375  p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0)),
3376  p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0)),
3377  p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0)),
3378  p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0)),
3379  p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0)),
3380  };
3381 
3382  // one needs to draw the seven cubes to
3383  // understand what's going on here
3384  const unsigned int n_cells = 7;
3385  const int cell_vertices[n_cells][8] = {
3386  {0, 1, 4, 5, 3, 2, 7, 6}, // center
3387  {8, 9, 12, 13, 0, 1, 4, 5}, // bottom
3388  {9, 13, 1, 5, 10, 14, 2, 6}, // right
3389  {11, 10, 3, 2, 15, 14, 7, 6}, // top
3390  {8, 0, 12, 4, 11, 3, 15, 7}, // left
3391  {8, 9, 0, 1, 11, 10, 3, 2}, // front
3392  {12, 4, 13, 5, 15, 7, 14, 6}}; // back
3393 
3394  std::vector<CellData<3>> cells(n_cells, CellData<3>());
3395 
3396  for (unsigned int i = 0; i < n_cells; ++i)
3397  {
3398  for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j)
3399  cells[i].vertices[j] = cell_vertices[i][j];
3400  cells[i].material_id = 0;
3401  cells[i].manifold_id = i == 0 ? numbers::flat_manifold_id : 1;
3402  }
3403 
3404  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
3405  std::end(vertices)),
3406  cells,
3407  SubCellData()); // no boundary information
3408  tria.set_all_manifold_ids_on_boundary(0);
3409  tria.set_manifold(0, SphericalManifold<3>(p));
3410  if (internal_manifold)
3411  tria.set_manifold(1, SphericalManifold<3>(p));
3412  }
3413 
3414 
3415 
3416  template <int spacedim>
3417  void hyper_sphere(Triangulation<spacedim - 1, spacedim> &tria,
3418  const Point<spacedim> & p,
3419  const double radius)
3420  {
3421  Triangulation<spacedim> volume_mesh;
3422  GridGenerator::hyper_ball(volume_mesh, p, radius);
3423  std::set<types::boundary_id> boundary_ids;
3424  boundary_ids.insert(0);
3425  GridGenerator::extract_boundary_mesh(volume_mesh, tria, boundary_ids);
3426  tria.set_all_manifold_ids(0);
3427  tria.set_manifold(0, SphericalManifold<spacedim - 1, spacedim>(p));
3428  }
3429 
3430 
3431 
3432  // Implementation for 3D only
3433  template <>
3434  void cylinder(Triangulation<3> &tria,
3435  const double radius,
3436  const double half_length)
3437  {
3438  // Copy the base from hyper_ball<3>
3439  // and transform it to yz
3440  const double d = radius / std::sqrt(2.0);
3441  const double a = d / (1 + std::sqrt(2.0));
3442  Point<3> vertices[24] = {
3443  Point<3>(-d, -half_length, -d),
3444  Point<3>(d, -half_length, -d),
3445  Point<3>(-a, -half_length, -a),
3446  Point<3>(a, -half_length, -a),
3447  Point<3>(-a, -half_length, a),
3448  Point<3>(a, -half_length, a),
3449  Point<3>(-d, -half_length, d),
3450  Point<3>(d, -half_length, d),
3451  Point<3>(-d, 0, -d),
3452  Point<3>(d, 0, -d),
3453  Point<3>(-a, 0, -a),
3454  Point<3>(a, 0, -a),
3455  Point<3>(-a, 0, a),
3456  Point<3>(a, 0, a),
3457  Point<3>(-d, 0, d),
3458  Point<3>(d, 0, d),
3459  Point<3>(-d, half_length, -d),
3460  Point<3>(d, half_length, -d),
3461  Point<3>(-a, half_length, -a),
3462  Point<3>(a, half_length, -a),
3463  Point<3>(-a, half_length, a),
3464  Point<3>(a, half_length, a),
3465  Point<3>(-d, half_length, d),
3466  Point<3>(d, half_length, d),
3467  };
3468  // Turn cylinder such that y->x
3469  for (unsigned int i = 0; i < 24; ++i)
3470  {
3471  const double h = vertices[i](1);
3472  vertices[i](1) = -vertices[i](0);
3473  vertices[i](0) = h;
3474  }
3475 
3476  int cell_vertices[10][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
3477  {0, 2, 8, 10, 6, 4, 14, 12},
3478  {2, 3, 10, 11, 4, 5, 12, 13},
3479  {1, 7, 9, 15, 3, 5, 11, 13},
3480  {6, 4, 14, 12, 7, 5, 15, 13}};
3481  for (unsigned int i = 0; i < 5; ++i)
3482  for (unsigned int j = 0; j < 8; ++j)
3483  cell_vertices[i + 5][j] = cell_vertices[i][j] + 8;
3484 
3485  std::vector<CellData<3>> cells(10, CellData<3>());
3486 
3487  for (unsigned int i = 0; i < 10; ++i)
3488  {
3489  for (unsigned int j = 0; j < 8; ++j)
3490  cells[i].vertices[j] = cell_vertices[i][j];
3491  cells[i].material_id = 0;
3492  }
3493 
3494  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
3495  std::end(vertices)),
3496  cells,
3497  SubCellData()); // no boundary information
3498 
3499  // set boundary indicators for the
3500  // faces at the ends to 1 and 2,
3501  // respectively. note that we also
3502  // have to deal with those lines
3503  // that are purely in the interior
3504  // of the ends. we determine whether
3505  // an edge is purely in the
3506  // interior if one of its vertices
3507  // is at coordinates '+-a' as set
3508  // above
3509  Triangulation<3>::cell_iterator cell = tria.begin();
3510  Triangulation<3>::cell_iterator end = tria.end();
3511 
3512  tria.set_all_manifold_ids_on_boundary(0);
3513 
3514  for (; cell != end; ++cell)
3515  for (unsigned int i = 0; i < GeometryInfo<3>::faces_per_cell; ++i)
3516  if (cell->at_boundary(i))
3517  {
3518  if (cell->face(i)->center()(0) > half_length - 1.e-5)
3519  {
3520  cell->face(i)->set_boundary_id(2);
3521  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
3522 
3523  for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
3524  ++e)
3525  if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
3526  (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
3527  (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
3528  (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
3529  {
3530  cell->face(i)->line(e)->set_boundary_id(2);
3531  cell->face(i)->line(e)->set_manifold_id(
3532  numbers::flat_manifold_id);
3533  }
3534  }
3535  else if (cell->face(i)->center()(0) < -half_length + 1.e-5)
3536  {
3537  cell->face(i)->set_boundary_id(1);
3538  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
3539 
3540  for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
3541  ++e)
3542  if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
3543  (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
3544  (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
3545  (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
3546  {
3547  cell->face(i)->line(e)->set_boundary_id(1);
3548  cell->face(i)->line(e)->set_manifold_id(
3549  numbers::flat_manifold_id);
3550  }
3551  }
3552  }
3553  tria.set_manifold(0, CylindricalManifold<3>());
3554  }
3555 
3556 
3557  template <>
3558  void quarter_hyper_ball(Triangulation<3> &tria,
3559  const Point<3> & center,
3560  const double radius)
3561  {
3562  const unsigned int dim = 3;
3563 
3564  // equilibrate cell sizes at
3565  // transition from the inner part
3566  // to the radial cells
3567  const Point<dim> vertices[15] = {
3568  center + Point<dim>(0, 0, 0) * radius,
3569  center + Point<dim>(+1, 0, 0) * radius,
3570  center + Point<dim>(+1, 0, 0) * (radius / 2.),
3571  center + Point<dim>(0, +1, 0) * (radius / 2.),
3572  center + Point<dim>(+1, +1, 0) * (radius / (2 * sqrt(2.0))),
3573  center + Point<dim>(0, +1, 0) * radius,
3574  center + Point<dim>(+1, +1, 0) * (radius / std::sqrt(2.0)),
3575  center + Point<dim>(0, 0, 1) * radius / 2.,
3576  center + Point<dim>(+1, 0, 1) * radius / std::sqrt(2.0),
3577  center + Point<dim>(+1, 0, 1) * (radius / (2 * std::sqrt(2.0))),
3578  center + Point<dim>(0, +1, 1) * (radius / (2 * std::sqrt(2.0))),
3579  center + Point<dim>(+1, +1, 1) * (radius / (2 * std::sqrt(3.0))),
3580  center + Point<dim>(0, +1, 1) * radius / std::sqrt(2.0),
3581  center + Point<dim>(+1, +1, 1) * (radius / (std::sqrt(3.0))),
3582  center + Point<dim>(0, 0, 1) * radius};
3583  const int cell_vertices[4][8] = {{0, 2, 3, 4, 7, 9, 10, 11},
3584  {1, 6, 2, 4, 8, 13, 9, 11},
3585  {5, 3, 6, 4, 12, 10, 13, 11},
3586  {7, 9, 10, 11, 14, 8, 12, 13}};
3587 
3588  std::vector<CellData<dim>> cells(4, CellData<dim>());
3589 
3590  for (unsigned int i = 0; i < 4; ++i)
3591  {
3592  for (unsigned int j = 0; j < 8; ++j)
3593  cells[i].vertices[j] = cell_vertices[i][j];
3594  cells[i].material_id = 0;
3595  }
3596 
3597  tria.create_triangulation(std::vector<Point<dim>>(std::begin(vertices),
3598  std::end(vertices)),
3599  cells,
3600  SubCellData()); // no boundary information
3601 
3602  Triangulation<dim>::cell_iterator cell = tria.begin();
3603  Triangulation<dim>::cell_iterator end = tria.end();
3604 
3605  tria.set_all_manifold_ids_on_boundary(0);
3606  while (cell != end)
3607  {
3608  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
3609  {
3610  if (cell->face(i)->boundary_id() ==
3611  numbers::internal_face_boundary_id)
3612  continue;
3613 
3614  // If x,y or z is zero, then this is part of the plane
3615  if (cell->face(i)->center()(0) < center(0) + 1.e-5 * radius ||
3616  cell->face(i)->center()(1) < center(1) + 1.e-5 * radius ||
3617  cell->face(i)->center()(2) < center(2) + 1.e-5 * radius)
3618  {
3619  cell->face(i)->set_boundary_id(1);
3620  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
3621  // also set the boundary indicators of the bounding lines,
3622  // unless both vertices are on the perimeter
3623  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
3624  ++j)
3625  {
3626  const Point<3> line_vertices[2] = {
3627  cell->face(i)->line(j)->vertex(0),
3628  cell->face(i)->line(j)->vertex(1)};
3629  if ((std::fabs(line_vertices[0].distance(center) - radius) >
3630  1e-5 * radius) ||
3631  (std::fabs(line_vertices[1].distance(center) - radius) >
3632  1e-5 * radius))
3633  {
3634  cell->face(i)->line(j)->set_boundary_id(1);
3635  cell->face(i)->line(j)->set_manifold_id(
3636  numbers::flat_manifold_id);
3637  }
3638  }
3639  }
3640  }
3641  ++cell;
3642  }
3643  tria.set_manifold(0, SphericalManifold<3>(center));
3644  }
3645 
3646 
3647  // Implementation for 3D only
3648  template <>
3649  void half_hyper_ball(Triangulation<3> &tria,
3650  const Point<3> & center,
3651  const double radius)
3652  {
3653  // These are for the two lower squares
3654  const double d = radius / std::sqrt(2.0);
3655  const double a = d / (1 + std::sqrt(2.0));
3656  // These are for the two upper square
3657  const double b = a / 2.0;
3658  const double c = d / 2.0;
3659  // And so are these
3660  const double hb = radius * std::sqrt(3.0) / 4.0;
3661  const double hc = radius * std::sqrt(3.0) / 2.0;
3662 
3663  Point<3> vertices[16] = {
3664  center + Point<3>(0, d, -d),
3665  center + Point<3>(0, -d, -d),
3666  center + Point<3>(0, a, -a),
3667  center + Point<3>(0, -a, -a),
3668  center + Point<3>(0, a, a),
3669  center + Point<3>(0, -a, a),
3670  center + Point<3>(0, d, d),
3671  center + Point<3>(0, -d, d),
3672 
3673  center + Point<3>(hc, c, -c),
3674  center + Point<3>(hc, -c, -c),
3675  center + Point<3>(hb, b, -b),
3676  center + Point<3>(hb, -b, -b),
3677  center + Point<3>(hb, b, b),
3678  center + Point<3>(hb, -b, b),
3679  center + Point<3>(hc, c, c),
3680  center + Point<3>(hc, -c, c),
3681  };
3682 
3683  int cell_vertices[6][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
3684  {0, 2, 8, 10, 6, 4, 14, 12},
3685  {2, 3, 10, 11, 4, 5, 12, 13},
3686  {1, 7, 9, 15, 3, 5, 11, 13},
3687  {6, 4, 14, 12, 7, 5, 15, 13},
3688  {8, 10, 9, 11, 14, 12, 15, 13}};
3689 
3690  std::vector<CellData<3>> cells(6, CellData<3>());
3691 
3692  for (unsigned int i = 0; i < 6; ++i)
3693  {
3694  for (unsigned int j = 0; j < 8; ++j)
3695  cells[i].vertices[j] = cell_vertices[i][j];
3696  cells[i].material_id = 0;
3697  }
3698 
3699  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
3700  std::end(vertices)),
3701  cells,
3702  SubCellData()); // no boundary information
3703 
3704  Triangulation<3>::cell_iterator cell = tria.begin();
3705  Triangulation<3>::cell_iterator end = tria.end();
3706 
3707  tria.set_all_manifold_ids_on_boundary(0);
3708 
3709  // go over all faces. for the ones on the flat face, set boundary
3710  // indicator for face and edges to one; the rest will remain at
3711  // zero but we have to pay attention to those edges that are
3712  // at the perimeter of the flat face since they should not be
3713  // set to one
3714  while (cell != end)
3715  {
3716  for (unsigned int i = 0; i < GeometryInfo<3>::faces_per_cell; ++i)
3717  {
3718  if (!cell->at_boundary(i))
3719  continue;
3720 
3721  // If the center is on the plane x=0, this is a planar element. set
3722  // its boundary indicator. also set the boundary indicators of the
3723  // bounding faces unless both vertices are on the perimeter
3724  if (cell->face(i)->center()(0) < center(0) + 1.e-5 * radius)
3725  {
3726  cell->face(i)->set_boundary_id(1);
3727  cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
3728  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
3729  ++j)
3730  {
3731  const Point<3> line_vertices[2] = {
3732  cell->face(i)->line(j)->vertex(0),
3733  cell->face(i)->line(j)->vertex(1)};
3734  if ((std::fabs(line_vertices[0].distance(center) - radius) >
3735  1e-5 * radius) ||
3736  (std::fabs(line_vertices[1].distance(center) - radius) >
3737  1e-5 * radius))
3738  {
3739  cell->face(i)->line(j)->set_boundary_id(1);
3740  cell->face(i)->line(j)->set_manifold_id(
3741  numbers::flat_manifold_id);
3742  }
3743  }
3744  }
3745  }
3746  ++cell;
3747  }
3748  tria.set_manifold(0, SphericalManifold<3>(center));
3749  }
3750 
3751 
3752  template <>
3753  void hyper_shell(Triangulation<3> & tria,
3754  const Point<3> & p,
3755  const double inner_radius,
3756  const double outer_radius,
3757  const unsigned int n_cells,
3758  const bool colorize)
3759  {
3760  Assert((inner_radius > 0) && (inner_radius < outer_radius),
3761  ExcInvalidRadii());
3762 
3763  const unsigned int n = (n_cells == 0) ? 6 : n_cells;
3764 
3765  const double irad = inner_radius / std::sqrt(3.0);
3766  const double orad = outer_radius / std::sqrt(3.0);
3767  std::vector<Point<3>> vertices;
3768  std::vector<CellData<3>> cells;
3769 
3770  // Start with the shell bounded by
3771  // two nested cubes
3772  if (n == 6)
3773  {
3774  for (unsigned int i = 0; i < 8; ++i)
3775  vertices.push_back(p + hexahedron[i] * irad);
3776  for (unsigned int i = 0; i < 8; ++i)
3777  vertices.push_back(p + hexahedron[i] * orad);
3778 
3779  const unsigned int n_cells = 6;
3780  const int cell_vertices[n_cells][8] = {
3781  {8, 9, 10, 11, 0, 1, 2, 3}, // bottom
3782  {9, 11, 1, 3, 13, 15, 5, 7}, // right
3783  {12, 13, 4, 5, 14, 15, 6, 7}, // top
3784  {8, 0, 10, 2, 12, 4, 14, 6}, // left
3785  {8, 9, 0, 1, 12, 13, 4, 5}, // front
3786  {10, 2, 11, 3, 14, 6, 15, 7}}; // back
3787 
3788  cells.resize(n_cells, CellData<3>());
3789 
3790  for (unsigned int i = 0; i < n_cells; ++i)
3791  {
3792  for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell;
3793  ++j)
3794  cells[i].vertices[j] = cell_vertices[i][j];
3795  cells[i].material_id = 0;
3796  }
3797 
3798  tria.create_triangulation(vertices, cells, SubCellData());
3799  }
3800  // A more regular subdivision can
3801  // be obtained by two nested
3802  // rhombic dodecahedra
3803  else if (n == 12)
3804  {
3805  for (unsigned int i = 0; i < 8; ++i)
3806  vertices.push_back(p + hexahedron[i] * irad);
3807  for (unsigned int i = 0; i < 6; ++i)
3808  vertices.push_back(p + octahedron[i] * inner_radius);
3809  for (unsigned int i = 0; i < 8; ++i)
3810  vertices.push_back(p + hexahedron[i] * orad);
3811  for (unsigned int i = 0; i < 6; ++i)
3812  vertices.push_back(p + octahedron[i] * outer_radius);
3813 
3814  const unsigned int n_cells = 12;
3815  const unsigned int rhombi[n_cells][4] = {{10, 4, 0, 8},
3816  {4, 13, 8, 6},
3817  {10, 5, 4, 13},
3818  {1, 9, 10, 5},
3819  {9, 7, 5, 13},
3820  {7, 11, 13, 6},
3821  {9, 3, 7, 11},
3822  {1, 12, 9, 3},
3823  {12, 2, 3, 11},
3824  {2, 8, 11, 6},
3825  {12, 0, 2, 8},
3826  {1, 10, 12, 0}};
3827 
3828  cells.resize(n_cells, CellData<3>());
3829 
3830  for (unsigned int i = 0; i < n_cells; ++i)
3831  {
3832  for (unsigned int j = 0; j < 4; ++j)
3833  {
3834  cells[i].vertices[j] = rhombi[i][j];
3835  cells[i].vertices[j + 4] = rhombi[i][j] + 14;
3836  }
3837  cells[i].material_id = 0;
3838  }
3839 
3840  tria.create_triangulation(vertices, cells, SubCellData());
3841  }
3842  else if (n == 96)
3843  {
3844  // create a triangulation based on the 12-cell version. This function
3845  // was needed before SphericalManifold was written: it manually
3846  // adjusted the interior vertices to lie along concentric
3847  // spheres. Nowadays we can just refine globally:
3848  Triangulation<3> tmp;
3849  hyper_shell(tmp, p, inner_radius, outer_radius, 12);
3850  tmp.refine_global(1);
3851 
3852  // now copy the resulting level 1 cells into the new triangulation,
3853  cells.resize(tmp.n_active_cells(), CellData<3>());
3854  for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active();
3855  cell != tmp.end();
3856  ++cell)
3857  {
3858  const unsigned int cell_index = cell->active_cell_index();
3859  for (unsigned int v = 0; v < GeometryInfo<3>::vertices_per_cell;
3860  ++v)
3861  cells[cell_index].vertices[v] = cell->vertex_index(v);
3862  cells[cell_index].material_id = 0;
3863  }
3864 
3865  tria.create_triangulation(tmp.get_vertices(), cells, SubCellData());
3866  }
3867  else
3868  {
3869  Assert(false, ExcMessage("Invalid number of coarse mesh cells."));
3870  }
3871 
3872  if (colorize)
3873  colorize_hyper_shell(tria, p, inner_radius, outer_radius);
3874  tria.set_all_manifold_ids(0);
3875  tria.set_manifold(0, SphericalManifold<3>(p));
3876  }
3877 
3878 
3879 
3880  // Implementation for 3D only
3881  template <>
3882  void half_hyper_shell(Triangulation<3> & tria,
3883  const Point<3> & center,
3884  const double inner_radius,
3885  const double outer_radius,
3886  const unsigned int n,
3887  const bool colorize)
3888  {
3889  Assert((inner_radius > 0) && (inner_radius < outer_radius),
3890  ExcInvalidRadii());
3891 
3892  if (n <= 5)
3893  {
3894  // These are for the two lower squares
3895  const double d = outer_radius / std::sqrt(2.0);
3896  const double a = inner_radius / std::sqrt(2.0);
3897  // These are for the two upper square
3898  const double b = a / 2.0;
3899  const double c = d / 2.0;
3900  // And so are these
3901  const double hb = inner_radius * std::sqrt(3.0) / 2.0;
3902  const double hc = outer_radius * std::sqrt(3.0) / 2.0;
3903 
3904  Point<3> vertices[16] = {
3905  center + Point<3>(0, d, -d),
3906  center + Point<3>(0, -d, -d),
3907  center + Point<3>(0, a, -a),
3908  center + Point<3>(0, -a, -a),
3909  center + Point<3>(0, a, a),
3910  center + Point<3>(0, -a, a),
3911  center + Point<3>(0, d, d),
3912  center + Point<3>(0, -d, d),
3913 
3914  center + Point<3>(hc, c, -c),
3915  center + Point<3>(hc, -c, -c),
3916  center + Point<3>(hb, b, -b),
3917  center + Point<3>(hb, -b, -b),
3918  center + Point<3>(hb, b, b),
3919  center + Point<3>(hb, -b, b),
3920  center + Point<3>(hc, c, c),
3921  center + Point<3>(hc, -c, c),
3922  };
3923 
3924  int cell_vertices[5][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
3925  {0, 2, 8, 10, 6, 4, 14, 12},
3926  {1, 7, 9, 15, 3, 5, 11, 13},
3927  {6, 4, 14, 12, 7, 5, 15, 13},
3928  {8, 10, 9, 11, 14, 12, 15, 13}};
3929 
3930  std::vector<CellData<3>> cells(5, CellData<3>());
3931 
3932  for (unsigned int i = 0; i < 5; ++i)
3933  {
3934  for (unsigned int j = 0; j < 8; ++j)
3935  cells[i].vertices[j] = cell_vertices[i][j];
3936  cells[i].material_id = 0;
3937  }
3938 
3939  tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
3940  std::end(vertices)),
3941  cells,
3942  SubCellData()); // no boundary information
3943  }
3944  else
3945  {
3946  Assert(false, ExcIndexRange(n, 0, 5));
3947  }
3948  if (colorize)
3949  {
3950  // We want to use a standard boundary description where
3951  // the boundary is not curved. Hence set boundary id 2 to
3952  // to all faces in a first step.
3953  Triangulation<3>::cell_iterator cell = tria.begin();
3954  for (; cell != tria.end(); ++cell)
3955  for (unsigned int i = 0; i < GeometryInfo<3>::faces_per_cell; ++i)
3956  if (cell->at_boundary(i))
3957  cell->face(i)->set_all_boundary_ids(2);
3958 
3959  // Next look for the curved boundaries. If the x value of the
3960  // center of the face is not equal to center(0), we're on a curved
3961  // boundary. Then decide whether the center is nearer to the inner
3962  // or outer boundary to set the correct boundary id.
3963  for (cell = tria.begin(); cell != tria.end(); ++cell)
3964  for (unsigned int i = 0; i < GeometryInfo<3>::faces_per_cell; ++i)
3965  if (cell->at_boundary(i))
3966  {
3967  const Triangulation<3>::face_iterator face = cell->face(i);
3968 
3969  const Point<3> face_center(face->center());
3970  if (std::abs(face_center(0) - center(0)) >
3971  1.e-6 * face_center.norm())
3972  {
3973  if (std::abs((face_center - center).norm() - inner_radius) <
3974  std::abs((face_center - center).norm() - outer_radius))
3975  face->set_all_boundary_ids(0);
3976  else
3977  face->set_all_boundary_ids(1);
3978  }
3979  }
3980  }
3981  tria.set_all_manifold_ids(0);
3982  tria.set_manifold(0, SphericalManifold<3>(center));
3983  }
3984 
3985 
3986  // Implementation for 3D only
3987  template <>
3988  void quarter_hyper_shell(Triangulation<3> & tria,
3989  const Point<3> & center,
3990  const double inner_radius,
3991  const double outer_radius,
3992  const unsigned int n,
3993  const bool colorize)
3994  {
3995  Assert((inner_radius > 0) && (inner_radius < outer_radius),
3996  ExcInvalidRadii());
3997  if (n == 0 || n == 3)
3998  {
3999  const double a = inner_radius * std::sqrt(2.0) / 2e0;
4000  const double b = outer_radius * std::sqrt(2.0) / 2e0;
4001  const double c = a * std::sqrt(3.0) / 2e0;
4002  const double d = b * std::sqrt(3.0) / 2e0;
4003  const double e = outer_radius / 2e0;
4004  const double h = inner_radius / 2e0;
4005 
4006  std::vector<Point<3>> vertices;
4007 
4008  vertices.push_back(center + Point<3>(0, inner_radius, 0)); // 0
4009  vertices.push_back(center + Point<3>(a, a, 0)); // 1
4010  vertices.push_back(center + Point<3>(b, b, 0)); // 2
4011  vertices.push_back(center + Point<3>(0, outer_radius, 0)); // 3
4012  vertices.push_back(center + Point<3>(0, a, a)); // 4
4013  vertices.push_back(center + Point<3>(c, c, h)); // 5
4014  vertices.push_back(center + Point<3>(d, d, e)); // 6
4015  vertices.push_back(center + Point<3>(0, b, b)); // 7
4016  vertices.push_back(center + Point<3>(inner_radius, 0, 0)); // 8
4017  vertices.push_back(center + Point<3>(outer_radius, 0, 0)); // 9
4018  vertices.push_back(center + Point<3>(a, 0, a)); // 10
4019  vertices.push_back(center + Point<3>(b, 0, b)); // 11
4020  vertices.push_back(center + Point<3>(0, 0, inner_radius)); // 12
4021  vertices.push_back(center + Point<3>(0, 0, outer_radius)); // 13
4022 
4023  const int cell_vertices[3][8] = {
4024  {0, 1, 3, 2, 4, 5, 7, 6},
4025  {1, 8, 2, 9, 5, 10, 6, 11},
4026  {4, 5, 7, 6, 12, 10, 13, 11},
4027  };
4028  std::vector<CellData<3>> cells(3);
4029 
4030  for (unsigned int i = 0; i < 3; ++i)
4031  {
4032  for (unsigned int j = 0; j < 8; ++j)
4033  cells[i].vertices[j] = cell_vertices[i][j];
4034  cells[i].material_id = 0;
4035  }
4036 
4037  tria.create_triangulation(vertices,
4038  cells,
4039  SubCellData()); // no boundary information
4040  }
4041  else
4042  {
4043  AssertThrow(false, ExcNotImplemented());
4044  }
4045 
4046  if (colorize)
4047  colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius);
4048 
4049  tria.set_all_manifold_ids(0);
4050  tria.set_manifold(0, SphericalManifold<3>(center));
4051  }
4052 
4053 
4054  // Implementation for 3D only
4055  template <>
4056  void cylinder_shell(Triangulation<3> & tria,
4057  const double length,
4058  const double inner_radius,
4059  const double outer_radius,
4060  const unsigned int n_radial_cells,
4061  const unsigned int n_axial_cells)
4062  {
4063  Assert((inner_radius > 0) && (inner_radius < outer_radius),
4064  ExcInvalidRadii());
4065 
4066  const double pi = numbers::PI;
4067 
4068  // determine the number of cells
4069  // for the grid. if not provided by
4070  // the user determine it such that
4071  // the length of each cell on the
4072  // median (in the middle between
4073  // the two circles) is equal to its
4074  // radial extent (which is the
4075  // difference between the two
4076  // radii)
4077  const unsigned int N_r =
4078  (n_radial_cells == 0 ? static_cast<unsigned int>(std::ceil(
4079  (2 * pi * (outer_radius + inner_radius) / 2) /
4080  (outer_radius - inner_radius))) :
4081  n_radial_cells);
4082  const unsigned int N_z =
4083  (n_axial_cells == 0 ?
4084  static_cast<unsigned int>(std::ceil(
4085  length / (2 * pi * (outer_radius + inner_radius) / 2 / N_r))) :
4086  n_axial_cells);
4087 
4088  // set up N vertices on the
4089  // outer and N vertices on
4090  // the inner circle. the
4091  // first N ones are on the
4092  // outer one, and all are
4093  // numbered counter-clockwise
4094  std::vector<Point<2>> vertices_2d(2 * N_r);
4095  for (unsigned int i = 0; i < N_r; ++i)
4096  {
4097  vertices_2d[i] =
4098  Point<2>(std::cos(2 * pi * i / N_r), std::sin(2 * pi * i / N_r)) *
4099  outer_radius;
4100  vertices_2d[i + N_r] = vertices_2d[i] * (inner_radius / outer_radius);
4101  }
4102 
4103  std::vector<Point<3>> vertices_3d;
4104  vertices_3d.reserve(2 * N_r * (N_z + 1));
4105  for (unsigned int j = 0; j <= N_z; ++j)
4106  for (unsigned int i = 0; i < 2 * N_r; ++i)
4107  {
4108  const Point<3> v(vertices_2d[i][0],
4109  vertices_2d[i][1],
4110  j * length / N_z);
4111  vertices_3d.push_back(v);
4112  }
4113 
4114  std::vector<CellData<3>> cells(N_r * N_z, CellData<3>());
4115 
4116  for (unsigned int j = 0; j < N_z; ++j)
4117  for (unsigned int i = 0; i < N_r; ++i)
4118  {
4119  cells[i + j * N_r].vertices[0] = i + (j + 1) * 2 * N_r;
4120  cells[i + j * N_r].vertices[1] = (i + 1) % N_r + (j + 1) * 2 * N_r;
4121  cells[i + j * N_r].vertices[2] = i + j * 2 * N_r;
4122  cells[i + j * N_r].vertices[3] = (i + 1) % N_r + j * 2 * N_r;
4123 
4124  cells[i + j * N_r].vertices[4] = N_r + i + (j + 1) * 2 * N_r;
4125  cells[i + j * N_r].vertices[5] =
4126  N_r + ((i + 1) % N_r) + (j + 1) * 2 * N_r;
4127  cells[i + j * N_r].vertices[6] = N_r + i + j * 2 * N_r;
4128  cells[i + j * N_r].vertices[7] = N_r + ((i + 1) % N_r) + j * 2 * N_r;
4129 
4130  cells[i + j * N_r].material_id = 0;
4131  }
4132 
4133  tria.create_triangulation(vertices_3d, cells, SubCellData());
4134  tria.set_all_manifold_ids(0);
4135  tria.set_manifold(0, CylindricalManifold<3>(2));
4136  }
4137 
4138 
4139 
4140  template <int dim, int spacedim>
4141  void
4142  merge_triangulations(
4143  const std::initializer_list<const Triangulation<dim, spacedim> *const>
4144  triangulations,
4145  Triangulation<dim, spacedim> &result,
4146  const double duplicated_vertex_tolerance,
4147  const bool copy_manifold_ids)
4148  {
4149  for (const auto triangulation : triangulations)
4150  {
4151  (void)triangulation;
4152  Assert(triangulation->n_levels() == 1,
4153  ExcMessage("The input triangulations must be non-empty "
4154  "and must not be refined."));
4155  }
4156 
4157  // get the union of the set of vertices
4158  unsigned int n_vertices = 0;
4159  for (const auto triangulation : triangulations)
4160  n_vertices += triangulation->n_vertices();
4161 
4162  std::vector<Point<spacedim>> vertices;
4163  vertices.reserve(n_vertices);
4164  for (const auto triangulation : triangulations)
4165  vertices.insert(vertices.end(),
4166  triangulation->get_vertices().begin(),
4167  triangulation->get_vertices().end());
4168 
4169  // now form the union of the set of cells
4170  std::vector<CellData<dim>> cells;
4171  unsigned int n_cells = 0;
4172  for (const auto triangulation : triangulations)
4173  n_cells += triangulation->n_cells();
4174 
4175  cells.reserve(n_cells);
4176  unsigned int n_accumulated_vertices = 0;
4177  for (const auto triangulation : triangulations)
4178  {
4179  for (const auto &cell : triangulation->cell_iterators())
4180  {
4181  CellData<dim> this_cell;
4182  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
4183  ++v)
4184  this_cell.vertices[v] =
4185  cell->vertex_index(v) + n_accumulated_vertices;
4186  this_cell.material_id = cell->material_id();
4187  this_cell.manifold_id = cell->manifold_id();
4188  cells.push_back(std::move(this_cell));
4189  }
4190  n_accumulated_vertices += triangulation->n_vertices();
4191  }
4192 
4193  // Now for the SubCellData
4194  SubCellData subcell_data;
4195  if (copy_manifold_ids)
4196  {
4197  switch (dim)
4198  {
4199  case 1:
4200  break;
4201 
4202  case 2:
4203  {
4204  unsigned int n_boundary_lines = 0;
4205  for (const auto triangulation : triangulations)
4206  n_boundary_lines += triangulation->n_lines();
4207 
4208  subcell_data.boundary_lines.reserve(n_boundary_lines);
4209 
4210  auto copy_line_manifold_ids =
4211  [&](const Triangulation<dim, spacedim> &tria,
4212  const unsigned int offset) {
4213  for (typename Triangulation<dim, spacedim>::line_iterator
4214  line = tria.begin_face();
4215  line != tria.end_face();
4216  ++line)
4217  if (line->manifold_id() != numbers::flat_manifold_id)
4218  {
4219  CellData<1> boundary_line;
4220  boundary_line.vertices[0] =
4221  line->vertex_index(0) + offset;
4222  boundary_line.vertices[1] =
4223  line->vertex_index(1) + offset;
4224  boundary_line.manifold_id = line->manifold_id();
4225  subcell_data.boundary_lines.push_back(
4226  boundary_line); // trivially-copyable
4227  }
4228  };
4229 
4230  unsigned int n_accumulated_vertices = 0;
4231  for (const auto triangulation : triangulations)
4232  {
4233  copy_line_manifold_ids(*triangulation,
4234  n_accumulated_vertices);
4235  n_accumulated_vertices += triangulation->n_vertices();
4236  }
4237  break;
4238  }
4239 
4240  case 3:
4241  {
4242  unsigned int n_boundary_quads = 0;
4243  for (const auto triangulation : triangulations)
4244  n_boundary_quads += triangulation->n_quads();
4245 
4246  subcell_data.boundary_quads.reserve(n_boundary_quads);
4247  // we can't do better here than to loop over all the lines
4248  // bounding a face. For regular meshes an (interior) line in
4249  // 3D is part of four cells. So this should be an appropriate
4250  // guess.
4251  subcell_data.boundary_lines.reserve(4 * n_boundary_quads);
4252 
4253  auto copy_face_and_line_manifold_ids =
4254  [&](const Triangulation<dim, spacedim> &tria,
4255  const unsigned int offset) {
4256  for (typename Triangulation<dim, spacedim>::face_iterator
4257  face = tria.begin_face();
4258  face != tria.end_face();
4259  ++face)
4260  if (face->manifold_id() != numbers::flat_manifold_id)
4261  {
4262  CellData<2> boundary_quad;
4263  boundary_quad.vertices[0] =
4264  face->vertex_index(0) + offset;
4265  boundary_quad.vertices[1] =
4266  face->vertex_index(1) + offset;
4267  boundary_quad.vertices[2] =
4268  face->vertex_index(2) + offset;
4269  boundary_quad.vertices[3] =
4270  face->vertex_index(3) + offset;
4271 
4272  boundary_quad.manifold_id = face->manifold_id();
4273  subcell_data.boundary_quads.push_back(
4274  boundary_quad); // trivially-copyable
4275  }
4276  for (const auto &cell : tria.cell_iterators())
4277  for (unsigned int l = 0; l < 12; ++l)
4278  if (cell->line(l)->manifold_id() !=
4279  numbers::flat_manifold_id)
4280  {
4281  CellData<1> boundary_line;
4282  boundary_line.vertices[0] =
4283  cell->line(l)->vertex_index(0) + offset;
4284  boundary_line.vertices[1] =
4285  cell->line(l)->vertex_index(1) + offset;
4286  boundary_line.manifold_id =
4287  cell->line(l)->manifold_id();
4288  subcell_data.boundary_lines.push_back(
4289  boundary_line); // trivially_copyable
4290  }
4291  };
4292  unsigned int n_accumulated_vertices = 0;
4293  for (const auto triangulation : triangulations)
4294  {
4295  copy_face_and_line_manifold_ids(*triangulation,
4296  n_accumulated_vertices);
4297  n_accumulated_vertices += triangulation->n_vertices();
4298  }
4299  break;
4300  }
4301  default:
4302  Assert(false, ExcNotImplemented());
4303  }
4304  }
4305 
4306  // throw out duplicated vertices
4307  std::vector<unsigned int> considered_vertices;
4308  GridTools::delete_duplicated_vertices(vertices,
4309  cells,
4310  subcell_data,
4311  considered_vertices,
4312  duplicated_vertex_tolerance);
4313 
4314  // reorder the cells to ensure that they satisfy the convention for
4315  // edge and face directions
4316  GridReordering<dim, spacedim>::reorder_cells(cells, true);
4317  result.clear();
4318  result.create_triangulation(vertices, cells, subcell_data);
4319  }
4320 
4321 
4322 
4323  template <int dim, int spacedim>
4324  void
4325  merge_triangulations(const Triangulation<dim, spacedim> &triangulation_1,
4326  const Triangulation<dim, spacedim> &triangulation_2,
4327  Triangulation<dim, spacedim> & result,
4328  const double duplicated_vertex_tolerance,
4329  const bool copy_manifold_ids)
4330  {
4331  // if either Triangulation is empty then merging is just a copy.
4332  if (triangulation_1.n_cells() == 0)
4333  {
4334  result.copy_triangulation(triangulation_2);
4335  return;
4336  }
4337  if (triangulation_2.n_cells() == 0)
4338  {
4339  result.copy_triangulation(triangulation_1);
4340  return;
4341  }
4342  merge_triangulations({&triangulation_1, &triangulation_2},
4343  result,
4344  duplicated_vertex_tolerance,
4345  copy_manifold_ids);
4346  }
4347 
4348 
4349 
4350  template <int dim, int spacedim>
4351  void
4352  create_union_triangulation(
4353  const Triangulation<dim, spacedim> &triangulation_1,
4354  const Triangulation<dim, spacedim> &triangulation_2,
4355  Triangulation<dim, spacedim> & result)
4356  {
4357  Assert(GridTools::have_same_coarse_mesh(triangulation_1, triangulation_2),
4358  ExcMessage("The two input triangulations are not derived from "
4359  "the same coarse mesh as required."));
4360  Assert((dynamic_cast<
4361  const parallel::distributed::Triangulation<dim, spacedim> *>(
4362  &triangulation_1) == nullptr) &&
4363  (dynamic_cast<
4364  const parallel::distributed::Triangulation<dim, spacedim> *>(
4365  &triangulation_2) == nullptr),
4366  ExcMessage("The source triangulations for this function must both "
4367  "be available entirely locally, and not be distributed "
4368  "triangulations."));
4369 
4370  // first copy triangulation_1, and
4371  // then do as many iterations as
4372  // there are levels in
4373  // triangulation_2 to refine
4374  // additional cells. since this is
4375  // the maximum number of
4376  // refinements to get from the
4377  // coarse grid to triangulation_2,
4378  // it is clear that this is also
4379  // the maximum number of
4380  // refinements to get from any cell
4381  // on triangulation_1 to
4382  // triangulation_2
4383  result.clear();
4384  result.copy_triangulation(triangulation_1);
4385  for (unsigned int iteration = 0; iteration < triangulation_2.n_levels();
4386  ++iteration)
4387  {
4388  InterGridMap<Triangulation<dim, spacedim>> intergrid_map;
4389  intergrid_map.make_mapping(result, triangulation_2);
4390 
4391  bool any_cell_flagged = false;
4392  for (typename Triangulation<dim, spacedim>::active_cell_iterator
4393  result_cell = result.begin_active();
4394  result_cell != result.end();
4395  ++result_cell)
4396  if (intergrid_map[result_cell]->has_children())
4397  {
4398  any_cell_flagged = true;
4399  result_cell->set_refine_flag();
4400  }
4401 
4402  if (any_cell_flagged == false)
4403  break;
4404  else
4405  result.execute_coarsening_and_refinement();
4406  }
4407  }
4408 
4409 
4410 
4411  template <int dim, int spacedim>
4412  void
4413  create_triangulation_with_removed_cells(
4414  const Triangulation<dim, spacedim> &input_triangulation,
4415  const std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>
4416  & cells_to_remove,
4417  Triangulation<dim, spacedim> &result)
4418  {
4419  // simply copy the vertices; we will later strip those
4420  // that turn out to be unused
4421  std::vector<Point<spacedim>> vertices = input_triangulation.get_vertices();
4422 
4423  // the loop through the cells and copy stuff, excluding
4424  // the ones we are to remove
4425  std::vector<CellData<dim>> cells;
4426  for (typename Triangulation<dim, spacedim>::active_cell_iterator cell =
4427  input_triangulation.begin_active();
4428  cell != input_triangulation.end();
4429  ++cell)
4430  if (cells_to_remove.find(cell) == cells_to_remove.end())
4431  {
4432  Assert(static_cast<unsigned int>(cell->level()) ==
4433  input_triangulation.n_levels() - 1,
4434  ExcMessage(
4435  "Your input triangulation appears to have "
4436  "adaptively refined cells. This is not allowed. You can "
4437  "only call this function on a triangulation in which "
4438  "all cells are on the same refinement level."));
4439 
4440  CellData<dim> this_cell;
4441  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
4442  ++v)
4443  this_cell.vertices[v] = cell->vertex_index(v);
4444  this_cell.material_id = cell->material_id();
4445  cells.push_back(this_cell);
4446  }
4447 
4448  // throw out duplicated vertices from the two meshes, reorder vertices as
4449  // necessary and create the triangulation
4450  SubCellData subcell_data;
4451  std::vector<unsigned int> considered_vertices;
4452  GridTools::delete_duplicated_vertices(vertices,
4453  cells,
4454  subcell_data,
4455  considered_vertices);
4456 
4457  // then clear the old triangulation and create the new one
4458  result.clear();
4459  result.create_triangulation(vertices, cells, subcell_data);
4460  }
4461 
4462 
4463 
4464  void
4465  extrude_triangulation(const Triangulation<2, 2> &input,
4466  const unsigned int n_slices,
4467  const double height,
4468  Triangulation<3, 3> & result)
4469  {
4470  Assert(input.n_levels() == 1,
4471  ExcMessage(
4472  "The input triangulation must be a coarse mesh, i.e., it must "
4473  "not have been refined."));
4474  Assert(result.n_cells() == 0,
4475  ExcMessage("The output triangulation object needs to be empty."));
4476  Assert(height > 0,
4477  ExcMessage("The given height for extrusion must be positive."));
4478  Assert(n_slices >= 2,
4479  ExcMessage(
4480  "The number of slices for extrusion must be at least 2."));
4481 
4482  const double delta_h = height / (n_slices - 1);
4483  std::vector<double> slices_z_values;
4484  for (unsigned int i = 0; i < n_slices; ++i)
4485  slices_z_values.push_back(i * delta_h);
4486  extrude_triangulation(input, slices_z_values, result);
4487  }
4488 
4489  void
4490  extrude_triangulation(const Triangulation<2, 2> &input,
4491  const std::vector<double> &slice_coordinates,
4492  Triangulation<3, 3> & result)
4493  {
4494  Assert(input.n_levels() == 1,
4495  ExcMessage(
4496  "The input triangulation must be a coarse mesh, i.e., it must "
4497  "not have been refined."));
4498  Assert(result.n_cells() == 0,
4499  ExcMessage("The output triangulation object needs to be empty."));
4500  Assert(slice_coordinates.size() >= 2,
4501  ExcMessage(
4502  "The number of slices for extrusion must be at least 2."));
4503  const unsigned int n_slices = slice_coordinates.size();
4504  Assert(std::is_sorted(slice_coordinates.begin(), slice_coordinates.end()),
4505  ExcMessage("Slice z-coordinates should be in ascending order"));
4506  std::vector<Point<3>> points(n_slices * input.n_vertices());
4507  std::vector<CellData<3>> cells;
4508  cells.reserve((n_slices - 1) * input.n_active_cells());
4509 
4510  // copy the array of points as many times as there will be slices,
4511  // one slice at a time. The z-axis value are defined in slices_coordinates
4512  for (unsigned int slice = 0; slice < n_slices; ++slice)
4513  {
4514  for (unsigned int i = 0; i < input.n_vertices(); ++i)
4515  {
4516  const Point<2> &v = input.get_vertices()[i];
4517  points[slice * input.n_vertices() + i](0) = v(0);
4518  points[slice * input.n_vertices() + i](1) = v(1);
4519  points[slice * input.n_vertices() + i](2) =
4520  slice_coordinates[slice];
4521  }
4522  }
4523 
4524  // then create the cells of each of the slices, one stack at a
4525  // time
4526  for (Triangulation<2, 2>::cell_iterator cell = input.begin();
4527  cell != input.end();
4528  ++cell)
4529  {
4530  for (unsigned int slice = 0; slice < n_slices - 1; ++slice)
4531  {
4532  CellData<3> this_cell;
4533  for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell;
4534  ++v)
4535  {
4536  this_cell.vertices[v] =
4537  cell->vertex_index(v) + slice * input.n_vertices();
4538  this_cell.vertices[v + GeometryInfo<2>::vertices_per_cell] =
4539  cell->vertex_index(v) + (slice + 1) * input.n_vertices();
4540  }
4541 
4542  this_cell.material_id = cell->material_id();
4543  cells.push_back(this_cell);
4544  }
4545  }
4546 
4547  // next, create face data for all of the outer faces for which the
4548  // boundary indicator will not be equal to zero (where we would
4549  // explicitly set it to something that is already the default --
4550  // no need to do that)
4551  SubCellData s;
4552  types::boundary_id max_boundary_id = 0;
4553  s.boundary_quads.reserve(input.n_active_lines() * (n_slices - 1) +
4554  input.n_active_cells() * 2);
4555  for (Triangulation<2, 2>::cell_iterator cell = input.begin();
4556  cell != input.end();
4557  ++cell)
4558  {
4559  CellData<2> quad;
4560  for (unsigned int f = 0; f < 4; ++f)
4561  if (cell->at_boundary(f) && (cell->face(f)->boundary_id() != 0))
4562  {
4563  quad.boundary_id = cell->face(f)->boundary_id();
4564  max_boundary_id = std::max(max_boundary_id, quad.boundary_id);
4565  for (unsigned int slice = 0; slice < n_slices - 1; ++slice)
4566  {
4567  quad.vertices[0] =
4568  cell->face(f)->vertex_index(0) + slice * input.n_vertices();
4569  quad.vertices[1] =
4570  cell->face(f)->vertex_index(1) + slice * input.n_vertices();
4571  quad.vertices[2] = cell->face(f)->vertex_index(0) +
4572  (slice + 1) * input.n_vertices();
4573  quad.vertices[3] = cell->face(f)->vertex_index(1) +
4574  (slice + 1) * input.n_vertices();
4575  s.boundary_quads.push_back(quad);
4576  }
4577  }
4578  }
4579 
4580  // then mark the bottom and top boundaries of the extruded mesh
4581  // with max_boundary_id+1 and max_boundary_id+2. check that this
4582  // remains valid
4583  Assert((max_boundary_id != numbers::invalid_boundary_id) &&
4584  (max_boundary_id + 1 != numbers::invalid_boundary_id) &&
4585  (max_boundary_id + 2 != numbers::invalid_boundary_id),
4586  ExcMessage(
4587  "The input triangulation to this function is using boundary "
4588  "indicators in a range that do not allow using "
4589  "max_boundary_id+1 and max_boundary_id+2 as boundary "
4590  "indicators for the bottom and top faces of the "
4591  "extruded triangulation."));
4592  for (Triangulation<2, 2>::cell_iterator cell = input.begin();
4593  cell != input.end();
4594  ++cell)
4595  {
4596  CellData<2> quad;
4597  quad.boundary_id = max_boundary_id + 1;
4598  quad.vertices[0] = cell->vertex_index(0);
4599  quad.vertices[1] = cell->vertex_index(1);
4600  quad.vertices[2] = cell->vertex_index(2);
4601  quad.vertices[3] = cell->vertex_index(3);
4602  s.boundary_quads.push_back(quad);
4603 
4604  quad.boundary_id = max_boundary_id + 2;
4605  for (int i = 0; i < 4; ++i)
4606  quad.vertices[i] += (n_slices - 1) * input.n_vertices();
4607  s.boundary_quads.push_back(quad);
4608  }
4609 
4610  // use all of this to finally create the extruded 3d
4611  // triangulation. it is not necessary to call
4612  // GridReordering<3,3>::reorder_cells because the cells we have
4613  // constructed above are automatically correctly oriented. this is
4614  // because the 2d base mesh is always correctly oriented, and
4615  // extruding it automatically yields a correctly oriented 3d mesh,
4616  // as discussed in the edge orientation paper mentioned in the
4617  // introduction to the GridReordering class.
4618  result.create_triangulation(points, cells, s);
4619  }
4620 
4621 
4622  template <>
4623  void hyper_cube_with_cylindrical_hole(Triangulation<1> &,
4624  const double,
4625  const double,
4626  const double,
4627  const unsigned int,
4628  const bool)
4629  {
4630  Assert(false, ExcNotImplemented());
4631  }
4632 
4633 
4634 
4635  template <>
4636  void hyper_cube_with_cylindrical_hole(Triangulation<2> &triangulation,
4637  const double inner_radius,
4638  const double outer_radius,
4639  const double, // width,
4640  const unsigned int, // width_repetition,
4641  const bool colorize)
4642  {
4643  const int dim = 2;
4644 
4645  Assert(inner_radius < outer_radius,
4646  ExcMessage("outer_radius has to be bigger than inner_radius."));
4647 
4648  Point<dim> center;
4649  // We create an hyper_shell in two dimensions, and then we modify it.
4650  hyper_shell(triangulation, center, inner_radius, outer_radius, 8);
4651  triangulation.set_all_manifold_ids(numbers::flat_manifold_id);
4652  Triangulation<dim>::active_cell_iterator cell =
4653  triangulation.begin_active(),
4654  endc = triangulation.end();
4655  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
4656  for (; cell != endc; ++cell)
4657  {
4658  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
4659  if (cell->face(f)->at_boundary())
4660  {
4661  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
4662  ++v)
4663  {
4664  unsigned int vv = cell->face(f)->vertex_index(v);
4665  if (treated_vertices[vv] == false)
4666  {
4667  treated_vertices[vv] = true;
4668  switch (vv)
4669  {
4670  case 1:
4671  cell->face(f)->vertex(v) =
4672  center + Point<dim>(outer_radius, outer_radius);
4673  break;
4674  case 3:
4675  cell->face(f)->vertex(v) =
4676  center + Point<dim>(-outer_radius, outer_radius);
4677  break;
4678  case 5:
4679  cell->face(f)->vertex(v) =
4680  center + Point<dim>(-outer_radius, -outer_radius);
4681  break;
4682  case 7:
4683  cell->face(f)->vertex(v) =
4684  center + Point<dim>(outer_radius, -outer_radius);
4685  break;
4686  default:
4687  break;
4688  }
4689  }
4690  }
4691  }
4692  }
4693  double eps = 1e-3 * outer_radius;
4694  cell = triangulation.begin_active();
4695  for (; cell != endc; ++cell)
4696  {
4697  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
4698  if (cell->face(f)->at_boundary())
4699  {
4700  double dx = cell->face(f)->center()(0) - center(0);
4701  double dy = cell->face(f)->center()(1) - center(1);
4702  if (colorize)
4703  {
4704  if (std::abs(dx + outer_radius) < eps)
4705  cell->face(f)->set_boundary_id(0);
4706  else if (std::abs(dx - outer_radius) < eps)
4707  cell->face(f)->set_boundary_id(1);
4708  else if (std::abs(dy + outer_radius) < eps)
4709  cell->face(f)->set_boundary_id(2);
4710  else if (std::abs(dy - outer_radius) < eps)
4711  cell->face(f)->set_boundary_id(3);
4712  else
4713  {
4714  cell->face(f)->set_boundary_id(4);
4715  cell->face(f)->set_manifold_id(0);
4716  }
4717  }
4718  else
4719  {
4720  double d = (cell->face(f)->center() - center).norm();
4721  if (d - inner_radius < 0)
4722  {
4723  cell->face(f)->set_boundary_id(1);
4724  cell->face(f)->set_manifold_id(0);
4725  }
4726  else
4727  cell->face(f)->set_boundary_id(0);
4728  }
4729  }
4730  }
4731  triangulation.set_manifold(0, PolarManifold<2>(center));
4732  }
4733 
4734 
4735 
4736  template <int dim>
4737  void
4738  concentric_hyper_shells(Triangulation<dim> &triangulation,
4739  const Point<dim> & center,
4740  const double inner_radius,
4741  const double outer_radius,
4742  const unsigned int n_shells,
4743  const double skewness,
4744  const unsigned int n_cells,
4745  const bool colorize)
4746  {
4747  Assert(dim == 2 || dim == 3, ExcNotImplemented());
4748  (void)colorize;
4749  (void)n_cells;
4750  Assert(inner_radius < outer_radius,
4751  ExcMessage("outer_radius has to be bigger than inner_radius."));
4752  if (n_shells == 0)
4753  return; // empty Triangulation
4754 
4755  std::vector<double> radii;
4756  radii.push_back(inner_radius);
4757  for (unsigned int shell_n = 1; shell_n < n_shells; ++shell_n)
4758  if (skewness == 0.0)
4759  // same as below, but works in the limiting case of zero skewness
4760  radii.push_back(inner_radius +
4761  (outer_radius - inner_radius) *
4762  (1.0 - (1.0 - double(shell_n) / n_shells)));
4763  else
4764  radii.push_back(
4765  inner_radius +
4766  (outer_radius - inner_radius) *
4767  (1.0 - std::tanh(skewness * (1.0 - double(shell_n) / n_shells)) /
4768  std::tanh(skewness)));
4769  radii.push_back(outer_radius);
4770 
4771  auto min_line_length = [](const Triangulation<dim> &tria) -> double {
4772  double length = std::numeric_limits<double>::max();
4773  for (const auto &cell : tria.active_cell_iterators())
4774  for (unsigned int n = 0; n < GeometryInfo<dim>::lines_per_cell; ++n)
4775  length = std::min(length, cell->line(n)->diameter());
4776  return length;
4777  };
4778 
4779  double grid_vertex_tolerance = 0.0;
4780  for (unsigned int shell_n = 0; shell_n < radii.size() - 1; ++shell_n)
4781  {
4782  Triangulation<dim> current_shell;
4783  GridGenerator::hyper_shell(current_shell,
4784  center,
4785  radii[shell_n],
4786  radii[shell_n + 1],
4787  n_cells == 0 ? (dim == 2 ? 8 : 12) :
4788  n_cells);
4789 
4790  // The innermost shell has the smallest cells: use that to set the
4791  // vertex merging tolerance
4792  if (grid_vertex_tolerance == 0.0)
4793  grid_vertex_tolerance = 0.5 * min_line_length(current_shell);
4794 
4795  Triangulation<dim> temp(std::move(triangulation));
4796  triangulation.clear();
4797  GridGenerator::merge_triangulations(current_shell,
4798  temp,
4799  triangulation,
4800  grid_vertex_tolerance);
4801  }
4802 
4803  const types::manifold_id manifold_id = 0;
4804  triangulation.set_all_manifold_ids(manifold_id);
4805  if (dim == 2)
4806  triangulation.set_manifold(manifold_id, PolarManifold<dim>(center));
4807  else if (dim == 3)
4808  triangulation.set_manifold(manifold_id, SphericalManifold<dim>(center));
4809 
4810  // We use boundary vertex positions to see if things are on the inner or
4811  // outer boundary.
4812  constexpr double radial_vertex_tolerance =
4813  100.0 * std::numeric_limits<double>::epsilon();
4814  auto assert_vertex_distance_within_tolerance =
4815  [center, radial_vertex_tolerance](
4816  const TriaIterator<TriaAccessor<dim - 1, dim, dim>> face,
4817  const double radius) {
4818  (void)center;
4819  (void)radial_vertex_tolerance;
4820  (void)face;
4821  (void)radius;
4822  for (unsigned int vertex_n = 0;
4823  vertex_n < GeometryInfo<dim>::vertices_per_face;
4824  ++vertex_n)
4825  {
4826  Assert(std::abs((face->vertex(vertex_n) - center).norm() - radius) <
4827  (center.norm() + radius) * radial_vertex_tolerance,
4828  ExcInternalError());
4829  }
4830  };
4831  if (colorize)
4832  for (const auto &cell : triangulation.active_cell_iterators())
4833  for (unsigned int face_n = 0;
4834  face_n < GeometryInfo<dim>::faces_per_cell;
4835  ++face_n)
4836  {
4837  auto face = cell->face(face_n);
4838  if (face->at_boundary())
4839  {
4840  if (((face->vertex(0) - center).norm() - inner_radius) <
4841  (center.norm() + inner_radius) * radial_vertex_tolerance)
4842  {
4843  // we must be at an inner face, but check
4844  assert_vertex_distance_within_tolerance(face, inner_radius);
4845  face->set_all_boundary_ids(0);
4846  }
4847  else
4848  {
4849  // we must be at an outer face, but check
4850  assert_vertex_distance_within_tolerance(face, outer_radius);
4851  face->set_all_boundary_ids(1);
4852  }
4853  }
4854  }
4855  }
4856 
4857 
4858 
4859  template <>
4860  void hyper_cube_with_cylindrical_hole(Triangulation<3> & triangulation,
4861  const double inner_radius,
4862  const double outer_radius,
4863  const double L,
4864  const unsigned int Nz,
4865  const bool colorize)
4866  {
4867  const int dim = 3;
4868 
4869  Assert(inner_radius < outer_radius,
4870  ExcMessage("outer_radius has to be bigger than inner_radius."));
4871  Assert(L > 0, ExcMessage("Must give positive extension L"));
4872  Assert(Nz >= 1, ExcLowerRange(1, Nz));
4873 
4874  cylinder_shell(triangulation, L, inner_radius, outer_radius, 8, Nz);
4875  triangulation.set_all_manifold_ids(numbers::flat_manifold_id);
4876 
4877  Triangulation<dim>::active_cell_iterator cell =
4878  triangulation.begin_active(),
4879  endc = triangulation.end();
4880  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
4881  for (; cell != endc; ++cell)
4882  {
4883  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
4884  if (cell->face(f)->at_boundary())
4885  {
4886  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
4887  ++v)
4888  {
4889  unsigned int vv = cell->face(f)->vertex_index(v);
4890  if (treated_vertices[vv] == false)
4891  {
4892  treated_vertices[vv] = true;
4893  for (unsigned int i = 0; i <= Nz; ++i)
4894  {
4895  double d = ((double)i) * L / ((double)Nz);
4896  switch (vv - i * 16)
4897  {
4898  case 1:
4899  cell->face(f)->vertex(v) =
4900  Point<dim>(outer_radius, outer_radius, d);
4901  break;
4902  case 3:
4903  cell->face(f)->vertex(v) =
4904  Point<dim>(-outer_radius, outer_radius, d);
4905  break;
4906  case 5:
4907  cell->face(f)->vertex(v) =
4908  Point<dim>(-outer_radius, -outer_radius, d);
4909  break;
4910  case 7:
4911  cell->face(f)->vertex(v) =
4912  Point<dim>(outer_radius, -outer_radius, d);
4913  break;
4914  default:
4915  break;
4916  }
4917  }
4918  }
4919  }
4920  }
4921  }
4922  double eps = 1e-3 * outer_radius;
4923  cell = triangulation.begin_active();
4924  for (; cell != endc; ++cell)
4925  {
4926  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
4927  if (cell->face(f)->at_boundary())
4928  {
4929  double dx = cell->face(f)->center()(0);
4930  double dy = cell->face(f)->center()(1);
4931  double dz = cell->face(f)->center()(2);
4932 
4933  if (colorize)
4934  {
4935  if (std::abs(dx + outer_radius) < eps)
4936  cell->face(f)->set_boundary_id(0);
4937 
4938  else if (std::abs(dx - outer_radius) < eps)
4939  cell->face(f)->set_boundary_id(1);
4940 
4941  else if (std::abs(dy + outer_radius) < eps)
4942  cell->face(f)->set_boundary_id(2);
4943 
4944  else if (std::abs(dy - outer_radius) < eps)
4945  cell->face(f)->set_boundary_id(3);
4946 
4947  else if (std::abs(dz) < eps)
4948  cell->face(f)->set_boundary_id(4);
4949 
4950  else if (std::abs(dz - L) < eps)
4951  cell->face(f)->set_boundary_id(5);
4952 
4953  else
4954  {
4955  cell->face(f)->set_all_boundary_ids(6);
4956  cell->face(f)->set_all_manifold_ids(0);
4957  }
4958  }
4959  else
4960  {
4961  Point<dim> c = cell->face(f)->center();
4962  c(2) = 0;
4963  double d = c.norm();
4964  if (d - inner_radius < 0)
4965  {
4966  cell->face(f)->set_all_boundary_ids(1);
4967  cell->face(f)->set_all_manifold_ids(0);
4968  }
4969  else
4970  cell->face(f)->set_boundary_id(0);
4971  }
4972  }
4973  }
4974  triangulation.set_manifold(0, CylindricalManifold<3>(2));
4975  }
4976 
4977  template <int dim, int spacedim1, int spacedim2>
4978  void
4979  flatten_triangulation(const Triangulation<dim, spacedim1> &in_tria,
4980  Triangulation<dim, spacedim2> & out_tria)
4981  {
4982  const parallel::distributed::Triangulation<dim, spacedim1> *pt =
4983  dynamic_cast<
4984  const parallel::distributed::Triangulation<dim, spacedim1> *>(&in_tria);
4985 
4986  (void)pt;
4987  Assert(
4988  pt == nullptr,
4989  ExcMessage(
4990  "Cannot use this function on parallel::distributed::Triangulation."));
4991 
4992  std::vector<Point<spacedim2>> v;
4993  std::vector<CellData<dim>> cells;
4994  SubCellData subcelldata;
4995 
4996  const unsigned int spacedim = std::min(spacedim1, spacedim2);
4997  const std::vector<Point<spacedim1>> &in_vertices = in_tria.get_vertices();
4998 
4999  v.resize(in_vertices.size());
5000  for (unsigned int i = 0; i < in_vertices.size(); ++i)
5001  for (unsigned int d = 0; d < spacedim; ++d)
5002  v[i][d] = in_vertices[i][d];
5003 
5004  cells.resize(in_tria.n_active_cells());
5005  typename Triangulation<dim, spacedim1>::active_cell_iterator
5006  cell = in_tria.begin_active(),
5007  endc = in_tria.end();
5008 
5009  for (unsigned int id = 0; cell != endc; ++cell, ++id)
5010  {
5011  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
5012  cells[id].vertices[i] = cell->vertex_index(i);
5013  cells[id].material_id = cell->material_id();
5014  cells[id].manifold_id = cell->manifold_id();
5015  }
5016 
5017  if (dim > 1)
5018  {
5019  typename Triangulation<dim, spacedim1>::active_face_iterator
5020  face = in_tria.begin_active_face(),
5021  endf = in_tria.end_face();
5022 
5023  // Face counter for both dim == 2 and dim == 3
5024  unsigned int f = 0;
5025  switch (dim)
5026  {
5027  case 2:
5028  {
5029  subcelldata.boundary_lines.resize(in_tria.n_active_faces());
5030  for (; face != endf; ++face)
5031  if (face->at_boundary())
5032  {
5033  for (unsigned int i = 0;
5034  i < GeometryInfo<dim>::vertices_per_face;
5035  ++i)
5036  subcelldata.boundary_lines[f].vertices[i] =
5037  face->vertex_index(i);
5038  subcelldata.boundary_lines[f].boundary_id =
5039  face->boundary_id();
5040  subcelldata.boundary_lines[f].manifold_id =
5041  face->manifold_id();
5042  ++f;
5043  }
5044  subcelldata.boundary_lines.resize(f);
5045  }
5046  break;
5047  case 3:
5048  {
5049  subcelldata.boundary_quads.resize(in_tria.n_active_faces());
5050  for (; face != endf; ++face)
5051  if (face->at_boundary())
5052  {
5053  for (unsigned int i = 0;
5054  i < GeometryInfo<dim>::vertices_per_face;
5055  ++i)
5056  subcelldata.boundary_quads[f].vertices[i] =
5057  face->vertex_index(i);
5058  subcelldata.boundary_quads[f].boundary_id =
5059  face->boundary_id();
5060  subcelldata.boundary_quads[f].manifold_id =
5061  face->manifold_id();
5062  ++f;
5063  }
5064  subcelldata.boundary_quads.resize(f);
5065  }
5066  break;
5067  default:
5068  Assert(false, ExcInternalError());
5069  }
5070  }
5071  out_tria.create_triangulation(v, cells, subcelldata);
5072  }
5073 
5074 
5075 
5076  template <template <int, int> class MeshType, int dim, int spacedim>
5077 #ifndef _MSC_VER
5078  std::map<typename MeshType<dim - 1, spacedim>::cell_iterator,
5079  typename MeshType<dim, spacedim>::face_iterator>
5080 #else
5081  typename ExtractBoundaryMesh<MeshType, dim, spacedim>::return_type
5082 #endif
5083  extract_boundary_mesh(const MeshType<dim, spacedim> & volume_mesh,
5084  MeshType<dim - 1, spacedim> & surface_mesh,
5085  const std::set<types::boundary_id> &boundary_ids)
5086  {
5087  Assert((dynamic_cast<
5088  const parallel::distributed::Triangulation<dim, spacedim> *>(
5089  &volume_mesh.get_triangulation()) == nullptr),
5090  ExcNotImplemented());
5091 
5092  // This function works using the following assumption:
5093  // Triangulation::create_triangulation(...) will create cells that
5094  // preserve the order of cells passed in using the CellData argument;
5095  // also, that it will not reorder the vertices.
5096 
5097  std::map<typename MeshType<dim - 1, spacedim>::cell_iterator,
5098  typename MeshType<dim, spacedim>::face_iterator>
5099  surface_to_volume_mapping;
5100 
5101  const unsigned int boundary_dim = dim - 1; // dimension of the boundary mesh
5102 
5103  // First create surface mesh and mapping
5104  // from only level(0) cells of volume_mesh
5105  std::vector<typename MeshType<dim, spacedim>::face_iterator>
5106  mapping; // temporary map for level==0
5107 
5108 
5109  std::vector<bool> touched(volume_mesh.get_triangulation().n_vertices(),
5110  false);
5111  std::vector<CellData<boundary_dim>> cells;
5112  SubCellData subcell_data;
5113  std::vector<Point<spacedim>> vertices;
5114 
5115  std::map<unsigned int, unsigned int>
5116  map_vert_index; // volume vertex indices to surf ones
5117 
5118  for (typename MeshType<dim, spacedim>::cell_iterator cell =
5119  volume_mesh.begin(0);
5120  cell != volume_mesh.end(0);
5121  ++cell)
5122  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
5123  {
5124  const typename MeshType<dim, spacedim>::face_iterator face =
5125  cell->face(i);
5126 
5127  if (face->at_boundary() &&
5128  (boundary_ids.empty() ||
5129  (boundary_ids.find(face->boundary_id()) != boundary_ids.end())))
5130  {
5131  CellData<boundary_dim> c_data;
5132 
5133  for (unsigned int j = 0;
5134  j < GeometryInfo<boundary_dim>::vertices_per_cell;
5135  ++j)
5136  {
5137  const unsigned int v_index = face->vertex_index(j);
5138 
5139  if (!touched[v_index])
5140  {
5141  vertices.push_back(face->vertex(j));
5142  map_vert_index[v_index] = vertices.size() - 1;
5143  touched[v_index] = true;
5144  }
5145 
5146  c_data.vertices[j] = map_vert_index[v_index];
5147  c_data.material_id =
5148  static_cast<types::material_id>(face->boundary_id());
5149  c_data.manifold_id = face->manifold_id();
5150  }
5151 
5152  // if we start from a 3d mesh, then we have copied the
5153  // vertex information in the same order in which they
5154  // appear in the face; however, this means that we
5155  // impart a coordinate system that is right-handed when
5156  // looked at *from the outside* of the cell if the
5157  // current face has index 0, 2, 4 within a 3d cell, but
5158  // right-handed when looked at *from the inside* for the
5159  // other faces. we fix this by flipping opposite
5160  // vertices if we are on a face 1, 3, 5
5161  if (dim == 3)
5162  if (i % 2 == 1)
5163  std::swap(c_data.vertices[1], c_data.vertices[2]);
5164 
5165  // in 3d, we also need to make sure we copy the manifold
5166  // indicators from the edges of the volume mesh to the
5167  // edges of the surface mesh
5168  //
5169  // one might think that we we can also prescribe
5170  // boundary indicators for edges, but this is only
5171  // possible for edges that aren't just on the boundary
5172  // of the domain (all of the edges we consider are!) but
5173  // that would actually end up at the boundary of the
5174  // surface mesh. there is no easy way to check this, so
5175  // we simply don't do it and instead set it to an
5176  // invalid value that makes sure
5177  // Triangulation::create_triangulation doesn't copy it
5178  if (dim == 3)
5179  for (unsigned int e = 0; e < 4; ++e)
5180  {
5181  // see if we already saw this edge from a
5182  // neighboring face, either in this or the reverse
5183  // orientation. if so, skip it.
5184  {
5185  bool edge_found = false;
5186  for (unsigned int i = 0;
5187  i < subcell_data.boundary_lines.size();
5188  ++i)
5189  if (((subcell_data.boundary_lines[i].vertices[0] ==
5190  map_vert_index[face->line(e)->vertex_index(0)]) &&
5191  (subcell_data.boundary_lines[i].vertices[1] ==
5192  map_vert_index[face->line(e)->vertex_index(
5193  1)])) ||
5194  ((subcell_data.boundary_lines[i].vertices[0] ==
5195  map_vert_index[face->line(e)->vertex_index(1)]) &&
5196  (subcell_data.boundary_lines[i].vertices[1] ==
5197  map_vert_index[face->line(e)->vertex_index(0)])))
5198  {
5199  edge_found = true;
5200  break;
5201  }
5202  if (edge_found == true)
5203  continue; // try next edge of current face
5204  }
5205 
5206  CellData<1> edge;
5207  edge.vertices[0] =
5208  map_vert_index[face->line(e)->vertex_index(0)];
5209  edge.vertices[1] =
5210  map_vert_index[face->line(e)->vertex_index(1)];
5211  edge.boundary_id = numbers::internal_face_boundary_id;
5212  edge.manifold_id = face->line(e)->manifold_id();
5213 
5214  subcell_data.boundary_lines.push_back(edge);
5215  }
5216 
5217 
5218  cells.push_back(c_data);
5219  mapping.push_back(face);
5220  }
5221  }
5222 
5223  // create level 0 surface triangulation
5224  Assert(cells.size() > 0, ExcMessage("No boundary faces selected"));
5225  const_cast<Triangulation<dim - 1, spacedim> &>(
5226  surface_mesh.get_triangulation())
5227  .create_triangulation(vertices, cells, subcell_data);
5228 
5229  // Make the actual mapping
5230  for (typename MeshType<dim - 1, spacedim>::active_cell_iterator cell =
5231  surface_mesh.begin(0);
5232  cell != surface_mesh.end(0);
5233  ++cell)
5234  surface_to_volume_mapping[cell] = mapping.at(cell->index());
5235 
5236  do
5237  {
5238  bool changed = false;
5239 
5240  for (typename MeshType<dim - 1, spacedim>::active_cell_iterator cell =
5241  surface_mesh.begin_active();
5242  cell != surface_mesh.end();
5243  ++cell)
5244  if (surface_to_volume_mapping[cell]->has_children() == true)
5245  {
5246  cell->set_refine_flag();
5247  changed = true;
5248  }
5249 
5250  if (changed)
5251  {
5252  const_cast<Triangulation<dim - 1, spacedim> &>(
5253  surface_mesh.get_triangulation())
5254  .execute_coarsening_and_refinement();
5255 
5256  for (typename MeshType<dim - 1, spacedim>::cell_iterator
5257  surface_cell = surface_mesh.begin();
5258  surface_cell != surface_mesh.end();
5259  ++surface_cell)
5260  for (unsigned int c = 0; c < surface_cell->n_children(); c++)
5261  if (surface_to_volume_mapping.find(surface_cell->child(c)) ==
5262  surface_to_volume_mapping.end())
5263  surface_to_volume_mapping[surface_cell->child(c)] =
5264  surface_to_volume_mapping[surface_cell]->child(c);
5265  }
5266  else
5267  break;
5268  }
5269  while (true);
5270 
5271  return surface_to_volume_mapping;
5272  }
5273 
5274 } // namespace GridGenerator
5275 
5276 // explicit instantiations
5277 #include "grid_generator.inst"
5278 
5279 DEAL_II_NAMESPACE_CLOSE
SubCellData::boundary_lines
std::vector< CellData< 1 > > boundary_lines
Definition: tria.h:258
Triangulation::copy_triangulation
virtual void copy_triangulation(const Triangulation< dim, spacedim > &other_tria)
Definition: tria.cc:10473
GridGenerator::create_union_triangulation
void create_union_triangulation(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result)
Definition: grid_generator.cc:4352
GridReordering::reorder_cells
static void reorder_cells(std::vector< CellData< dim >> &original_cells, const bool use_new_style_ordering=false)
Definition: grid_reordering.cc:1080
numbers::flat_manifold_id
const types::manifold_id flat_manifold_id
Definition: types.h:243
types::manifold_id
unsigned int manifold_id
Definition: types.h:123
FlatManifold::normal_vector
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
Definition: manifold.cc:853
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1366
determinant
Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:2676
Triangulation::get_vertices
const std::vector< Point< spacedim > > & get_vertices() const
Triangulation::create_triangulation_compatibility
virtual void create_triangulation_compatibility(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
Definition: tria.cc:10536
Triangulation::end
cell_iterator end() const
Definition: tria.cc:12004
GeometryInfo::face_to_cell_vertices
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
GridTools::have_same_coarse_mesh
bool have_same_coarse_mesh(const Triangulation< dim, spacedim > &mesh_1, const Triangulation< dim, spacedim > &mesh_2)
Definition: grid_tools_dof_handlers.cc:1220
GridGenerator::simplex
void simplex(Triangulation< dim, dim > &tria, const std::vector< Point< dim >> &vertices)
Triangulation of a d-simplex with (d+1) vertices and mesh cells.
Triangulation::n_active_lines
unsigned int n_active_lines() const
Definition: tria.cc:12743
Triangulation::begin_active_face
active_face_iterator begin_active_face() const
Definition: tria.cc:12117
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Triangulation::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12060
GridGenerator::concentric_hyper_shells
void concentric_hyper_shells(Triangulation< dim > &triangulation, const Point< dim > &center, const double inner_radius=0.125, const double outer_radius=0.25, const unsigned int n_shells=1, const double skewness=0.1, const unsigned int n_cells_per_shell=0, const bool colorize=false)
Definition: grid_generator.cc:4738
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:133
GridGenerator::truncated_cone
void truncated_cone(Triangulation< dim > &tria, const double radius_0=1.0, const double radius_1=0.5, const double half_length=1.0)
GridGenerator::hyper_rectangle
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
Definition: grid_generator.cc:356
GridGenerator::moebius
void moebius(Triangulation< 3, 3 > &tria, const unsigned int n_cells, const unsigned int n_rotations, const double R, const double r)
GridGenerator::flatten_triangulation
void flatten_triangulation(const Triangulation< dim, spacedim1 > &in_tria, Triangulation< dim, spacedim2 > &out_tria)
Definition: grid_generator.cc:4979
GridGenerator::quarter_hyper_shell
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > &center, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
GridGenerator::extrude_triangulation
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result)
Definition: grid_generator.cc:4465
GridGenerator::enclosed_hyper_cube
void enclosed_hyper_cube(Triangulation< dim > &tria, const double left=0., const double right=1., const double thickness=1., const bool colorize=false)
Triangulation::begin_face
face_iterator begin_face() const
Definition: tria.cc:12096
GridGenerator::subdivided_parallelepiped
void subdivided_parallelepiped(Triangulation< dim > &tria, const unsigned int n_subdivisions, const Point< dim >(&corners)[dim], const bool colorize=false)
Definition: grid_generator.cc:849
GridGenerator::ExcInvalidRepetitionsDimension
static::ExceptionBase & ExcInvalidRepetitionsDimension(int arg1)
types::material_id
unsigned int material_id
Definition: types.h:134
GridGenerator::parallelogram
void parallelogram(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
Triangulation::last
cell_iterator last() const
Definition: tria.cc:11955
Tensor::norm
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1301
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1329
CellData::boundary_id
types::boundary_id boundary_id
Definition: tria.h:171
GridGenerator::hyper_shell
void hyper_shell(Triangulation< dim > &tria, const Point< dim > &center, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, bool colorize=false)
StandardExceptions::ExcIndexRange
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
GridGenerator::hyper_sphere
void hyper_sphere(Triangulation< spacedim-1, spacedim > &tria, const Point< spacedim > &center=Point< spacedim >(), const double radius=1.)
Definition: grid_generator.cc:3417
GridGenerator::ExcInvalidInputOrientation
static::ExceptionBase & ExcInvalidInputOrientation()
Triangulation::begin
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:11915
Triangulation::n_active_cells
unsigned int n_active_cells() const
Definition: tria.cc:12509
GridGenerator::cylinder_shell
void cylinder_shell(Triangulation< dim > &tria, const double length, const double inner_radius, const double outer_radius, const unsigned int n_radial_cells=0, const unsigned int n_axial_cells=0)
TransfiniteInterpolationManifold::initialize
void initialize(const Triangulation< dim, spacedim > &triangulation)
Definition: manifold_lib.cc:1439
Triangulation::n_levels
unsigned int n_levels() const
GridGenerator::hyper_cross
void hyper_cross(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &sizes, const bool colorize_cells=false)
A center cell with stacks of cell protruding from each surface.
Definition: grid_generator.cc:2264
Triangulation::set_manifold
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
Definition: tria.cc:10268
Triangulation::cell_iterators
IteratorRange< cell_iterator > cell_iterators() const
Definition: tria.cc:12051
GridGenerator::hyper_L
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
GridTools::delete_unused_vertices
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:434
numbers::PI
static const double PI
Definition: numbers.h:143
Triangulation::execute_coarsening_and_refinement
virtual void execute_coarsening_and_refinement()
Definition: tria.cc:13229
GridGenerator::half_hyper_ball
void half_hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1.)
StandardExceptions::ExcMessage
static::ExceptionBase & ExcMessage(std::string arg1)
Triangulation::n_cells
unsigned int n_cells() const
Definition: tria.cc:12501
TableBase::size
size_type size(const unsigned int i) const
GridGenerator::cylinder
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
Triangulation::create_triangulation
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
Definition: tria.cc:10556
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1227
Triangulation::set_all_manifold_ids
void set_all_manifold_ids(const types::manifold_id number)
Definition: tria.cc:10310
numbers::invalid_boundary_id
const types::boundary_id invalid_boundary_id
Definition: types.h:207
GridGenerator::quarter_hyper_ball
void quarter_hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1.)
CellData::material_id
types::material_id material_id
Definition: tria.h:160
GridGenerator::subdivided_hyper_cube
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1.)
Definition: grid_generator.cc:1139
GridGenerator::create_triangulation_with_removed_cells
void create_triangulation_with_removed_cells(const Triangulation< dim, spacedim > &input_triangulation, const std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_remove, Triangulation< dim, spacedim > &result)
Definition: grid_generator.cc:4413
InterGridMap::make_mapping
void make_mapping(const MeshType &source_grid, const MeshType &destination_grid)
Definition: intergrid_map.cc:46
Triangulation::n_lines
unsigned int n_lines() const
Definition: tria.cc:12651
GridGenerator::extract_boundary_mesh
std::map< typename MeshType< dim-1, spacedim >::cell_iterator, typename MeshType< dim, spacedim >::face_iterator > extract_boundary_mesh(const MeshType< dim, spacedim > &volume_mesh, MeshType< dim-1, spacedim > &surface_mesh, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())
Definition: grid_generator.cc:5083
GridTools::rotate
void rotate(const double angle, Triangulation< 2 > &triangulation)
Definition: grid_tools.cc:707
GridGenerator::hyper_cube_slit
void hyper_cube_slit(Triangulation< dim > &tria, const double left=0., const double right=1., const bool colorize=false)
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:96
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
StandardExceptions::ExcLowerRange
static::ExceptionBase & ExcLowerRange(int arg1, int arg2)
GridGenerator::subdivided_hyper_rectangle
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
Definition: grid_generator.cc:1163
GridGenerator::merge_triangulations
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false)
Definition: grid_generator.cc:4325
CellData::manifold_id
types::manifold_id manifold_id
Definition: tria.h:182
GridGenerator::ExcInvalidRadii
static::ExceptionBase & ExcInvalidRadii()
swap
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1353
GridGenerator::cheese
void cheese(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &holes)
Rectangular domain with rectangular pattern of holes.
Definition: grid_generator.cc:1796
Triangulation::n_active_faces
unsigned int n_active_faces() const
Definition: tria.cc:12561
Triangulation::end_face
face_iterator end_face() const
Definition: tria.cc:12138
GridGenerator::parallelepiped
void parallelepiped(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
Definition: grid_generator.cc:835
GridGenerator::general_cell
void general_cell(Triangulation< dim > &tria, const std::vector< Point< dim >> &vertices, const bool colorize=false)
Definition: grid_generator.cc:775
GridGenerator::half_hyper_shell
void half_hyper_shell(Triangulation< dim > &tria, const Point< dim > &center, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
SubCellData::boundary_quads
std::vector< CellData< 2 > > boundary_quads
Definition: tria.h:266
GridGenerator::hyper_cube
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
Definition: grid_generator.cc:420
Triangulation::n_quads
unsigned int n_quads() const
Definition: tria.cc:12909
Triangulation::refine_global
void refine_global(const unsigned int times=1)
Definition: tria.cc:10773
GridGenerator::hyper_ball
void hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
StandardExceptions::ExcNotImplemented
static::ExceptionBase & ExcNotImplemented()
Table
Definition: table.h:37
GridReordering::invert_all_cells_of_negative_grid
static void invert_all_cells_of_negative_grid(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &original_cells)
Triangulation::set_all_manifold_ids_on_boundary
void set_all_manifold_ids_on_boundary(const types::manifold_id number)
Definition: tria.cc:10329
numbers::internal_face_boundary_id
const types::boundary_id internal_face_boundary_id
Definition: types.h:223
GridGenerator::plate_with_a_hole
void plate_with_a_hole(Triangulation< dim > &tria, const double inner_radius=0.4, const double outer_radius=1., const double pad_bottom=2., const double pad_top=2., const double pad_left=1., const double pad_right=1., const Point< dim > center=Point< dim >(), const types::manifold_id polar_manifold_id=0, const types::manifold_id tfi_manifold_id=1, const double L=1., const unsigned int n_slices=2, const bool colorize=false)
Rectangular plate with an (offset) cylindrical hole.
Triangulation::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11935
GridGenerator::hyper_cube_with_cylindrical_hole
void hyper_cube_with_cylindrical_hole(Triangulation< dim > &triangulation, const double inner_radius=.25, const double outer_radius=.5, const double L=.5, const unsigned int repetitions=1, const bool colorize=false)
numbers::invalid_material_id
const types::material_id invalid_material_id
Definition: types.h:196
Triangulation::n_vertices
unsigned int n_vertices() const
types::boundary_id
unsigned int boundary_id
Definition: types.h:111
GridTools::shift
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:698
Triangulation::clear
virtual void clear()
Definition: tria.cc:10232
CellData::vertices
unsigned int vertices[GeometryInfo< structdim >::vertices_per_cell]
Definition: tria.h:141
GridGenerator::ExcInvalidRepetitions
static::ExceptionBase & ExcInvalidRepetitions(int arg1)
StandardExceptions::ExcInternalError
static::ExceptionBase & ExcInternalError()
GridTools::delete_duplicated_vertices
void delete_duplicated_vertices(std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata, std::vector< unsigned int > &considered_vertices, const double tol=1e-12)
Definition: grid_tools.cc:533
Generated by   doxygen 1.8.11