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