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