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