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