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