Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_generator.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17
21
28#include <deal.II/grid/tria.h>
31
33
34#include <array>
35#include <cmath>
36#include <limits>
37
38
40
41// work around the problem that doxygen for some reason lists all template
42// specializations in this file
43#ifndef DOXYGEN
44
45namespace GridGenerator
46{
47 namespace Airfoil
48 {
50 // airfoil configuration
51 : airfoil_type("NACA")
52 , naca_id("2412")
53 , joukowski_center(-0.1, 0.14)
54 , airfoil_length(1.0)
55 // far field
56 , height(30.0)
57 , length_b2(15.0)
58 // mesh
59 , incline_factor(0.35)
60 , bias_factor(2.5)
61 , refinements(2)
62 , n_subdivision_x_0(3)
63 , n_subdivision_x_1(2)
64 , n_subdivision_x_2(5)
65 , n_subdivision_y(3)
66 , airfoil_sampling_factor(2)
67 {
68 Assert(
69 airfoil_length <= height,
71 "Mesh is to small to enclose airfoil! Choose larger field or smaller"
72 " chord length!"));
73 Assert(incline_factor < 1.0 && incline_factor >= 0.0,
74 ExcMessage("incline_factor has to be in [0,1)!"));
75 }
76
77
78
79 void
80 AdditionalData::add_parameters(ParameterHandler &prm)
81 {
82 prm.enter_subsection("FarField");
83 {
84 prm.add_parameter(
85 "Height",
86 height,
87 "Mesh height measured from airfoil nose to horizontal boundaries");
88 prm.add_parameter(
89 "LengthB2",
90 length_b2,
91 "Length measured from airfoil leading edge to vertical outlet boundary");
92 prm.add_parameter(
93 "InclineFactor",
94 incline_factor,
95 "Define obliqueness of the vertical mesh around the airfoil");
96 }
97 prm.leave_subsection();
98
99 prm.enter_subsection("AirfoilType");
100 {
101 prm.add_parameter(
102 "Type",
103 airfoil_type,
104 "Type of airfoil geometry, either NACA or Joukowski airfoil",
105 Patterns::Selection("NACA|Joukowski"));
106 }
107 prm.leave_subsection();
108
109 prm.enter_subsection("NACA");
110 {
111 prm.add_parameter("NacaId", naca_id, "Naca serial number");
112 }
113 prm.leave_subsection();
114
115 prm.enter_subsection("Joukowski");
116 {
117 prm.add_parameter("Center",
118 joukowski_center,
119 "Joukowski circle center coordinates");
120 prm.add_parameter("AirfoilLength",
121 airfoil_length,
122 "Joukowski airfoil length leading to trailing edge");
123 }
124 prm.leave_subsection();
125
126 prm.enter_subsection("Mesh");
127 {
128 prm.add_parameter("Refinements",
129 refinements,
130 "Number of global refinements");
131 prm.add_parameter(
132 "NumberSubdivisionX0",
133 n_subdivision_x_0,
134 "Number of subdivisions along the airfoil in blocks with material ID 1 and 4");
135 prm.add_parameter(
136 "NumberSubdivisionX1",
137 n_subdivision_x_1,
138 "Number of subdivisions along the airfoil in blocks with material ID 2 and 5");
139 prm.add_parameter(
140 "NumberSubdivisionX2",
141 n_subdivision_x_2,
142 "Number of subdivisions in horizontal direction on the right of the trailing edge, i.e., blocks with material ID 3 and 6");
143 prm.add_parameter("NumberSubdivisionY",
144 n_subdivision_y,
145 "Number of subdivisions normal to airfoil");
146 prm.add_parameter(
147 "BiasFactor",
148 bias_factor,
149 "Factor to obtain a finer mesh at the airfoil surface");
150 }
151 prm.leave_subsection();
152 }
153
154
155 namespace
156 {
160 class MeshGenerator
161 {
162 public:
163 // IDs of the mesh blocks
164 static const unsigned int id_block_1 = 1;
165 static const unsigned int id_block_2 = 2;
166 static const unsigned int id_block_3 = 3;
167 static const unsigned int id_block_4 = 4;
168 static const unsigned int id_block_5 = 5;
169 static const unsigned int id_block_6 = 6;
170
174 MeshGenerator(const AdditionalData &data)
175 : refinements(data.refinements)
176 , n_subdivision_x_0(data.n_subdivision_x_0)
177 , n_subdivision_x_1(data.n_subdivision_x_1)
178 , n_subdivision_x_2(data.n_subdivision_x_2)
179 , n_subdivision_y(data.n_subdivision_y)
180 , height(data.height)
181 , length_b2(data.length_b2)
182 , incline_factor(data.incline_factor)
183 , bias_factor(data.bias_factor)
184 , edge_length(1.0)
185 , n_cells_x_0(Utilities::pow(2, refinements) * n_subdivision_x_0)
186 , n_cells_x_1(Utilities::pow(2, refinements) * n_subdivision_x_1)
187 , n_cells_x_2(Utilities::pow(2, refinements) * n_subdivision_x_2)
188 , n_cells_y(Utilities::pow(2, refinements) * n_subdivision_y)
189 , n_points_on_each_side(n_cells_x_0 + n_cells_x_1 + 1)
190 // create points on the airfoil
191 , airfoil_1D(set_airfoil_length(
192 // call either the 'joukowski' or 'naca' static member function
193 data.airfoil_type == "Joukowski" ?
194 joukowski(data.joukowski_center,
195 n_points_on_each_side,
196 data.airfoil_sampling_factor) :
197 (data.airfoil_type == "NACA" ?
198 naca(data.naca_id,
199 n_points_on_each_side,
200 data.airfoil_sampling_factor) :
201 std::array<std::vector<Point<2>>, 2>{
202 {std::vector<Point<2>>{Point<2>(0), Point<2>(1)},
203 std::vector<Point<2>>{
204 Point<2>(0),
205 Point<2>(
206 1)}}} /* dummy vector since we are asserting later*/),
207 data.airfoil_length))
208 , end_b0_x_u(airfoil_1D[0][n_cells_x_0](0))
209 , end_b0_x_l(airfoil_1D[1][n_cells_x_0](0))
210 , nose_x(airfoil_1D[0].front()(0))
211 , tail_x(airfoil_1D[0].back()(0))
212 , tail_y(airfoil_1D[0].back()(1))
213 , center_mesh(0.5 * std::abs(end_b0_x_u + end_b0_x_l))
214 , length_b1_x(tail_x - center_mesh)
215 , gamma(std::atan(height /
216 (edge_length + std::abs(nose_x - center_mesh))))
217 // points on coarse grid
218 // coarse grid has to be symmetric in respect to x-axis to allow
219 // periodic BC and make sure that interpolate() works
220 , A(nose_x - edge_length, 0)
221 , B(nose_x, 0)
222 , C(center_mesh, +std::abs(nose_x - center_mesh) * std::tan(gamma))
223 , D(center_mesh, height)
224 , E(center_mesh, -std::abs(nose_x - center_mesh) * std::tan(gamma))
225 , F(center_mesh, -height)
226 , G(tail_x, height)
227 , H(tail_x, 0)
228 , I(tail_x, -height)
229 , J(tail_x + length_b2, 0)
230 , K(J(0), G(1))
231 , L(J(0), I(1))
232 {
233 Assert(data.airfoil_type == "Joukowski" ||
234 data.airfoil_type == "NACA",
235 ExcMessage("Unknown airfoil type."));
236 }
237
241 void
243 Triangulation<2> & tria_grid,
244 std::vector<GridTools::PeriodicFacePair<
245 typename Triangulation<2>::cell_iterator>> *periodic_faces) const
246 {
247 make_coarse_grid(tria_grid);
248
249 set_boundary_ids(tria_grid);
250
251 if (periodic_faces != nullptr)
252 {
254 tria_grid, 5, 4, 1, *periodic_faces);
255 tria_grid.add_periodicity(*periodic_faces);
256 }
257
258 tria_grid.refine_global(refinements);
259 interpolate(tria_grid);
260 }
261
265 void
268 std::vector<GridTools::PeriodicFacePair<
269 typename Triangulation<2>::cell_iterator>> *periodic_faces) const
270 {
271 (void)parallel_grid;
272 (void)periodic_faces;
273
274 AssertThrow(false, ExcMessage("Not implemented, yet!")); // TODO [PM]
275 }
276
277 private:
278 // number of global refinements
279 const unsigned int refinements;
280
281 // number of subdivisions of coarse grid in blocks 1 and 4
282 const unsigned int n_subdivision_x_0;
283
284 // number of subdivisions of coarse grid in blocks 2 and 5
285 const unsigned int n_subdivision_x_1;
286
287 // number of subdivisions of coarse grid in blocks 3 and 6
288 const unsigned int n_subdivision_x_2;
289
290 // number of subdivisions of coarse grid in all blocks (normal to
291 // airfoil or in y-direction, respectively)
292 const unsigned int n_subdivision_y;
293
294 // height of mesh, i.e. length JK or JL and radius of semicircle
295 // (C-Mesh) that arises after interpolation in blocks 1 and 4
296 const double height;
297
298 // length block 3 and 6
299 const double length_b2;
300
301 // factor to move points G and I horizontal to the right, i.e. make
302 // faces HG and HI inclined instead of vertical
303 const double incline_factor;
304
305 // bias factor (if factor goes to zero than equal y = x)
306 const double bias_factor;
307
308 // x-distance between coarse grid vertices A and B, i.e. used only once;
309 const double edge_length;
310
311 // number of cells (after refining) in block 1 and 4 along airfoil
312 const unsigned int n_cells_x_0;
313
314 // number of cells (after refining) in block 2 and 5 along airfoil
315 const unsigned int n_cells_x_1;
316
317 // number of cells (after refining) in block 3 and 6 in x-direction
318 const unsigned int n_cells_x_2;
319
320 // number of cells (after refining) in all blocks normal to airfoil or
321 // in y-direction, respectively
322 const unsigned int n_cells_y;
323
324 // number of airfoil points on each side
325 const unsigned int n_points_on_each_side;
326
327 // vector containing upper/lower airfoil points. First and last point
328 // are identical
329 const std::array<std::vector<Point<2>>, 2> airfoil_1D;
330
331 // x-coordinate of n-th airfoilpoint where n indicates number of cells
332 // in block 1. end_b0_x_u = end_b0_x_l for symmetric airfoils
333 const double end_b0_x_u;
334
335 // x-coordinate of n-th airfoilpoint where n indicates number of cells
336 // in block 4. end_b0_x_u = end_b0_x_l for symmetric airfoils
337 const double end_b0_x_l;
338
339 // x-coordinate of first airfoil point in airfoil_1D[0] and
340 // airfoil_1D[1]
341 const double nose_x;
342
343 // x-coordinate of last airfoil point in airfoil_1D[0] and airfoil_1D[1]
344 const double tail_x;
345
346 // y-coordinate of last airfoil point in airfoil_1D[0] and airfoil_1D[1]
347 const double tail_y;
348
349 // x-coordinate of C,D,E,F indicating ending of blocks 1 and 4 or
350 // beginning of blocks 2 and 5, respectively
351 const double center_mesh;
352
353 // length of blocks 2 and 5
354 const double length_b1_x;
355
356 // angle enclosed between faces DAB and FAB
357 const double gamma;
358
359
360
381 const Point<2> A, B, C, D, E, F, G, H, I, J, K, L;
382
383
384
420 static std::array<std::vector<Point<2>>, 2>
421 joukowski(const Point<2> & centerpoint,
422 const unsigned int number_points,
423 const unsigned int factor)
424 {
425 std::array<std::vector<Point<2>>, 2> airfoil_1D;
426 const unsigned int total_points = 2 * number_points - 2;
427 const unsigned int n_airfoilpoints = factor * total_points;
428 // joukowski points on the entire airfoil, i.e. upper and lower side
429 const auto jouk_points =
430 joukowski_transform(joukowski_circle(centerpoint, n_airfoilpoints));
431
432 // vectors to collect airfoil points on either upper or lower side
433 std::vector<Point<2>> upper_points;
434 std::vector<Point<2>> lower_points;
435
436 {
437 // find point on nose and point on tail
438 unsigned int nose_index = 0;
439 unsigned int tail_index = 0;
440 double nose_x_coordinate = 0;
441 double tail_x_coordinate = 0;
442
443
444 // find index in vector to nose point (min) and tail point (max)
445 for (unsigned int i = 0; i < jouk_points.size(); ++i)
446 {
447 if (jouk_points[i](0) < nose_x_coordinate)
448 {
449 nose_x_coordinate = jouk_points[i](0);
450 nose_index = i;
451 }
452 if (jouk_points[i](0) > tail_x_coordinate)
453 {
454 tail_x_coordinate = jouk_points[i](0);
455 tail_index = i;
456 }
457 }
458
459 // copy point on upper side of airfoil
460 for (unsigned int i = tail_index; i < jouk_points.size(); ++i)
461 upper_points.emplace_back(jouk_points[i]);
462 for (unsigned int i = 0; i <= nose_index; ++i)
463 upper_points.emplace_back(jouk_points[i]);
464 std::reverse(upper_points.begin(), upper_points.end());
465
466 // copy point on lower side of airfoil
467 lower_points.insert(lower_points.end(),
468 jouk_points.begin() + nose_index,
469 jouk_points.begin() + tail_index + 1);
470 }
471
472 airfoil_1D[0] = make_points_equidistant(upper_points, number_points);
473 airfoil_1D[1] = make_points_equidistant(lower_points, number_points);
474
475 // move nose to origin
476 auto move_nose_to_origin = [](std::vector<Point<2>> &vector) {
477 const double nose_x_pos = vector.front()(0);
478 for (auto &i : vector)
479 i(0) -= nose_x_pos;
480 };
481
482 move_nose_to_origin(airfoil_1D[1]);
483 move_nose_to_origin(airfoil_1D[0]);
484
485 return airfoil_1D;
486 }
487
512 static std::vector<Point<2>>
513 joukowski_circle(const Point<2> & center,
514 const unsigned int number_points)
515 {
516 std::vector<Point<2>> circle_points;
517
518 // Create Circle with number_points - points
519 // unsigned int number_points = 2 * points_per_side - 2;
520
521 // Calculate radius so that point (x=1|y=0) is enclosed - requirement
522 // for Joukowski transform
523 const double radius = std::sqrt(center(1) * center(1) +
524 (1 - center(0)) * (1 - center(0)));
525 const double radius_test = std::sqrt(
526 center(1) * center(1) + (1 + center(0)) * (1 + center(0)));
527 // Make sure point (x=-1|y=0) is enclosed by the circle
528 (void)radius_test;
530 radius_test < radius,
532 "Error creating lower circle: Circle for Joukowski-transform does"
533 " not enclose point zeta = -1! Choose different center "
534 "coordinate."));
535 // Create a full circle with radius 'radius' around Point 'center' of
536 // (number_points) equidistant points.
537 const double theta = 2 * numbers::PI / number_points;
538 // first point is leading edge then counterclockwise
539 for (unsigned int i = 0; i < number_points; ++i)
540 circle_points.emplace_back(center[0] - radius * cos(i * theta),
541 center[1] - radius * sin(i * theta));
542
543 return circle_points;
544 }
545
554 static std::vector<Point<2>>
555 joukowski_transform(const std::vector<Point<2>> &circle_points)
556 {
557 std::vector<Point<2>> joukowski_points(circle_points.size());
558
559 // transform each point
560 for (unsigned int i = 0; i < circle_points.size(); ++i)
561 {
562 const double chi = circle_points[i](0);
563 const double eta = circle_points[i](1);
564 const std::complex<double> zeta(chi, eta);
565 const std::complex<double> z = zeta + 1. / zeta;
566
567 joukowski_points[i] = {real(z), imag(z)};
568 }
569 return joukowski_points;
570 }
571
588 static std::array<std::vector<Point<2>>, 2>
589 naca(const std::string &serialnumber,
590 const unsigned int number_points,
591 const unsigned int factor)
592 {
593 // number of non_equidistant airfoilpoints among which will be
594 // interpolated
595 const unsigned int n_airfoilpoints = factor * number_points;
596
597 // create equidistant airfoil points for upper and lower side
598 return {{make_points_equidistant(
599 naca_create_points(serialnumber, n_airfoilpoints, true),
600 number_points),
601 make_points_equidistant(
602 naca_create_points(serialnumber, n_airfoilpoints, false),
603 number_points)}};
604 }
605
617 static std::vector<Point<2>>
618 naca_create_points(const std::string &serialnumber,
619 const unsigned int number_points,
620 const bool is_upper)
621 {
622 Assert(serialnumber.size() == 4,
623 ExcMessage("This NACA-serial number is not implemented!"));
624
625 return naca_create_points_4_digits(serialnumber,
626 number_points,
627 is_upper);
628 }
629
644 static std::vector<Point<2>>
645 naca_create_points_4_digits(const std::string &serialnumber,
646 const unsigned int number_points,
647 const bool is_upper)
648 {
649 // conversion string (char * ) to int
650 const unsigned int digit_0 = (serialnumber[0] - '0');
651 const unsigned int digit_1 = (serialnumber[1] - '0');
652 const unsigned int digit_2 = (serialnumber[2] - '0');
653 const unsigned int digit_3 = (serialnumber[3] - '0');
654
655 const unsigned int digit_23 = 10 * digit_2 + digit_3;
656
657 // maximum thickness in percentage of the cord
658 const double t = static_cast<double>(digit_23) / 100.0;
659
660 std::vector<Point<2>> naca_points;
661
662 if (digit_0 == 0 && digit_1 == 0) // is symmetric
663 for (unsigned int i = 0; i < number_points; ++i)
664 {
665 const double x = i * 1 / (1.0 * number_points - 1);
666 const double y_t =
667 5 * t *
668 (0.2969 * std::pow(x, 0.5) - 0.126 * x -
669 0.3516 * std::pow(x, 2) + 0.2843 * std::pow(x, 3) -
670 0.1036 * std::pow(x, 4)); // half thickness at a position x
671
672 if (is_upper)
673 naca_points.emplace_back(x, +y_t);
674 else
675 naca_points.emplace_back(x, -y_t);
676 }
677 else // is asymmetric
678 for (unsigned int i = 0; i < number_points; ++i)
679 {
680 const double m = 1.0 * digit_0 / 100; // max. chamber
681 const double p = 1.0 * digit_1 / 10; // location of max. chamber
682 const double x = i * 1 / (1.0 * number_points - 1);
683
684 const double y_c =
685 (x <= p) ? m / std::pow(p, 2) * (2 * p * x - std::pow(x, 2)) :
686 m / std::pow(1 - p, 2) *
687 ((1 - 2 * p) + 2 * p * x - std::pow(x, 2));
688
689 const double dy_c = (x <= p) ?
690 2 * m / std::pow(p, 2) * (p - x) :
691 2 * m / std::pow(1 - p, 2) * (p - x);
692
693 const double y_t =
694 5 * t *
695 (0.2969 * std::pow(x, 0.5) - 0.126 * x -
696 0.3516 * std::pow(x, 2) + 0.2843 * std::pow(x, 3) -
697 0.1036 * std::pow(x, 4)); // half thickness at a position x
698
699 const double theta = std::atan(dy_c);
700
701 if (is_upper)
702 naca_points.emplace_back(x - y_t * std::sin(theta),
703 y_c + y_t * std::cos(theta));
704 else
705 naca_points.emplace_back(x + y_t * std::sin(theta),
706 y_c - y_t * std::cos(theta));
707 }
708
709 return naca_points;
710 }
711
712
713
722 static std::array<std::vector<Point<2>>, 2>
723 set_airfoil_length(const std::array<std::vector<Point<2>>, 2> &input,
724 const double desired_len)
725 {
726 std::array<std::vector<Point<2>>, 2> output;
727 output[0] = set_airfoil_length(input[0], desired_len);
728 output[1] = set_airfoil_length(input[1], desired_len);
729
730 return output;
731 }
732
740 static std::vector<Point<2>>
741 set_airfoil_length(const std::vector<Point<2>> &input,
742 const double desired_len)
743 {
744 std::vector<Point<2>> output = input;
745
746 const double scale =
747 desired_len / input.front().distance(input.back());
748
749 for (auto &x : output)
750 x *= scale;
751
752 return output;
753 }
754
765 static std::vector<Point<2>>
766 make_points_equidistant(
767 const std::vector<Point<2>> &non_equidistant_points,
768 const unsigned int number_points)
769 {
770 const unsigned int n_points =
771 non_equidistant_points
772 .size(); // number provided airfoilpoints to interpolate
773
774 // calculate arclength
775 std::vector<double> arclength_L(non_equidistant_points.size(), 0);
776 for (unsigned int i = 0; i < non_equidistant_points.size() - 1; ++i)
777 arclength_L[i + 1] =
778 arclength_L[i] +
779 non_equidistant_points[i + 1].distance(non_equidistant_points[i]);
780
781
782 const auto airfoil_length =
783 arclength_L.back(); // arclength upper or lower side
784 const auto deltaX = airfoil_length / (number_points - 1);
785
786 // Create equidistant points: keep the first (and last) point
787 // unchanged
788 std::vector<Point<2>> equidist(
789 number_points); // number_points is required points on each side for
790 // mesh
791 equidist[0] = non_equidistant_points[0];
792 equidist[number_points - 1] = non_equidistant_points[n_points - 1];
793
794
795 // loop over all subsections
796 for (unsigned int j = 0, i = 1; j < n_points - 1; ++j)
797 {
798 // get reference left and right end of this section
799 const auto Lj = arclength_L[j];
800 const auto Ljp = arclength_L[j + 1];
801
802 while (Lj <= i * deltaX && i * deltaX <= Ljp &&
803 i < number_points - 1)
804 {
805 equidist[i] = Point<2>((i * deltaX - Lj) / (Ljp - Lj) *
806 (non_equidistant_points[j + 1] -
807 non_equidistant_points[j]) +
808 non_equidistant_points[j]);
809 ++i;
810 }
811 }
812 return equidist;
813 }
814
815
816
823 void
824 make_coarse_grid(Triangulation<2> &tria) const
825 {
826 // create vector of serial triangulations for each block and
827 // temporary storage for merging them
828 std::vector<Triangulation<2>> trias(10);
829
830 // helper function to create a subdivided quadrilateral
831 auto make = [](Triangulation<2> & tria,
832 const std::vector<Point<2>> & corner_vertices,
833 const std::vector<unsigned int> &repetitions,
834 const unsigned int material_id) {
835 // create subdivided rectangle with corner points (-1,-1)
836 // and (+1, +1). It serves as reference system
838 repetitions,
839 {-1, -1},
840 {+1, +1});
841
842 // move all vertices to the correct position
843 for (auto it = tria.begin_vertex(); it < tria.end_vertex(); ++it)
844 {
845 auto & point = it->vertex();
846 const double xi = point(0);
847 const double eta = point(1);
848
849 // bilinear mapping
850 point = 0.25 * ((1 - xi) * (1 - eta) * corner_vertices[0] +
851 (1 + xi) * (1 - eta) * corner_vertices[1] +
852 (1 - xi) * (1 + eta) * corner_vertices[2] +
853 (1 + xi) * (1 + eta) * corner_vertices[3]);
854 }
855
856 // set material id of block
857 for (auto cell : tria.active_cell_iterators())
858 cell->set_material_id(material_id);
859 };
860
861 // create a subdivided quadrilateral for each block (see last number
862 // of block id)
863 make(trias[0],
864 {A, B, D, C},
865 {n_subdivision_y, n_subdivision_x_0},
866 id_block_1);
867 make(trias[1],
868 {F, E, A, B},
869 {n_subdivision_y, n_subdivision_x_0},
870 id_block_4);
871 make(trias[2],
872 {C, H, D, G},
873 {n_subdivision_x_1, n_subdivision_y},
874 id_block_2);
875 make(trias[3],
876 {F, I, E, H},
877 {n_subdivision_x_1, n_subdivision_y},
878 id_block_5);
879 make(trias[4],
880 {H, J, G, K},
881 {n_subdivision_x_2, n_subdivision_y},
882 id_block_3);
883 make(trias[5],
884 {I, L, H, J},
885 {n_subdivision_x_2, n_subdivision_y},
886 id_block_6);
887
888
889 // merge triangulation (warning: do not change the order here since
890 // this might change the face ids)
891 GridGenerator::merge_triangulations(trias[0], trias[1], trias[6]);
892 GridGenerator::merge_triangulations(trias[2], trias[3], trias[7]);
893 GridGenerator::merge_triangulations(trias[4], trias[5], trias[8]);
894 GridGenerator::merge_triangulations(trias[6], trias[7], trias[9]);
895 GridGenerator::merge_triangulations(trias[8], trias[9], tria);
896 }
897
898 /*
899 * Loop over all (cells and) boundary faces of a given triangulation
900 * and set the boundary_ids depending on the material_id of the cell and
901 * the face number. The resulting boundary_ids are:
902 * - 0: inlet
903 * - 1: outlet
904 * - 2: upper airfoil surface (aka. suction side)
905 * - 3, lower airfoil surface (aka. pressure side),
906 * - 4: upper far-field side
907 * - 5: lower far-field side
908 */
909 static void
910 set_boundary_ids(Triangulation<2> &tria)
911 {
912 for (auto cell : tria.active_cell_iterators())
913 for (unsigned int f : GeometryInfo<2>::face_indices())
914 {
915 if (cell->face(f)->at_boundary() == false)
916 continue;
917
918 const auto mid = cell->material_id();
919
920 if ((mid == id_block_1 && f == 0) ||
921 (mid == id_block_4 && f == 0))
922 cell->face(f)->set_boundary_id(0); // inlet
923 else if ((mid == id_block_3 && f == 0) ||
924 (mid == id_block_6 && f == 2))
925 cell->face(f)->set_boundary_id(1); // outlet
926 else if ((mid == id_block_1 && f == 1) ||
927 (mid == id_block_2 && f == 1))
928 cell->face(f)->set_boundary_id(2); // upper airfoil side
929 else if ((mid == id_block_4 && f == 1) ||
930 (mid == id_block_5 && f == 3))
931 cell->face(f)->set_boundary_id(3); // lower airfoil side
932 else if ((mid == id_block_2 && f == 0) ||
933 (mid == id_block_3 && f == 2))
934 cell->face(f)->set_boundary_id(4); // upper far-field side
935 else if ((mid == id_block_5 && f == 2) ||
936 (mid == id_block_6 && f == 0))
937 cell->face(f)->set_boundary_id(5); // lower far-field side
938 else
939 Assert(false, ExcIndexRange(mid, id_block_1, id_block_6));
940 }
941 }
942
943 /*
944 * Interpolate all vertices of the given triangulation onto the airfoil
945 * geometry, depending on the material_id of the block.
946 * Due to symmetry of coarse grid in respect to
947 * x-axis (by definition of points A-L), blocks 1&4, 2&4 and 3&6 can be
948 * interpolated with the same geometric computations Consider a
949 * bias_factor and incline_factor during interpolation to obtain a more
950 * dense mesh next to airfoil geometry and receive an inclined boundary
951 * between block 2&3 and 5&6, respectively
952 */
953 void
955 {
956 // array storing the information if a vertex was processed
957 std::vector<bool> vertex_processed(tria.n_vertices(), false);
958
959 // rotation matrix for clockwise rotation of block 1 by angle gamma
960 const Tensor<2, 2, double> rotation_matrix_1 =
962 const Tensor<2, 2, double> rotation_matrix_2 =
963 transpose(rotation_matrix_1);
964
965 // horizontal offset in order to place coarse-grid node A in the
966 // origin
967 const Point<2, double> horizontal_offset(A(0), 0.0);
968
969 // Move block 1 so that face BC coincides the x-axis
970 const Point<2, double> trapeze_offset(0.0,
971 std::sin(gamma) * edge_length);
972
973 // loop over vertices of all cells
974 for (auto &cell : tria)
975 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
976 {
977 // vertex has been already processed: nothing to do
978 if (vertex_processed[cell.vertex_index(v)])
979 continue;
980
981 // mark vertex as processed
982 vertex_processed[cell.vertex_index(v)] = true;
983
984 auto &node = cell.vertex(v);
985
986 // distinguish blocks
987 if (cell.material_id() == id_block_1 ||
988 cell.material_id() == id_block_4) // block 1 and 4
989 {
990 // step 1: rotate block 1 clockwise by gamma and move block
991 // 1 so that A(0) is on y-axis so that faces AD and BC are
992 // horizontal. This simplifies the computation of the
993 // required indices for interpolation (all x-nodes are
994 // positive) Move trapeze to be in first quadrant by adding
995 // trapeze_offset
996 Point<2, double> node_;
997 if (cell.material_id() == id_block_1)
998 {
999 node_ = Point<2, double>(rotation_matrix_1 *
1000 (node - horizontal_offset) +
1001 trapeze_offset);
1002 }
1003 // step 1: rotate block 4 counterclockwise and move down so
1004 // that trapeze is located in fourth quadrant (subtracting
1005 // trapeze_offset)
1006 else if (cell.material_id() == id_block_4)
1007 {
1008 node_ = Point<2, double>(rotation_matrix_2 *
1009 (node - horizontal_offset) -
1010 trapeze_offset);
1011 }
1012 // step 2: compute indices ix and iy and interpolate
1013 // trapezoid to a rectangle of length pi/2.
1014 {
1015 const double trapeze_height =
1016 std::sin(gamma) * edge_length;
1017 const double L = height / std::sin(gamma);
1018 const double l_a = std::cos(gamma) * edge_length;
1019 const double l_b = trapeze_height * std::tan(gamma);
1020 const double x1 = std::abs(node_(1)) / std::tan(gamma);
1021 const double x2 = L - l_a - l_b;
1022 const double x3 = std::abs(node_(1)) * std::tan(gamma);
1023 const double Dx = x1 + x2 + x3;
1024 const double deltax =
1025 (trapeze_height - std::abs(node_(1))) / std::tan(gamma);
1026 const double dx = Dx / n_cells_x_0;
1027 const double dy = trapeze_height / n_cells_y;
1028 const int ix =
1029 static_cast<int>(std::round((node_(0) - deltax) / dx));
1030 const int iy =
1031 static_cast<int>(std::round(std::abs(node_(1)) / dy));
1032
1033 node_(0) = numbers::PI / 2 * (1.0 * ix) / n_cells_x_0;
1034 node_(1) = height * (1.0 * iy) / n_cells_y;
1035 }
1036
1037 // step 3: Interpolation between semicircle (of C-Mesh) and
1038 // airfoil contour
1039 {
1040 const double dx = numbers::PI / 2 / n_cells_x_0;
1041 const double dy = height / n_cells_y;
1042 const int ix =
1043 static_cast<int>(std::round(node_(0) / dx));
1044 const int iy =
1045 static_cast<int>(std::round(node_(1) / dy));
1046 const double alpha =
1047 bias_alpha(1 - (1.0 * iy) / n_cells_y);
1048 const double theta = node_(0);
1049 const Point<2> p(-height * std::cos(theta) + center_mesh,
1050 ((cell.material_id() == id_block_1) ?
1051 (height) :
1052 (-height)) *
1053 std::sin(theta));
1054 node =
1055 airfoil_1D[(
1056 (cell.material_id() == id_block_1) ? (0) : (1))][ix] *
1057 alpha +
1058 p * (1 - alpha);
1059 }
1060 }
1061 else if (cell.material_id() == id_block_2 ||
1062 cell.material_id() == id_block_5) // block 2 and 5
1063 {
1064 // geometric parameters and indices for interpolation
1065 Assert(
1066 (std::abs(D(1) - C(1)) == std::abs(F(1) - E(1))) &&
1067 (std::abs(C(1)) == std::abs(E(1))) &&
1068 (std::abs(G(1)) == std::abs(I(1))),
1069 ExcMessage(
1070 "Points D,C,G and E,F,I are not defined symmetric to "
1071 "x-axis, which is required to interpolate block 2"
1072 " and 5 with same geometric computations."));
1073 const double l_y = D(1) - C(1);
1074 const double l_h = D(1) - l_y;
1075 const double by = -l_h / length_b1_x * (node(0) - H(0));
1076 const double dy = (height - by) / n_cells_y;
1077 const int iy = static_cast<int>(
1078 std::round((std::abs(node(1)) - by) / dy));
1079 const double dx = length_b1_x / n_cells_x_1;
1080 const int ix = static_cast<int>(
1081 std::round(std::abs(node(0) - center_mesh) / dx));
1082
1083 const double alpha = bias_alpha(1 - (1.0 * iy) / n_cells_y);
1084 // define points on upper/lower horizontal far field side,
1085 // i.e. face DG or FI. Incline factor to move points G and I
1086 // to the right by distance incline_facor*lenght_b2
1087 const Point<2> p(ix * dx + center_mesh +
1088 incline_factor * length_b2 * ix /
1089 n_cells_x_1,
1090 ((cell.material_id() == id_block_2) ?
1091 (height) :
1092 (-height)));
1093 // interpolate between y = height and upper airfoil points
1094 // (block2) or y = -height and lower airfoil points (block5)
1095 node = airfoil_1D[(
1096 (cell.material_id() == id_block_2) ? (0) : (1))]
1097 [n_cells_x_0 + ix] *
1098 alpha +
1099 p * (1 - alpha);
1100 }
1101 else if (cell.material_id() == id_block_3 ||
1102 cell.material_id() == id_block_6) // block 3 and 6
1103 {
1104 // compute indices ix and iy
1105 const double dx = length_b2 / n_cells_x_2;
1106 const double dy = height / n_cells_y;
1107 const int ix = static_cast<int>(
1108 std::round(std::abs(node(0) - H(0)) / dx));
1109 const int iy =
1110 static_cast<int>(std::round(std::abs(node(1)) / dy));
1111
1112 const double alpha_y = bias_alpha(1 - 1.0 * iy / n_cells_y);
1113 const double alpha_x =
1114 bias_alpha(1 - (static_cast<double>(ix)) / n_cells_x_2);
1115 // define on upper/lower horizontal far field side at y =
1116 // +/- height, i.e. face GK or IL incline factor to move
1117 // points G and H to the right
1118 const Point<2> p1(J(0) - (1 - incline_factor) * length_b2 *
1119 (alpha_x),
1120 ((cell.material_id() == id_block_3) ?
1121 (height) :
1122 (-height)));
1123 // define points on HJ but use tail_y as y-coordinate, in
1124 // case last airfoil point has y =/= 0
1125 const Point<2> p2(J(0) - alpha_x * length_b2, tail_y);
1126 node = p1 * (1 - alpha_y) + p2 * alpha_y;
1127 }
1128 else
1129 {
1130 Assert(false,
1131 ExcIndexRange(cell.material_id(),
1132 id_block_1,
1133 id_block_6));
1134 }
1135 }
1136 }
1137
1138
1139 /*
1140 * This function returns a bias factor 'alpha' which is used to make the
1141 * mesh more tight in close distance of the airfoil.
1142 * It is a bijective function mapping from [0,1] onto [0,1] where values
1143 * near 1 are made tighter.
1144 */
1145 double
1146 bias_alpha(double alpha) const
1147 {
1148 return std::tanh(bias_factor * alpha) / std::tanh(bias_factor);
1149 }
1150 };
1151 } // namespace
1152
1153
1154
1155 void
1156 internal_create_triangulation(
1158 std::vector<GridTools::PeriodicFacePair<
1159 typename Triangulation<2, 2>::cell_iterator>> *periodic_faces,
1160 const AdditionalData & additional_data)
1161 {
1162 MeshGenerator mesh_generator(additional_data);
1163 // Cast the triangulation to the right type so that the right
1164 // specialization of the function create_triangulation is picked up.
1165 if (auto parallel_tria =
1167 &tria))
1168 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1169 else if (auto parallel_tria = dynamic_cast<
1171 &tria))
1172 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1173 else
1174 mesh_generator.create_triangulation(tria, periodic_faces);
1175 }
1176
1177 template <>
1178 void
1179 create_triangulation(Triangulation<1, 1> &, const AdditionalData &)
1180 {
1181 Assert(false, ExcMessage("Airfoils only exist for 2D and 3D!"));
1182 }
1183
1184
1185
1186 template <>
1187 void
1189 std::vector<GridTools::PeriodicFacePair<
1191 const AdditionalData &)
1192 {
1193 Assert(false, ExcMessage("Airfoils only exist for 2D and 3D!"));
1194 }
1195
1196
1197
1198 template <>
1199 void
1201 const AdditionalData &additional_data)
1202 {
1203 internal_create_triangulation(tria, nullptr, additional_data);
1204 }
1205
1206
1207
1208 template <>
1209 void
1212 std::vector<GridTools::PeriodicFacePair<
1213 typename Triangulation<2, 2>::cell_iterator>> &periodic_faces,
1214 const AdditionalData & additional_data)
1215 {
1216 internal_create_triangulation(tria, &periodic_faces, additional_data);
1217 }
1218
1219
1220
1221 template <>
1222 void
1225 std::vector<GridTools::PeriodicFacePair<
1226 typename Triangulation<3, 3>::cell_iterator>> &periodic_faces,
1227 const AdditionalData & additional_data)
1228 {
1229 Assert(false, ExcMessage("3D airfoils are not implemented yet!"));
1230 (void)tria;
1231 (void)additional_data;
1232 (void)periodic_faces;
1233 }
1234 } // namespace Airfoil
1235
1236
1237 namespace
1238 {
1243 template <int dim, int spacedim>
1244 void
1245 colorize_hyper_rectangle(Triangulation<dim, spacedim> &tria)
1246 {
1247 // there is nothing to do in 1d
1248 if (dim > 1)
1249 {
1250 // there is only one cell, so
1251 // simple task
1253 tria.begin();
1254 for (auto f : GeometryInfo<dim>::face_indices())
1255 cell->face(f)->set_boundary_id(f);
1256 }
1257 }
1258
1259
1260
1261 template <int spacedim>
1262 void
1263 colorize_subdivided_hyper_rectangle(Triangulation<1, spacedim> &tria,
1264 const Point<spacedim> &,
1265 const Point<spacedim> &,
1266 const double)
1267 {
1269 tria.begin();
1270 cell != tria.end();
1271 ++cell)
1272 if (cell->center()(0) > 0)
1273 cell->set_material_id(1);
1274 // boundary indicators are set to
1275 // 0 (left) and 1 (right) by default.
1276 }
1277
1278
1279
1280 template <int dim, int spacedim>
1281 void
1282 colorize_subdivided_hyper_rectangle(Triangulation<dim, spacedim> &tria,
1283 const Point<spacedim> & p1,
1284 const Point<spacedim> & p2,
1285 const double epsilon)
1286 {
1287 // run through all faces and check
1288 // if one of their center coordinates matches
1289 // one of the corner points. Comparisons
1290 // are made using an epsilon which
1291 // should be smaller than the smallest cell
1292 // diameter.
1293
1295 tria.begin_face(),
1296 endface =
1297 tria.end_face();
1298 for (; face != endface; ++face)
1299 if (face->at_boundary())
1300 if (face->boundary_id() == 0)
1301 {
1302 const Point<spacedim> center(face->center());
1303
1304 if (std::abs(center(0) - p1[0]) < epsilon)
1305 face->set_boundary_id(0);
1306 else if (std::abs(center(0) - p2[0]) < epsilon)
1307 face->set_boundary_id(1);
1308 else if (dim > 1 && std::abs(center(1) - p1[1]) < epsilon)
1309 face->set_boundary_id(2);
1310 else if (dim > 1 && std::abs(center(1) - p2[1]) < epsilon)
1311 face->set_boundary_id(3);
1312 else if (dim > 2 && std::abs(center(2) - p1[2]) < epsilon)
1313 face->set_boundary_id(4);
1314 else if (dim > 2 && std::abs(center(2) - p2[2]) < epsilon)
1315 face->set_boundary_id(5);
1316 else
1317 // triangulation says it
1318 // is on the boundary,
1319 // but we could not find
1320 // on which boundary.
1321 Assert(false, ExcInternalError());
1322 }
1323
1324 for (const auto &cell : tria.cell_iterators())
1325 {
1326 types::material_id id = 0;
1327 for (unsigned int d = 0; d < dim; ++d)
1328 if (cell->center()(d) > 0)
1329 id += (1 << d);
1330 cell->set_material_id(id);
1331 }
1332 }
1333
1334
1339 void
1340 colorize_hyper_shell(Triangulation<2> &tria,
1341 const Point<2> &,
1342 const double,
1343 const double)
1344 {
1345 // In spite of receiving geometrical
1346 // data, we do this only based on
1347 // topology.
1348
1349 // For the mesh based on cube,
1350 // this is highly irregular
1352 cell != tria.end();
1353 ++cell)
1354 {
1355 Assert(cell->face(2)->at_boundary(), ExcInternalError());
1356 cell->face(2)->set_all_boundary_ids(1);
1357 }
1358 }
1359
1360
1365 void
1366 colorize_hyper_shell(Triangulation<3> &tria,
1367 const Point<3> &,
1368 const double,
1369 const double)
1370 {
1371 // the following uses a good amount
1372 // of knowledge about the
1373 // orientation of cells. this is
1374 // probably not good style...
1375 if (tria.n_cells() == 6)
1376 {
1378
1379 Assert(cell->face(4)->at_boundary(), ExcInternalError());
1380 cell->face(4)->set_all_boundary_ids(1);
1381
1382 ++cell;
1383 Assert(cell->face(2)->at_boundary(), ExcInternalError());
1384 cell->face(2)->set_all_boundary_ids(1);
1385
1386 ++cell;
1387 Assert(cell->face(2)->at_boundary(), ExcInternalError());
1388 cell->face(2)->set_all_boundary_ids(1);
1389
1390 ++cell;
1391 Assert(cell->face(0)->at_boundary(), ExcInternalError());
1392 cell->face(0)->set_all_boundary_ids(1);
1393
1394 ++cell;
1395 Assert(cell->face(2)->at_boundary(), ExcInternalError());
1396 cell->face(2)->set_all_boundary_ids(1);
1397
1398 ++cell;
1399 Assert(cell->face(0)->at_boundary(), ExcInternalError());
1400 cell->face(0)->set_all_boundary_ids(1);
1401 }
1402 else if (tria.n_cells() == 12)
1403 {
1404 // again use some internal
1405 // knowledge
1407 cell != tria.end();
1408 ++cell)
1409 {
1410 Assert(cell->face(5)->at_boundary(), ExcInternalError());
1411 cell->face(5)->set_all_boundary_ids(1);
1412 }
1413 }
1414 else if (tria.n_cells() == 96)
1415 {
1416 // the 96-cell hypershell is
1417 // based on a once refined
1418 // 12-cell mesh. consequently,
1419 // since the outer faces all
1420 // are face_no==5 above, so
1421 // they are here (unless they
1422 // are in the interior). Use
1423 // this to assign boundary
1424 // indicators, but also make
1425 // sure that we encounter
1426 // exactly 48 such faces
1427 unsigned int count = 0;
1429 cell != tria.end();
1430 ++cell)
1431 if (cell->face(5)->at_boundary())
1432 {
1433 cell->face(5)->set_all_boundary_ids(1);
1434 ++count;
1435 }
1436 Assert(count == 48, ExcInternalError());
1437 }
1438 else
1439 Assert(false, ExcNotImplemented());
1440 }
1441
1442
1443
1449 void
1450 colorize_quarter_hyper_shell(Triangulation<3> &tria,
1451 const Point<3> & center,
1452 const double inner_radius,
1453 const double outer_radius)
1454 {
1455 if (tria.n_cells() != 3)
1457
1458 double middle = (outer_radius - inner_radius) / 2e0 + inner_radius;
1459 double eps = 1e-3 * middle;
1461
1462 for (; cell != tria.end(); ++cell)
1463 for (unsigned int f : GeometryInfo<3>::face_indices())
1464 {
1465 if (!cell->face(f)->at_boundary())
1466 continue;
1467
1468 double radius = cell->face(f)->center().norm() - center.norm();
1469 if (std::fabs(cell->face(f)->center()(0)) <
1470 eps) // x = 0 set boundary 2
1471 {
1472 cell->face(f)->set_boundary_id(2);
1473 for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1474 ++j)
1475 if (cell->face(f)->line(j)->at_boundary())
1476 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1477 cell->face(f)->line(j)->vertex(1).norm()) >
1478 eps)
1479 cell->face(f)->line(j)->set_boundary_id(2);
1480 }
1481 else if (std::fabs(cell->face(f)->center()(1)) <
1482 eps) // y = 0 set boundary 3
1483 {
1484 cell->face(f)->set_boundary_id(3);
1485 for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1486 ++j)
1487 if (cell->face(f)->line(j)->at_boundary())
1488 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1489 cell->face(f)->line(j)->vertex(1).norm()) >
1490 eps)
1491 cell->face(f)->line(j)->set_boundary_id(3);
1492 }
1493 else if (std::fabs(cell->face(f)->center()(2)) <
1494 eps) // z = 0 set boundary 4
1495 {
1496 cell->face(f)->set_boundary_id(4);
1497 for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1498 ++j)
1499 if (cell->face(f)->line(j)->at_boundary())
1500 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1501 cell->face(f)->line(j)->vertex(1).norm()) >
1502 eps)
1503 cell->face(f)->line(j)->set_boundary_id(4);
1504 }
1505 else if (radius < middle) // inner radius set boundary 0
1506 {
1507 cell->face(f)->set_boundary_id(0);
1508 for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1509 ++j)
1510 if (cell->face(f)->line(j)->at_boundary())
1511 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1512 cell->face(f)->line(j)->vertex(1).norm()) <
1513 eps)
1514 cell->face(f)->line(j)->set_boundary_id(0);
1515 }
1516 else if (radius > middle) // outer radius set boundary 1
1517 {
1518 cell->face(f)->set_boundary_id(1);
1519 for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1520 ++j)
1521 if (cell->face(f)->line(j)->at_boundary())
1522 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1523 cell->face(f)->line(j)->vertex(1).norm()) <
1524 eps)
1525 cell->face(f)->line(j)->set_boundary_id(1);
1526 }
1527 else
1528 Assert(false, ExcInternalError());
1529 }
1530 }
1531
1532 } // namespace
1533
1534
1535 template <int dim, int spacedim>
1536 void
1538 const Point<dim> & p_1,
1539 const Point<dim> & p_2,
1540 const bool colorize)
1541 {
1542 // First, extend dimensions from dim to spacedim and
1543 // normalize such that p1 is lower in all coordinate
1544 // directions. Additional entries will be 0.
1545 Point<spacedim> p1, p2;
1546 for (unsigned int i = 0; i < dim; ++i)
1547 {
1548 p1(i) = std::min(p_1(i), p_2(i));
1549 p2(i) = std::max(p_1(i), p_2(i));
1550 }
1551
1552 std::vector<Point<spacedim>> vertices(GeometryInfo<dim>::vertices_per_cell);
1553 switch (dim)
1554 {
1555 case 1:
1556 vertices[0] = p1;
1557 vertices[1] = p2;
1558 break;
1559 case 2:
1560 vertices[0] = vertices[1] = p1;
1561 vertices[2] = vertices[3] = p2;
1562
1563 vertices[1](0) = p2(0);
1564 vertices[2](0) = p1(0);
1565 break;
1566 case 3:
1567 vertices[0] = vertices[1] = vertices[2] = vertices[3] = p1;
1568 vertices[4] = vertices[5] = vertices[6] = vertices[7] = p2;
1569
1570 vertices[1](0) = p2(0);
1571 vertices[2](1) = p2(1);
1572 vertices[3](0) = p2(0);
1573 vertices[3](1) = p2(1);
1574
1575 vertices[4](0) = p1(0);
1576 vertices[4](1) = p1(1);
1577 vertices[5](1) = p1(1);
1578 vertices[6](0) = p1(0);
1579
1580 break;
1581 default:
1582 Assert(false, ExcNotImplemented());
1583 }
1584
1585 // Prepare cell data
1586 std::vector<CellData<dim>> cells(1);
1587 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1588 cells[0].vertices[i] = i;
1589 cells[0].material_id = 0;
1590
1592
1593 // Assign boundary indicators
1594 if (colorize)
1595 colorize_hyper_rectangle(tria);
1596 }
1597
1598
1599
1600 template <int dim, int spacedim>
1601 void
1603 const double left,
1604 const double right,
1605 const bool colorize)
1606 {
1607 Assert(left < right,
1608 ExcMessage("Invalid left-to-right bounds of hypercube"));
1609
1610 Point<dim> p1, p2;
1611 for (unsigned int i = 0; i < dim; ++i)
1612 {
1613 p1(i) = left;
1614 p2(i) = right;
1615 }
1616 hyper_rectangle(tria, p1, p2, colorize);
1617 }
1618
1619
1620
1621 template <int dim>
1622 void
1623 simplex(Triangulation<dim> &tria, const std::vector<Point<dim>> &vertices)
1624 {
1625 AssertDimension(vertices.size(), dim + 1);
1626 Assert(dim > 1, ExcNotImplemented());
1627 Assert(dim < 4, ExcNotImplemented());
1628
1629# ifdef DEBUG
1630 Tensor<2, dim> vector_matrix;
1631 for (unsigned int d = 0; d < dim; ++d)
1632 for (unsigned int c = 1; c <= dim; ++c)
1633 vector_matrix[c - 1][d] = vertices[c](d) - vertices[0](d);
1634 Assert(determinant(vector_matrix) > 0.,
1635 ExcMessage("Vertices of simplex must form a right handed system"));
1636# endif
1637
1638 // Set up the vertices by first copying into points.
1639 std::vector<Point<dim>> points = vertices;
1641 // Compute the edge midpoints and add up everything to compute the
1642 // center point.
1643 for (unsigned int i = 0; i <= dim; ++i)
1644 {
1645 points.push_back(0.5 * (points[i] + points[(i + 1) % (dim + 1)]));
1646 center += points[i];
1647 }
1648 if (dim > 2)
1649 {
1650 // In 3D, we have some more edges to deal with
1651 for (unsigned int i = 1; i < dim; ++i)
1652 points.push_back(0.5 * (points[i - 1] + points[i + 1]));
1653 // And we need face midpoints
1654 for (unsigned int i = 0; i <= dim; ++i)
1655 points.push_back(1. / 3. *
1656 (points[i] + points[(i + 1) % (dim + 1)] +
1657 points[(i + 2) % (dim + 1)]));
1658 }
1659 points.push_back((1. / (dim + 1)) * center);
1660
1661 std::vector<CellData<dim>> cells(dim + 1);
1662 switch (dim)
1663 {
1664 case 2:
1665 AssertDimension(points.size(), 7);
1666 cells[0].vertices[0] = 0;
1667 cells[0].vertices[1] = 3;
1668 cells[0].vertices[2] = 5;
1669 cells[0].vertices[3] = 6;
1670 cells[0].material_id = 0;
1671
1672 cells[1].vertices[0] = 3;
1673 cells[1].vertices[1] = 1;
1674 cells[1].vertices[2] = 6;
1675 cells[1].vertices[3] = 4;
1676 cells[1].material_id = 0;
1677
1678 cells[2].vertices[0] = 5;
1679 cells[2].vertices[1] = 6;
1680 cells[2].vertices[2] = 2;
1681 cells[2].vertices[3] = 4;
1682 cells[2].material_id = 0;
1683 break;
1684 case 3:
1685 AssertDimension(points.size(), 15);
1686 cells[0].vertices[0] = 0;
1687 cells[0].vertices[1] = 4;
1688 cells[0].vertices[2] = 8;
1689 cells[0].vertices[3] = 10;
1690 cells[0].vertices[4] = 7;
1691 cells[0].vertices[5] = 13;
1692 cells[0].vertices[6] = 12;
1693 cells[0].vertices[7] = 14;
1694 cells[0].material_id = 0;
1695
1696 cells[1].vertices[0] = 4;
1697 cells[1].vertices[1] = 1;
1698 cells[1].vertices[2] = 10;
1699 cells[1].vertices[3] = 5;
1700 cells[1].vertices[4] = 13;
1701 cells[1].vertices[5] = 9;
1702 cells[1].vertices[6] = 14;
1703 cells[1].vertices[7] = 11;
1704 cells[1].material_id = 0;
1705
1706 cells[2].vertices[0] = 8;
1707 cells[2].vertices[1] = 10;
1708 cells[2].vertices[2] = 2;
1709 cells[2].vertices[3] = 5;
1710 cells[2].vertices[4] = 12;
1711 cells[2].vertices[5] = 14;
1712 cells[2].vertices[6] = 6;
1713 cells[2].vertices[7] = 11;
1714 cells[2].material_id = 0;
1715
1716 cells[3].vertices[0] = 7;
1717 cells[3].vertices[1] = 13;
1718 cells[3].vertices[2] = 12;
1719 cells[3].vertices[3] = 14;
1720 cells[3].vertices[4] = 3;
1721 cells[3].vertices[5] = 9;
1722 cells[3].vertices[6] = 6;
1723 cells[3].vertices[7] = 11;
1724 cells[3].material_id = 0;
1725 break;
1726 default:
1727 Assert(false, ExcNotImplemented());
1728 }
1729 tria.create_triangulation(points, cells, SubCellData());
1730 }
1731
1732
1733
1734 template <int dim, int spacedim>
1735 void
1738 {
1739 AssertDimension(dim, reference_cell.get_dimension());
1740
1741 if (reference_cell == ReferenceCells::get_hypercube<dim>())
1742 {
1744 }
1745 else
1746 {
1747 // Create an array that contains the vertices of the reference cell.
1748 // We can query these points from ReferenceCell, but then we have
1749 // to embed them into the spacedim-dimensional space.
1750 std::vector<Point<spacedim>> vertices(reference_cell.n_vertices());
1751 for (const unsigned int v : reference_cell.vertex_indices())
1752 {
1753 const Point<dim> this_vertex = reference_cell.vertex<dim>(v);
1754 for (unsigned int d = 0; d < dim; ++d)
1755 vertices[v][d] = this_vertex[d];
1756 // Point<spacedim> initializes everything to zero, so any remaining
1757 // elements are left at zero and we don't have to explicitly pad
1758 // from 'dim' to 'spacedim' here.
1759 }
1760
1761 // Then make one cell out of these vertices. They are ordered correctly
1762 // already, so we just need to enumerate them
1763 std::vector<CellData<dim>> cells(1);
1764 cells[0].vertices.resize(reference_cell.n_vertices());
1765 for (const unsigned int v : reference_cell.vertex_indices())
1766 cells[0].vertices[v] = v;
1767
1768 // Turn all of this into a triangulation
1770 }
1771 }
1772
1773 void
1775 const unsigned int n_cells,
1776 const unsigned int n_rotations,
1777 const double R,
1778 const double r)
1779 {
1780 const unsigned int dim = 3;
1781 Assert(n_cells > 4,
1782 ExcMessage(
1783 "More than 4 cells are needed to create a moebius grid."));
1784 Assert(r > 0 && R > 0,
1785 ExcMessage("Outer and inner radius must be positive."));
1786 Assert(R > r,
1787 ExcMessage("Outer radius must be greater than inner radius."));
1788
1789
1790 std::vector<Point<dim>> vertices(4 * n_cells);
1791 double beta_step = n_rotations * numbers::PI / 2.0 / n_cells;
1792 double alpha_step = 2.0 * numbers::PI / n_cells;
1793
1794 for (unsigned int i = 0; i < n_cells; ++i)
1795 for (unsigned int j = 0; j < 4; ++j)
1796 {
1797 vertices[4 * i + j][0] =
1798 R * std::cos(i * alpha_step) +
1799 r * std::cos(i * beta_step + j * numbers::PI / 2.0) *
1800 std::cos(i * alpha_step);
1801 vertices[4 * i + j][1] =
1802 R * std::sin(i * alpha_step) +
1803 r * std::cos(i * beta_step + j * numbers::PI / 2.0) *
1804 std::sin(i * alpha_step);
1805 vertices[4 * i + j][2] =
1806 r * std::sin(i * beta_step + j * numbers::PI / 2.0);
1807 }
1808
1809 unsigned int offset = 0;
1810
1811 // This Triangulation is constructed using the UCD numbering scheme since,
1812 // in that numbering, the front face is first and the back face is second,
1813 // which is more convenient for creating a moebius
1814 std::vector<CellData<dim>> cells(n_cells);
1815 for (unsigned int i = 0; i < n_cells; ++i)
1816 {
1817 for (unsigned int j = 0; j < 2; ++j)
1818 {
1819 cells[i].vertices[GeometryInfo<3>::ucd_to_deal[0 + 4 * j]] =
1820 offset + 0 + 4 * j;
1821 cells[i].vertices[GeometryInfo<3>::ucd_to_deal[1 + 4 * j]] =
1822 offset + 3 + 4 * j;
1823 cells[i].vertices[GeometryInfo<3>::ucd_to_deal[2 + 4 * j]] =
1824 offset + 2 + 4 * j;
1825 cells[i].vertices[GeometryInfo<3>::ucd_to_deal[3 + 4 * j]] =
1826 offset + 1 + 4 * j;
1827 }
1828 offset += 4;
1829 cells[i].material_id = 0;
1830 }
1831
1832 // now correct the last four vertices
1833 cells[n_cells - 1].vertices[GeometryInfo<3>::ucd_to_deal[4]] =
1834 (0 + n_rotations) % 4;
1835 cells[n_cells - 1].vertices[GeometryInfo<3>::ucd_to_deal[5]] =
1836 (3 + n_rotations) % 4;
1837 cells[n_cells - 1].vertices[GeometryInfo<3>::ucd_to_deal[6]] =
1838 (2 + n_rotations) % 4;
1839 cells[n_cells - 1].vertices[GeometryInfo<3>::ucd_to_deal[7]] =
1840 (1 + n_rotations) % 4;
1841
1844 }
1845
1846
1847
1848 template <>
1849 void
1850 torus<2, 3>(Triangulation<2, 3> &tria,
1851 const double R,
1852 const double r,
1853 const unsigned int,
1854 const double)
1855 {
1856 Assert(R > r,
1857 ExcMessage("Outer radius R must be greater than the inner "
1858 "radius r."));
1859 Assert(r > 0.0, ExcMessage("The inner radius r must be positive."));
1860
1861 const unsigned int dim = 2;
1862 const unsigned int spacedim = 3;
1863 std::vector<Point<spacedim>> vertices(16);
1864
1865 vertices[0] = Point<spacedim>(R - r, 0, 0);
1866 vertices[1] = Point<spacedim>(R, -r, 0);
1867 vertices[2] = Point<spacedim>(R + r, 0, 0);
1868 vertices[3] = Point<spacedim>(R, r, 0);
1869 vertices[4] = Point<spacedim>(0, 0, R - r);
1870 vertices[5] = Point<spacedim>(0, -r, R);
1871 vertices[6] = Point<spacedim>(0, 0, R + r);
1872 vertices[7] = Point<spacedim>(0, r, R);
1873 vertices[8] = Point<spacedim>(-(R - r), 0, 0);
1874 vertices[9] = Point<spacedim>(-R, -r, 0);
1875 vertices[10] = Point<spacedim>(-(R + r), 0, 0);
1876 vertices[11] = Point<spacedim>(-R, r, 0);
1877 vertices[12] = Point<spacedim>(0, 0, -(R - r));
1878 vertices[13] = Point<spacedim>(0, -r, -R);
1879 vertices[14] = Point<spacedim>(0, 0, -(R + r));
1880 vertices[15] = Point<spacedim>(0, r, -R);
1881
1882 std::vector<CellData<dim>> cells(16);
1883 // Right Hand Orientation
1884 cells[0].vertices[0] = 0;
1885 cells[0].vertices[1] = 4;
1886 cells[0].vertices[2] = 3;
1887 cells[0].vertices[3] = 7;
1888 cells[0].material_id = 0;
1889
1890 cells[1].vertices[0] = 1;
1891 cells[1].vertices[1] = 5;
1892 cells[1].vertices[2] = 0;
1893 cells[1].vertices[3] = 4;
1894 cells[1].material_id = 0;
1895
1896 cells[2].vertices[0] = 2;
1897 cells[2].vertices[1] = 6;
1898 cells[2].vertices[2] = 1;
1899 cells[2].vertices[3] = 5;
1900 cells[2].material_id = 0;
1901
1902 cells[3].vertices[0] = 3;
1903 cells[3].vertices[1] = 7;
1904 cells[3].vertices[2] = 2;
1905 cells[3].vertices[3] = 6;
1906 cells[3].material_id = 0;
1907
1908 cells[4].vertices[0] = 4;
1909 cells[4].vertices[1] = 8;
1910 cells[4].vertices[2] = 7;
1911 cells[4].vertices[3] = 11;
1912 cells[4].material_id = 0;
1913
1914 cells[5].vertices[0] = 5;
1915 cells[5].vertices[1] = 9;
1916 cells[5].vertices[2] = 4;
1917 cells[5].vertices[3] = 8;
1918 cells[5].material_id = 0;
1919
1920 cells[6].vertices[0] = 6;
1921 cells[6].vertices[1] = 10;
1922 cells[6].vertices[2] = 5;
1923 cells[6].vertices[3] = 9;
1924 cells[6].material_id = 0;
1925
1926 cells[7].vertices[0] = 7;
1927 cells[7].vertices[1] = 11;
1928 cells[7].vertices[2] = 6;
1929 cells[7].vertices[3] = 10;
1930 cells[7].material_id = 0;
1931
1932 cells[8].vertices[0] = 8;
1933 cells[8].vertices[1] = 12;
1934 cells[8].vertices[2] = 11;
1935 cells[8].vertices[3] = 15;
1936 cells[8].material_id = 0;
1937
1938 cells[9].vertices[0] = 9;
1939 cells[9].vertices[1] = 13;
1940 cells[9].vertices[2] = 8;
1941 cells[9].vertices[3] = 12;
1942 cells[9].material_id = 0;
1943
1944 cells[10].vertices[0] = 10;
1945 cells[10].vertices[1] = 14;
1946 cells[10].vertices[2] = 9;
1947 cells[10].vertices[3] = 13;
1948 cells[10].material_id = 0;
1949
1950 cells[11].vertices[0] = 11;
1951 cells[11].vertices[1] = 15;
1952 cells[11].vertices[2] = 10;
1953 cells[11].vertices[3] = 14;
1954 cells[11].material_id = 0;
1955
1956 cells[12].vertices[0] = 12;
1957 cells[12].vertices[1] = 0;
1958 cells[12].vertices[2] = 15;
1959 cells[12].vertices[3] = 3;
1960 cells[12].material_id = 0;
1961
1962 cells[13].vertices[0] = 13;
1963 cells[13].vertices[1] = 1;
1964 cells[13].vertices[2] = 12;
1965 cells[13].vertices[3] = 0;
1966 cells[13].material_id = 0;
1967
1968 cells[14].vertices[0] = 14;
1969 cells[14].vertices[1] = 2;
1970 cells[14].vertices[2] = 13;
1971 cells[14].vertices[3] = 1;
1972 cells[14].material_id = 0;
1973
1974 cells[15].vertices[0] = 15;
1975 cells[15].vertices[1] = 3;
1976 cells[15].vertices[2] = 14;
1977 cells[15].vertices[3] = 2;
1978 cells[15].material_id = 0;
1979
1982
1985 }
1986
1987
1988
1989 template <>
1990 void
1991 torus<3, 3>(Triangulation<3, 3> &tria,
1992 const double R,
1993 const double r,
1994 const unsigned int n_cells_toroidal,
1995 const double phi)
1996 {
1997 Assert(R > r,
1998 ExcMessage("Outer radius R must be greater than the inner "
1999 "radius r."));
2000 Assert(r > 0.0, ExcMessage("The inner radius r must be positive."));
2001 Assert(n_cells_toroidal > 2,
2002 ExcMessage("Number of cells in toroidal direction has "
2003 "to be at least 3."));
2004 AssertThrow(phi > 0.0 && phi < 2.0 * numbers::PI + 1.0e-15,
2005 ExcMessage("Invalid angle phi specified."));
2006
2007 // the first 8 vertices are in the x-y-plane
2008 Point<3> const p = Point<3>(R, 0.0, 0.0);
2009 double const a = 1. / (1 + std::sqrt(2.0));
2010 // A value of 1 indicates "open" torus with angle < 2*pi, which
2011 // means that we need an additional layer of vertices
2012 const unsigned int additional_layer =
2013 (phi < 2.0 * numbers::PI - 1.0e-15) ?
2014 1 :
2015 0; // torus is closed (angle of 2*pi)
2016 const unsigned int n_point_layers_toroidal =
2017 n_cells_toroidal + additional_layer;
2018 std::vector<Point<3>> vertices(8 * n_point_layers_toroidal);
2019 vertices[0] = p + Point<3>(-1, -1, 0) * (r / std::sqrt(2.0)),
2020 vertices[1] = p + Point<3>(+1, -1, 0) * (r / std::sqrt(2.0)),
2021 vertices[2] = p + Point<3>(-1, -1, 0) * (r / std::sqrt(2.0) * a),
2022 vertices[3] = p + Point<3>(+1, -1, 0) * (r / std::sqrt(2.0) * a),
2023 vertices[4] = p + Point<3>(-1, +1, 0) * (r / std::sqrt(2.0) * a),
2024 vertices[5] = p + Point<3>(+1, +1, 0) * (r / std::sqrt(2.0) * a),
2025 vertices[6] = p + Point<3>(-1, +1, 0) * (r / std::sqrt(2.0)),
2026 vertices[7] = p + Point<3>(+1, +1, 0) * (r / std::sqrt(2.0));
2027
2028 // create remaining vertices by rotating around negative y-axis (the
2029 // direction is to ensure positive cell measures)
2030 double const phi_cell = phi / n_cells_toroidal;
2031 for (unsigned int c = 1; c < n_point_layers_toroidal; ++c)
2032 {
2033 for (unsigned int v = 0; v < 8; ++v)
2034 {
2035 double const r_2d = vertices[v][0];
2036 vertices[8 * c + v][0] = r_2d * std::cos(phi_cell * c);
2037 vertices[8 * c + v][1] = vertices[v][1];
2038 vertices[8 * c + v][2] = r_2d * std::sin(phi_cell * c);
2039 }
2040 }
2041
2042 // cell connectivity
2043 std::vector<CellData<3>> cells(5 * n_cells_toroidal);
2044 for (unsigned int c = 0; c < n_cells_toroidal; ++c)
2045 {
2046 for (unsigned int j = 0; j < 2; ++j)
2047 {
2048 const unsigned int offset =
2049 (8 * (c + j)) % (8 * n_point_layers_toroidal);
2050
2051 // cell 0 in x-y-plane
2052 cells[5 * c].vertices[0 + j * 4] = offset + 0;
2053 cells[5 * c].vertices[1 + j * 4] = offset + 1;
2054 cells[5 * c].vertices[2 + j * 4] = offset + 2;
2055 cells[5 * c].vertices[3 + j * 4] = offset + 3;
2056 // cell 1 in x-y-plane (cell on torus centerline)
2057 cells[5 * c + 1].vertices[0 + j * 4] = offset + 2;
2058 cells[5 * c + 1].vertices[1 + j * 4] = offset + 3;
2059 cells[5 * c + 1].vertices[2 + j * 4] = offset + 4;
2060 cells[5 * c + 1].vertices[3 + j * 4] = offset + 5;
2061 // cell 2 in x-y-plane
2062 cells[5 * c + 2].vertices[0 + j * 4] = offset + 4;
2063 cells[5 * c + 2].vertices[1 + j * 4] = offset + 5;
2064 cells[5 * c + 2].vertices[2 + j * 4] = offset + 6;
2065 cells[5 * c + 2].vertices[3 + j * 4] = offset + 7;
2066 // cell 3 in x-y-plane
2067 cells[5 * c + 3].vertices[0 + j * 4] = offset + 0;
2068 cells[5 * c + 3].vertices[1 + j * 4] = offset + 2;
2069 cells[5 * c + 3].vertices[2 + j * 4] = offset + 6;
2070 cells[5 * c + 3].vertices[3 + j * 4] = offset + 4;
2071 // cell 4 in x-y-plane
2072 cells[5 * c + 4].vertices[0 + j * 4] = offset + 3;
2073 cells[5 * c + 4].vertices[1 + j * 4] = offset + 1;
2074 cells[5 * c + 4].vertices[2 + j * 4] = offset + 5;
2075 cells[5 * c + 4].vertices[3 + j * 4] = offset + 7;
2076 }
2077
2078 cells[5 * c].material_id = 0;
2079 // mark cell on torus centerline
2080 cells[5 * c + 1].material_id = 1;
2081 cells[5 * c + 2].material_id = 0;
2082 cells[5 * c + 3].material_id = 0;
2083 cells[5 * c + 4].material_id = 0;
2084 }
2085
2087
2090
2091 for (auto &cell : tria.cell_iterators())
2092 {
2093 // identify faces on torus surface and set manifold to 1
2094 for (unsigned int f : GeometryInfo<3>::face_indices())
2095 {
2096 // faces 4 and 5 are those with normal vector aligned with torus
2097 // centerline
2098 if (cell->face(f)->at_boundary() && f != 4 && f != 5)
2099 {
2100 cell->face(f)->set_all_manifold_ids(1);
2101 }
2102 }
2103
2104 // set manifold id to 2 for those cells that are on the torus centerline
2105 if (cell->material_id() == 1)
2106 {
2107 cell->set_all_manifold_ids(2);
2108 // reset to 0
2109 cell->set_material_id(0);
2110 }
2111 }
2112
2116 Point<3>()));
2118 transfinite.initialize(tria);
2119 tria.set_manifold(0, transfinite);
2120 }
2121
2122
2123
2124 template <int dim, int spacedim>
2125 void
2127 const std::vector<Point<spacedim>> &vertices,
2128 const bool colorize)
2129 {
2131 ExcMessage("Wrong number of vertices."));
2132
2133 // First create a hyper_rectangle and then deform it.
2134 hyper_cube(tria, 0, 1, colorize);
2135
2138 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
2139 cell->vertex(i) = vertices[i];
2140
2141 // Check that the order of the vertices makes sense, i.e., the volume of the
2142 // cell is positive.
2144 ExcMessage(
2145 "The volume of the cell is not greater than zero. "
2146 "This could be due to the wrong ordering of the vertices."));
2147 }
2148
2149
2150
2151 template <>
2152 void
2154 const Point<3> (&/*corners*/)[3],
2155 const bool /*colorize*/)
2156 {
2157 Assert(false, ExcNotImplemented());
2158 }
2159
2160 template <>
2161 void
2163 const Point<1> (&/*corners*/)[1],
2164 const bool /*colorize*/)
2165 {
2166 Assert(false, ExcNotImplemented());
2167 }
2168
2169 // Implementation for 2D only
2170 template <>
2171 void
2173 const Point<2> (&corners)[2],
2174 const bool colorize)
2175 {
2176 Point<2> origin;
2177 std::array<Tensor<1, 2>, 2> edges;
2178 edges[0] = corners[0];
2179 edges[1] = corners[1];
2180 std::vector<unsigned int> subdivisions;
2181 subdivided_parallelepiped<2, 2>(
2182 tria, origin, edges, subdivisions, colorize);
2183 }
2184
2185
2186
2187 template <int dim>
2188 void
2190 const Point<dim> (&corners)[dim],
2191 const bool colorize)
2192 {
2193 unsigned int n_subdivisions[dim];
2194 for (unsigned int i = 0; i < dim; ++i)
2195 n_subdivisions[i] = 1;
2196
2197 // and call the function below
2198 subdivided_parallelepiped(tria, n_subdivisions, corners, colorize);
2199 }
2200
2201 template <int dim>
2202 void
2204 const unsigned int n_subdivisions,
2205 const Point<dim> (&corners)[dim],
2206 const bool colorize)
2207 {
2208 // Equalize number of subdivisions in each dim-direction, their
2209 // validity will be checked later
2210 unsigned int n_subdivisions_[dim];
2211 for (unsigned int i = 0; i < dim; ++i)
2212 n_subdivisions_[i] = n_subdivisions;
2213
2214 // and call the function below
2215 subdivided_parallelepiped(tria, n_subdivisions_, corners, colorize);
2216 }
2217
2218 template <int dim>
2219 void
2221# ifndef _MSC_VER
2222 const unsigned int (&n_subdivisions)[dim],
2223# else
2224 const unsigned int *n_subdivisions,
2225# endif
2226 const Point<dim> (&corners)[dim],
2227 const bool colorize)
2228 {
2229 Point<dim> origin;
2230 std::vector<unsigned int> subdivisions;
2231 std::array<Tensor<1, dim>, dim> edges;
2232 for (unsigned int i = 0; i < dim; ++i)
2233 {
2234 subdivisions.push_back(n_subdivisions[i]);
2235 edges[i] = corners[i];
2236 }
2237
2238 subdivided_parallelepiped<dim, dim>(
2239 tria, origin, edges, subdivisions, colorize);
2240 }
2241
2242 // Parallelepiped implementation in 1d, 2d, and 3d. @note The
2243 // implementation in 1d is similar to hyper_rectangle(), and in 2d is
2244 // similar to parallelogram().
2245 template <int dim, int spacedim>
2246 void
2248 const Point<spacedim> & origin,
2249 const std::array<Tensor<1, spacedim>, dim> &edges,
2250 const std::vector<unsigned int> &subdivisions,
2251 const bool colorize)
2252 {
2253 std::vector<unsigned int> compute_subdivisions = subdivisions;
2254 if (compute_subdivisions.size() == 0)
2255 {
2256 compute_subdivisions.resize(dim, 1);
2257 }
2258
2259 Assert(compute_subdivisions.size() == dim,
2260 ExcMessage("One subdivision must be provided for each dimension."));
2261 // check subdivisions
2262 for (unsigned int i = 0; i < dim; ++i)
2263 {
2264 Assert(compute_subdivisions[i] > 0,
2265 ExcInvalidRepetitions(subdivisions[i]));
2266 Assert(
2267 edges[i].norm() > 0,
2268 ExcMessage(
2269 "Edges in subdivided_parallelepiped() must not be degenerate."));
2270 }
2271
2272 /*
2273 * Verify that the edge points to the right in 1D, vectors are oriented in
2274 * a counter clockwise direction in 2D, or form a right handed system in
2275 * 3D.
2276 */
2277 bool twisted_data = false;
2278 switch (dim)
2279 {
2280 case 1:
2281 {
2282 twisted_data = (edges[0][0] < 0);
2283 break;
2284 }
2285 case 2:
2286 {
2287 if (spacedim == 2) // this check does not make sense otherwise
2288 {
2289 const double plane_normal =
2290 edges[0][0] * edges[1][1] - edges[0][1] * edges[1][0];
2291 twisted_data = (plane_normal < 0.0);
2292 }
2293 break;
2294 }
2295 case 3:
2296 {
2297 // Check that the first two vectors are not linear combinations to
2298 // avoid zero division later on.
2299 Assert(std::abs(edges[0] * edges[1] /
2300 (edges[0].norm() * edges[1].norm()) -
2301 1.0) > 1.0e-15,
2302 ExcMessage(
2303 "Edges in subdivided_parallelepiped() must point in"
2304 " different directions."));
2305 const Tensor<1, spacedim> plane_normal =
2306 cross_product_3d(edges[0], edges[1]);
2307
2308 /*
2309 * Ensure that edges 1, 2, and 3 form a right-handed set of
2310 * vectors. This works by applying the definition of the dot product
2311 *
2312 * cos(theta) = dot(x, y)/(norm(x)*norm(y))
2313 *
2314 * and then, since the normal vector and third edge should both
2315 * point away from the plane formed by the first two edges, the
2316 * angle between them must be between 0 and pi/2; hence we just need
2317 *
2318 * 0 < dot(x, y).
2319 */
2320 twisted_data = (plane_normal * edges[2] < 0.0);
2321 break;
2322 }
2323 default:
2324 Assert(false, ExcInternalError());
2325 }
2326 (void)twisted_data; // make the static analyzer happy
2327 Assert(
2328 !twisted_data,
2330 "The triangulation you are trying to create will consist of cells"
2331 " with negative measures. This is usually the result of input data"
2332 " that does not define a right-handed coordinate system. The usual"
2333 " fix for this is to ensure that in 1D the given point is to the"
2334 " right of the origin (or the given edge tensor is positive), in 2D"
2335 " that the two edges (and their cross product) obey the right-hand"
2336 " rule (which may usually be done by switching the order of the"
2337 " points or edge tensors), or in 3D that the edges form a"
2338 " right-handed coordinate system (which may also be accomplished by"
2339 " switching the order of the first two points or edge tensors)."));
2340
2341 // Check corners do not overlap (unique)
2342 for (unsigned int i = 0; i < dim; ++i)
2343 for (unsigned int j = i + 1; j < dim; ++j)
2344 Assert((edges[i] != edges[j]),
2345 ExcMessage(
2346 "Degenerate edges of subdivided_parallelepiped encountered."));
2347
2348 // Create a list of points
2349 std::vector<Point<spacedim>> points;
2350
2351 switch (dim)
2352 {
2353 case 1:
2354 for (unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2355 points.push_back(origin + edges[0] / compute_subdivisions[0] * x);
2356 break;
2357
2358 case 2:
2359 for (unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2360 for (unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2361 points.push_back(origin + edges[0] / compute_subdivisions[0] * x +
2362 edges[1] / compute_subdivisions[1] * y);
2363 break;
2364
2365 case 3:
2366 for (unsigned int z = 0; z <= compute_subdivisions[2]; ++z)
2367 for (unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2368 for (unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2369 points.push_back(origin +
2370 edges[0] / compute_subdivisions[0] * x +
2371 edges[1] / compute_subdivisions[1] * y +
2372 edges[2] / compute_subdivisions[2] * z);
2373 break;
2374
2375 default:
2376 Assert(false, ExcNotImplemented());
2377 }
2378
2379 // Prepare cell data
2380 unsigned int n_cells = 1;
2381 for (unsigned int i = 0; i < dim; ++i)
2382 n_cells *= compute_subdivisions[i];
2383 std::vector<CellData<dim>> cells(n_cells);
2384
2385 // Create fixed ordering of
2386 switch (dim)
2387 {
2388 case 1:
2389 for (unsigned int x = 0; x < compute_subdivisions[0]; ++x)
2390 {
2391 cells[x].vertices[0] = x;
2392 cells[x].vertices[1] = x + 1;
2393
2394 // wipe material id
2395 cells[x].material_id = 0;
2396 }
2397 break;
2398
2399 case 2:
2400 {
2401 // Shorthand
2402 const unsigned int n_dy = compute_subdivisions[1];
2403 const unsigned int n_dx = compute_subdivisions[0];
2404
2405 for (unsigned int y = 0; y < n_dy; ++y)
2406 for (unsigned int x = 0; x < n_dx; ++x)
2407 {
2408 const unsigned int c = y * n_dx + x;
2409 cells[c].vertices[0] = y * (n_dx + 1) + x;
2410 cells[c].vertices[1] = y * (n_dx + 1) + x + 1;
2411 cells[c].vertices[2] = (y + 1) * (n_dx + 1) + x;
2412 cells[c].vertices[3] = (y + 1) * (n_dx + 1) + x + 1;
2413
2414 // wipe material id
2415 cells[c].material_id = 0;
2416 }
2417 }
2418 break;
2419
2420 case 3:
2421 {
2422 // Shorthand
2423 const unsigned int n_dz = compute_subdivisions[2];
2424 const unsigned int n_dy = compute_subdivisions[1];
2425 const unsigned int n_dx = compute_subdivisions[0];
2426
2427 for (unsigned int z = 0; z < n_dz; ++z)
2428 for (unsigned int y = 0; y < n_dy; ++y)
2429 for (unsigned int x = 0; x < n_dx; ++x)
2430 {
2431 const unsigned int c = z * n_dy * n_dx + y * n_dx + x;
2432
2433 cells[c].vertices[0] =
2434 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2435 cells[c].vertices[1] =
2436 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x + 1;
2437 cells[c].vertices[2] =
2438 z * (n_dy + 1) * (n_dx + 1) + (y + 1) * (n_dx + 1) + x;
2439 cells[c].vertices[3] = z * (n_dy + 1) * (n_dx + 1) +
2440 (y + 1) * (n_dx + 1) + x + 1;
2441 cells[c].vertices[4] =
2442 (z + 1) * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2443 cells[c].vertices[5] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2444 y * (n_dx + 1) + x + 1;
2445 cells[c].vertices[6] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2446 (y + 1) * (n_dx + 1) + x;
2447 cells[c].vertices[7] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2448 (y + 1) * (n_dx + 1) + x + 1;
2449
2450 // wipe material id
2451 cells[c].material_id = 0;
2452 }
2453 break;
2454 }
2455
2456 default:
2457 Assert(false, ExcNotImplemented());
2458 }
2459
2460 // Create triangulation
2461 // reorder the cells to ensure that they satisfy the convention for
2462 // edge and face directions
2464 tria.create_triangulation(points, cells, SubCellData());
2465
2466 // Finally assign boundary indicators according to hyper_rectangle
2467 if (colorize)
2468 {
2471 endc = tria.end();
2472 for (; cell != endc; ++cell)
2473 {
2474 for (const unsigned int face : GeometryInfo<dim>::face_indices())
2475 {
2476 if (cell->face(face)->at_boundary())
2477 cell->face(face)->set_boundary_id(face);
2478 }
2479 }
2480 }
2481 }
2482
2483
2484 template <int dim, int spacedim>
2485 void
2487 const unsigned int repetitions,
2488 const double left,
2489 const double right,
2490 const bool colorize)
2491 {
2492 Assert(repetitions >= 1, ExcInvalidRepetitions(repetitions));
2493 Assert(left < right,
2494 ExcMessage("Invalid left-to-right bounds of hypercube"));
2495
2496 Point<dim> p0, p1;
2497 for (unsigned int i = 0; i < dim; ++i)
2498 {
2499 p0[i] = left;
2500 p1[i] = right;
2501 }
2502
2503 std::vector<unsigned int> reps(dim, repetitions);
2505 }
2506
2507
2508
2509 template <int dim, int spacedim>
2510 void
2512 const std::vector<unsigned int> &repetitions,
2513 const Point<dim> & p_1,
2514 const Point<dim> & p_2,
2515 const bool colorize)
2516 {
2517 Assert(repetitions.size() == dim, ExcInvalidRepetitionsDimension(dim));
2518
2519 // First, extend dimensions from dim to spacedim and
2520 // normalize such that p1 is lower in all coordinate
2521 // directions. Additional entries will be 0.
2522 Point<spacedim> p1, p2;
2523 for (unsigned int i = 0; i < dim; ++i)
2524 {
2525 p1(i) = std::min(p_1(i), p_2(i));
2526 p2(i) = std::max(p_1(i), p_2(i));
2527 }
2528
2529 // calculate deltas and validate input
2530 std::array<Point<spacedim>, dim> delta;
2531 for (unsigned int i = 0; i < dim; ++i)
2532 {
2533 Assert(repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
2534
2535 delta[i][i] = (p2[i] - p1[i]) / repetitions[i];
2536 Assert(
2537 delta[i][i] > 0.0,
2538 ExcMessage(
2539 "The first dim entries of coordinates of p1 and p2 need to be different."));
2540 }
2541
2542 // then generate the points
2543 std::vector<Point<spacedim>> points;
2544 switch (dim)
2545 {
2546 case 1:
2547 for (unsigned int x = 0; x <= repetitions[0]; ++x)
2548 points.push_back(p1 + x * delta[0]);
2549 break;
2550
2551 case 2:
2552 for (unsigned int y = 0; y <= repetitions[1]; ++y)
2553 for (unsigned int x = 0; x <= repetitions[0]; ++x)
2554 points.push_back(p1 + x * delta[0] + y * delta[1]);
2555 break;
2556
2557 case 3:
2558 for (unsigned int z = 0; z <= repetitions[2]; ++z)
2559 for (unsigned int y = 0; y <= repetitions[1]; ++y)
2560 for (unsigned int x = 0; x <= repetitions[0]; ++x)
2561 points.push_back(p1 + x * delta[0] + y * delta[1] +
2562 z * delta[2]);
2563 break;
2564
2565 default:
2566 Assert(false, ExcNotImplemented());
2567 }
2568
2569 // next create the cells
2570 std::vector<CellData<dim>> cells;
2571 switch (dim)
2572 {
2573 case 1:
2574 {
2575 cells.resize(repetitions[0]);
2576 for (unsigned int x = 0; x < repetitions[0]; ++x)
2577 {
2578 cells[x].vertices[0] = x;
2579 cells[x].vertices[1] = x + 1;
2580 cells[x].material_id = 0;
2581 }
2582 break;
2583 }
2584
2585 case 2:
2586 {
2587 cells.resize(repetitions[1] * repetitions[0]);
2588 for (unsigned int y = 0; y < repetitions[1]; ++y)
2589 for (unsigned int x = 0; x < repetitions[0]; ++x)
2590 {
2591 const unsigned int c = x + y * repetitions[0];
2592 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
2593 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2594 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2595 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
2596 cells[c].material_id = 0;
2597 }
2598 break;
2599 }
2600
2601 case 3:
2602 {
2603 const unsigned int n_x = (repetitions[0] + 1);
2604 const unsigned int n_xy =
2605 (repetitions[0] + 1) * (repetitions[1] + 1);
2606
2607 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
2608 for (unsigned int z = 0; z < repetitions[2]; ++z)
2609 for (unsigned int y = 0; y < repetitions[1]; ++y)
2610 for (unsigned int x = 0; x < repetitions[0]; ++x)
2611 {
2612 const unsigned int c = x + y * repetitions[0] +
2613 z * repetitions[0] * repetitions[1];
2614 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2615 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2616 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2617 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2618 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2619 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2620 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2621 cells[c].vertices[7] =
2622 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2623 cells[c].material_id = 0;
2624 }
2625 break;
2626 }
2627
2628 default:
2629 Assert(false, ExcNotImplemented());
2630 }
2631
2632 tria.create_triangulation(points, cells, SubCellData());
2633
2634 if (colorize)
2635 {
2636 // to colorize, run through all
2637 // faces of all cells and set
2638 // boundary indicator to the
2639 // correct value if it was 0.
2640
2641 // use a large epsilon to
2642 // compare numbers to avoid
2643 // roundoff problems.
2644 double epsilon = 10;
2645 for (unsigned int i = 0; i < dim; ++i)
2646 epsilon = std::min(epsilon, 0.01 * delta[i][i]);
2647 Assert(epsilon > 0,
2648 ExcMessage(
2649 "The distance between corner points must be positive."))
2650
2651 // actual code is external since
2652 // 1-D is different from 2/3D.
2653 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
2654 }
2655 }
2656
2657
2658
2659 template <int dim>
2660 void
2662 const std::vector<std::vector<double>> &step_sz,
2663 const Point<dim> & p_1,
2664 const Point<dim> & p_2,
2665 const bool colorize)
2666 {
2667 Assert(step_sz.size() == dim, ExcInvalidRepetitionsDimension(dim));
2668
2669 // First, normalize input such that
2670 // p1 is lower in all coordinate
2671 // directions and check the consistency of
2672 // step sizes, i.e. that they all
2673 // add up to the sizes specified by
2674 // p_1 and p_2
2675 Point<dim> p1(p_1);
2676 Point<dim> p2(p_2);
2677 std::vector<std::vector<double>> step_sizes(step_sz);
2678
2679 for (unsigned int i = 0; i < dim; ++i)
2680 {
2681 if (p1(i) > p2(i))
2682 {
2683 std::swap(p1(i), p2(i));
2684 std::reverse(step_sizes[i].begin(), step_sizes[i].end());
2685 }
2686
2687# ifdef DEBUG
2688 double x = 0;
2689 for (unsigned int j = 0; j < step_sizes.at(i).size(); ++j)
2690 x += step_sizes[i][j];
2691 Assert(std::fabs(x - (p2(i) - p1(i))) <= 1e-12 * std::fabs(x),
2692 ExcMessage(
2693 "The sequence of step sizes in coordinate direction " +
2695 " must be equal to the distance of the two given "
2696 "points in this coordinate direction."));
2697# endif
2698 }
2699
2700
2701 // then generate the necessary
2702 // points
2703 std::vector<Point<dim>> points;
2704 switch (dim)
2705 {
2706 case 1:
2707 {
2708 double x = 0;
2709 for (unsigned int i = 0;; ++i)
2710 {
2711 points.push_back(Point<dim>(p1[0] + x));
2712
2713 // form partial sums. in
2714 // the last run through
2715 // avoid accessing
2716 // non-existent values
2717 // and exit early instead
2718 if (i == step_sizes[0].size())
2719 break;
2720
2721 x += step_sizes[0][i];
2722 }
2723 break;
2724 }
2725
2726 case 2:
2727 {
2728 double y = 0;
2729 for (unsigned int j = 0;; ++j)
2730 {
2731 double x = 0;
2732 for (unsigned int i = 0;; ++i)
2733 {
2734 points.push_back(Point<dim>(p1[0] + x, p1[1] + y));
2735 if (i == step_sizes[0].size())
2736 break;
2737
2738 x += step_sizes[0][i];
2739 }
2740
2741 if (j == step_sizes[1].size())
2742 break;
2743
2744 y += step_sizes[1][j];
2745 }
2746 break;
2747 }
2748 case 3:
2749 {
2750 double z = 0;
2751 for (unsigned int k = 0;; ++k)
2752 {
2753 double y = 0;
2754 for (unsigned int j = 0;; ++j)
2755 {
2756 double x = 0;
2757 for (unsigned int i = 0;; ++i)
2758 {
2759 points.push_back(
2760 Point<dim>(p1[0] + x, p1[1] + y, p1[2] + z));
2761 if (i == step_sizes[0].size())
2762 break;
2763
2764 x += step_sizes[0][i];
2765 }
2766
2767 if (j == step_sizes[1].size())
2768 break;
2769
2770 y += step_sizes[1][j];
2771 }
2772
2773 if (k == step_sizes[2].size())
2774 break;
2775
2776 z += step_sizes[2][k];
2777 }
2778 break;
2779 }
2780
2781 default:
2782 Assert(false, ExcNotImplemented());
2783 }
2784
2785 // next create the cells
2786 // Prepare cell data
2787 std::vector<CellData<dim>> cells;
2788 switch (dim)
2789 {
2790 case 1:
2791 {
2792 cells.resize(step_sizes[0].size());
2793 for (unsigned int x = 0; x < step_sizes[0].size(); ++x)
2794 {
2795 cells[x].vertices[0] = x;
2796 cells[x].vertices[1] = x + 1;
2797 cells[x].material_id = 0;
2798 }
2799 break;
2800 }
2801
2802 case 2:
2803 {
2804 cells.resize(step_sizes[1].size() * step_sizes[0].size());
2805 for (unsigned int y = 0; y < step_sizes[1].size(); ++y)
2806 for (unsigned int x = 0; x < step_sizes[0].size(); ++x)
2807 {
2808 const unsigned int c = x + y * step_sizes[0].size();
2809 cells[c].vertices[0] = y * (step_sizes[0].size() + 1) + x;
2810 cells[c].vertices[1] = y * (step_sizes[0].size() + 1) + x + 1;
2811 cells[c].vertices[2] =
2812 (y + 1) * (step_sizes[0].size() + 1) + x;
2813 cells[c].vertices[3] =
2814 (y + 1) * (step_sizes[0].size() + 1) + x + 1;
2815 cells[c].material_id = 0;
2816 }
2817 break;
2818 }
2819
2820 case 3:
2821 {
2822 const unsigned int n_x = (step_sizes[0].size() + 1);
2823 const unsigned int n_xy =
2824 (step_sizes[0].size() + 1) * (step_sizes[1].size() + 1);
2825
2826 cells.resize(step_sizes[2].size() * step_sizes[1].size() *
2827 step_sizes[0].size());
2828 for (unsigned int z = 0; z < step_sizes[2].size(); ++z)
2829 for (unsigned int y = 0; y < step_sizes[1].size(); ++y)
2830 for (unsigned int x = 0; x < step_sizes[0].size(); ++x)
2831 {
2832 const unsigned int c =
2833 x + y * step_sizes[0].size() +
2834 z * step_sizes[0].size() * step_sizes[1].size();
2835 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2836 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2837 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2838 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2839 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2840 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2841 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2842 cells[c].vertices[7] =
2843 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2844 cells[c].material_id = 0;
2845 }
2846 break;
2847 }
2848
2849 default:
2850 Assert(false, ExcNotImplemented());
2851 }
2852
2853 tria.create_triangulation(points, cells, SubCellData());
2854
2855 if (colorize)
2856 {
2857 // to colorize, run through all
2858 // faces of all cells and set
2859 // boundary indicator to the
2860 // correct value if it was 0.
2861
2862 // use a large epsilon to
2863 // compare numbers to avoid
2864 // roundoff problems.
2865 double min_size =
2866 *std::min_element(step_sizes[0].begin(), step_sizes[0].end());
2867 for (unsigned int i = 1; i < dim; ++i)
2868 min_size = std::min(min_size,
2869 *std::min_element(step_sizes[i].begin(),
2870 step_sizes[i].end()));
2871 const double epsilon = 0.01 * min_size;
2872
2873 // actual code is external since
2874 // 1-D is different from 2/3D.
2875 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
2876 }
2877 }
2878
2879
2880
2881 template <>
2882 void
2884 const std::vector<std::vector<double>> &spacing,
2885 const Point<1> & p,
2887 const bool colorize)
2888 {
2889 Assert(spacing.size() == 1, ExcInvalidRepetitionsDimension(1));
2890
2891 const unsigned int n_cells = material_id.size(0);
2892
2893 Assert(spacing[0].size() == n_cells, ExcInvalidRepetitionsDimension(1));
2894
2895 double delta = std::numeric_limits<double>::max();
2896 for (unsigned int i = 0; i < n_cells; ++i)
2897 {
2898 Assert(spacing[0][i] >= 0, ExcInvalidRepetitions(-1));
2899 delta = std::min(delta, spacing[0][i]);
2900 }
2901
2902 // generate the necessary points
2903 std::vector<Point<1>> points;
2904 double ax = p[0];
2905 for (unsigned int x = 0; x <= n_cells; ++x)
2906 {
2907 points.emplace_back(ax);
2908 if (x < n_cells)
2909 ax += spacing[0][x];
2910 }
2911 // create the cells
2912 unsigned int n_val_cells = 0;
2913 for (unsigned int i = 0; i < n_cells; ++i)
2915 n_val_cells++;
2916
2917 std::vector<CellData<1>> cells(n_val_cells);
2918 unsigned int id = 0;
2919 for (unsigned int x = 0; x < n_cells; ++x)
2921 {
2922 cells[id].vertices[0] = x;
2923 cells[id].vertices[1] = x + 1;
2924 cells[id].material_id = material_id[x];
2925 id++;
2926 }
2927 // create triangulation
2928 SubCellData t;
2929 GridTools::delete_unused_vertices(points, cells, t);
2930
2931 tria.create_triangulation(points, cells, t);
2932
2933 // set boundary indicator
2934 if (colorize)
2935 Assert(false, ExcNotImplemented());
2936 }
2937
2938
2939 template <>
2940 void
2942 const std::vector<std::vector<double>> &spacing,
2943 const Point<2> & p,
2945 const bool colorize)
2946 {
2947 Assert(spacing.size() == 2, ExcInvalidRepetitionsDimension(2));
2948
2949 std::vector<unsigned int> repetitions(2);
2950 double delta = std::numeric_limits<double>::max();
2951 for (unsigned int i = 0; i < 2; ++i)
2952 {
2953 repetitions[i] = spacing[i].size();
2954 for (unsigned int j = 0; j < repetitions[i]; ++j)
2955 {
2956 Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
2957 delta = std::min(delta, spacing[i][j]);
2958 }
2959 Assert(material_id.size(i) == repetitions[i],
2961 }
2962
2963 // generate the necessary points
2964 std::vector<Point<2>> points;
2965 double ay = p[1];
2966 for (unsigned int y = 0; y <= repetitions[1]; ++y)
2967 {
2968 double ax = p[0];
2969 for (unsigned int x = 0; x <= repetitions[0]; ++x)
2970 {
2971 points.emplace_back(ax, ay);
2972 if (x < repetitions[0])
2973 ax += spacing[0][x];
2974 }
2975 if (y < repetitions[1])
2976 ay += spacing[1][y];
2977 }
2978
2979 // create the cells
2980 unsigned int n_val_cells = 0;
2981 for (unsigned int i = 0; i < material_id.size(0); ++i)
2982 for (unsigned int j = 0; j < material_id.size(1); ++j)
2984 n_val_cells++;
2985
2986 std::vector<CellData<2>> cells(n_val_cells);
2987 unsigned int id = 0;
2988 for (unsigned int y = 0; y < repetitions[1]; ++y)
2989 for (unsigned int x = 0; x < repetitions[0]; ++x)
2991 {
2992 cells[id].vertices[0] = y * (repetitions[0] + 1) + x;
2993 cells[id].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2994 cells[id].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2995 cells[id].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
2996 cells[id].material_id = material_id[x][y];
2997 id++;
2998 }
2999
3000 // create triangulation
3001 SubCellData t;
3002 GridTools::delete_unused_vertices(points, cells, t);
3003
3004 tria.create_triangulation(points, cells, t);
3005
3006 // set boundary indicator
3007 if (colorize)
3008 {
3009 double eps = 0.01 * delta;
3011 for (; cell != endc; ++cell)
3012 {
3013 Point<2> cell_center = cell->center();
3014 for (unsigned int f : GeometryInfo<2>::face_indices())
3015 if (cell->face(f)->boundary_id() == 0)
3016 {
3017 Point<2> face_center = cell->face(f)->center();
3018 for (unsigned int i = 0; i < 2; ++i)
3019 {
3020 if (face_center[i] < cell_center[i] - eps)
3021 cell->face(f)->set_boundary_id(i * 2);
3022 if (face_center[i] > cell_center[i] + eps)
3023 cell->face(f)->set_boundary_id(i * 2 + 1);
3024 }
3025 }
3026 }
3027 }
3028 }
3029
3030
3031 template <>
3032 void
3034 const std::vector<std::vector<double>> &spacing,
3035 const Point<3> & p,
3037 const bool colorize)
3038 {
3039 const unsigned int dim = 3;
3040
3041 Assert(spacing.size() == dim, ExcInvalidRepetitionsDimension(dim));
3042
3043 std::array<unsigned int, dim> repetitions;
3044 double delta = std::numeric_limits<double>::max();
3045 for (unsigned int i = 0; i < dim; ++i)
3046 {
3047 repetitions[i] = spacing[i].size();
3048 for (unsigned int j = 0; j < repetitions[i]; ++j)
3049 {
3050 Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
3051 delta = std::min(delta, spacing[i][j]);
3052 }
3053 Assert(material_id.size(i) == repetitions[i],
3055 }
3056
3057 // generate the necessary points
3058 std::vector<Point<dim>> points;
3059 double az = p[2];
3060 for (unsigned int z = 0; z <= repetitions[2]; ++z)
3061 {
3062 double ay = p[1];
3063 for (unsigned int y = 0; y <= repetitions[1]; ++y)
3064 {
3065 double ax = p[0];
3066 for (unsigned int x = 0; x <= repetitions[0]; ++x)
3067 {
3068 points.emplace_back(ax, ay, az);
3069 if (x < repetitions[0])
3070 ax += spacing[0][x];
3071 }
3072 if (y < repetitions[1])
3073 ay += spacing[1][y];
3074 }
3075 if (z < repetitions[2])
3076 az += spacing[2][z];
3077 }
3078
3079 // create the cells
3080 unsigned int n_val_cells = 0;
3081 for (unsigned int i = 0; i < material_id.size(0); ++i)
3082 for (unsigned int j = 0; j < material_id.size(1); ++j)
3083 for (unsigned int k = 0; k < material_id.size(2); ++k)
3085 n_val_cells++;
3086
3087 std::vector<CellData<dim>> cells(n_val_cells);
3088 unsigned int id = 0;
3089 const unsigned int n_x = (repetitions[0] + 1);
3090 const unsigned int n_xy = (repetitions[0] + 1) * (repetitions[1] + 1);
3091 for (unsigned int z = 0; z < repetitions[2]; ++z)
3092 for (unsigned int y = 0; y < repetitions[1]; ++y)
3093 for (unsigned int x = 0; x < repetitions[0]; ++x)
3095 {
3096 cells[id].vertices[0] = z * n_xy + y * n_x + x;
3097 cells[id].vertices[1] = z * n_xy + y * n_x + x + 1;
3098 cells[id].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3099 cells[id].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3100 cells[id].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3101 cells[id].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3102 cells[id].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3103 cells[id].vertices[7] = (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3104 cells[id].material_id = material_id[x][y][z];
3105 id++;
3106 }
3107
3108 // create triangulation
3109 SubCellData t;
3110 GridTools::delete_unused_vertices(points, cells, t);
3111
3112 tria.create_triangulation(points, cells, t);
3113
3114 // set boundary indicator
3115 if (colorize)
3116 {
3117 double eps = 0.01 * delta;
3119 endc = tria.end();
3120 for (; cell != endc; ++cell)
3121 {
3122 Point<dim> cell_center = cell->center();
3123 for (auto f : GeometryInfo<dim>::face_indices())
3124 if (cell->face(f)->boundary_id() == 0)
3125 {
3126 Point<dim> face_center = cell->face(f)->center();
3127 for (unsigned int i = 0; i < dim; ++i)
3128 {
3129 if (face_center[i] < cell_center[i] - eps)
3130 cell->face(f)->set_boundary_id(i * 2);
3131 if (face_center[i] > cell_center[i] + eps)
3132 cell->face(f)->set_boundary_id(i * 2 + 1);
3133 }
3134 }
3135 }
3136 }
3137 }
3138
3139 template <int dim, int spacedim>
3140 void
3142 const std::vector<unsigned int> &holes)
3143 {
3144 AssertDimension(holes.size(), dim);
3145 // The corner points of the first cell. If there is a desire at
3146 // some point to change the geometry of the cells, they can be
3147 // made an argument to the function.
3148
3149 Point<spacedim> p1;
3150 Point<spacedim> p2;
3151 for (unsigned int d = 0; d < dim; ++d)
3152 p2(d) = 1.;
3153
3154 // then check that all repetitions
3155 // are >= 1, and calculate deltas
3156 // convert repetitions from double
3157 // to int by taking the ceiling.
3158 std::array<Point<spacedim>, dim> delta;
3159 std::array<unsigned int, dim> repetitions;
3160 for (unsigned int i = 0; i < dim; ++i)
3161 {
3162 Assert(holes[i] >= 1,
3163 ExcMessage("At least one hole needed in each direction"));
3164 repetitions[i] = 2 * holes[i] + 1;
3165 delta[i][i] = (p2[i] - p1[i]);
3166 }
3167
3168 // then generate the necessary
3169 // points
3170 std::vector<Point<spacedim>> points;
3171 switch (dim)
3172 {
3173 case 1:
3174 for (unsigned int x = 0; x <= repetitions[0]; ++x)
3175 points.push_back(p1 + x * delta[0]);
3176 break;
3177
3178 case 2:
3179 for (unsigned int y = 0; y <= repetitions[1]; ++y)
3180 for (unsigned int x = 0; x <= repetitions[0]; ++x)
3181 points.push_back(p1 + x * delta[0] + y * delta[1]);
3182 break;
3183
3184 case 3:
3185 for (unsigned int z = 0; z <= repetitions[2]; ++z)
3186 for (unsigned int y = 0; y <= repetitions[1]; ++y)
3187 for (unsigned int x = 0; x <= repetitions[0]; ++x)
3188 points.push_back(p1 + x * delta[0] + y * delta[1] +
3189 z * delta[2]);
3190 break;
3191
3192 default:
3193 Assert(false, ExcNotImplemented());
3194 }
3195
3196 // next create the cells
3197 // Prepare cell data
3198 std::vector<CellData<dim>> cells;
3199 switch (dim)
3200 {
3201 case 2:
3202 {
3203 cells.resize(repetitions[1] * repetitions[0] - holes[1] * holes[0]);
3204 unsigned int c = 0;
3205 for (unsigned int y = 0; y < repetitions[1]; ++y)
3206 for (unsigned int x = 0; x < repetitions[0]; ++x)
3207 {
3208 if ((x % 2 == 1) && (y % 2 == 1))
3209 continue;
3210 Assert(c < cells.size(), ExcInternalError());
3211 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
3212 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
3213 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
3214 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
3215 cells[c].material_id = 0;
3216 ++c;
3217 }
3218 break;
3219 }
3220
3221 case 3:
3222 {
3223 const unsigned int n_x = (repetitions[0] + 1);
3224 const unsigned int n_xy =
3225 (repetitions[0] + 1) * (repetitions[1] + 1);
3226
3227 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
3228
3229 unsigned int c = 0;
3230 for (unsigned int z = 0; z < repetitions[2]; ++z)
3231 for (unsigned int y = 0; y < repetitions[1]; ++y)
3232 for (unsigned int x = 0; x < repetitions[0]; ++x)
3233 {
3234 Assert(c < cells.size(), ExcInternalError());
3235 cells[c].vertices[0] = z * n_xy + y * n_x + x;
3236 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
3237 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3238 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3239 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3240 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3241 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3242 cells[c].vertices[7] =
3243 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3244 cells[c].material_id = 0;
3245 ++c;
3246 }
3247 break;
3248 }
3249
3250 default:
3251 Assert(false, ExcNotImplemented());
3252 }
3253
3254 tria.create_triangulation(points, cells, SubCellData());
3255 }
3256
3257
3258
3259 template <>
3260 void
3262 const double /*inner_radius*/,
3263 const double /*outer_radius*/,
3264 const double /*pad_bottom*/,
3265 const double /*pad_top*/,
3266 const double /*pad_left*/,
3267 const double /*pad_right*/,
3268 const Point<1> & /*center*/,
3269 const types::manifold_id /*polar_manifold_id*/,
3270 const types::manifold_id /*tfi_manifold_id*/,
3271 const double /*L*/,
3272 const unsigned int /*n_slices*/,
3273 const bool /*colorize*/)
3274 {
3275 Assert(false, ExcNotImplemented());
3276 }
3277
3278
3279
3280 template <>
3281 void
3283 const double /*shell_region_width*/,
3284 const unsigned int /*n_shells*/,
3285 const double /*skewness*/,
3286 const bool /*colorize*/)
3287 {
3288 Assert(false, ExcNotImplemented());
3289 }
3290
3291
3292
3293 namespace internal
3294 {
3295 // helper function to check if point is in 2D box
3296 bool inline point_in_2d_box(const Point<2> &p,
3297 const Point<2> &c,
3298 const double radius)
3299 {
3300 return (std::abs(p[0] - c[0]) < radius) &&
3301 (std::abs(p[1] - c[1]) < radius);
3302 }
3303
3304
3305
3306 // Find the minimal distance between two vertices. This is useful for
3307 // computing a tolerance for merging vertices in
3308 // GridTools::merge_triangulations.
3309 template <int dim, int spacedim>
3310 double
3311 minimal_vertex_distance(const Triangulation<dim, spacedim> &triangulation)
3312 {
3313 double length = std::numeric_limits<double>::max();
3314 for (const auto &cell : triangulation.active_cell_iterators())
3315 for (unsigned int n = 0; n < GeometryInfo<dim>::lines_per_cell; ++n)
3316 length = std::min(length, cell->line(n)->diameter());
3317 return length;
3318 }
3319 } // namespace internal
3320
3321
3322
3323 template <>
3324 void
3326 const double inner_radius,
3327 const double outer_radius,
3328 const double pad_bottom,
3329 const double pad_top,
3330 const double pad_left,
3331 const double pad_right,
3332 const Point<2> & new_center,
3333 const types::manifold_id polar_manifold_id,
3334 const types::manifold_id tfi_manifold_id,
3335 const double L,
3336 const unsigned int /*n_slices*/,
3337 const bool colorize)
3338 {
3339 const bool with_padding =
3340 pad_bottom > 0 || pad_top > 0 || pad_left > 0 || pad_right > 0;
3341
3342 Assert(pad_bottom >= 0., ExcMessage("Negative bottom padding."));
3343 Assert(pad_top >= 0., ExcMessage("Negative top padding."));
3344 Assert(pad_left >= 0., ExcMessage("Negative left padding."));
3345 Assert(pad_right >= 0., ExcMessage("Negative right padding."));
3346
3347 const Point<2> center;
3348
3349 auto min_line_length = [](const Triangulation<2> &tria) -> double {
3350 double length = std::numeric_limits<double>::max();
3351 for (const auto &cell : tria.active_cell_iterators())
3352 for (unsigned int n = 0; n < GeometryInfo<2>::lines_per_cell; ++n)
3353 length = std::min(length, cell->line(n)->diameter());
3354 return length;
3355 };
3356
3357 // start by setting up the cylinder triangulation
3358 Triangulation<2> cylinder_tria_maybe;
3359 Triangulation<2> &cylinder_tria = with_padding ? cylinder_tria_maybe : tria;
3361 inner_radius,
3362 outer_radius,
3363 L,
3364 /*repetitions*/ 1,
3365 colorize);
3366
3367 // we will deal with face manifold ids after we merge triangulations
3368 for (const auto &cell : cylinder_tria.active_cell_iterators())
3369 cell->set_manifold_id(tfi_manifold_id);
3370
3371 const Point<2> bl(-outer_radius - pad_left, -outer_radius - pad_bottom);
3372 const Point<2> tr(outer_radius + pad_right, outer_radius + pad_top);
3373 if (with_padding)
3374 {
3375 // hyper_cube_with_cylindrical_hole will have 2 cells along
3376 // each face, so the element size is outer_radius
3377
3378 auto add_sizes = [](std::vector<double> &step_sizes,
3379 const double padding,
3380 const double h) -> void {
3381 // use std::round instead of std::ceil to improve aspect ratio
3382 // in case padding is only slightly larger than h.
3383 const auto rounded =
3384 static_cast<unsigned int>(std::round(padding / h));
3385 // in case padding is much smaller than h, make sure we
3386 // have at least 1 element
3387 const unsigned int num = (padding > 0. && rounded == 0) ? 1 : rounded;
3388 for (unsigned int i = 0; i < num; ++i)
3389 step_sizes.push_back(padding / num);
3390 };
3391
3392 std::vector<std::vector<double>> step_sizes(2);
3393 // x-coord
3394 // left:
3395 add_sizes(step_sizes[0], pad_left, outer_radius);
3396 // center
3397 step_sizes[0].push_back(outer_radius);
3398 step_sizes[0].push_back(outer_radius);
3399 // right
3400 add_sizes(step_sizes[0], pad_right, outer_radius);
3401 // y-coord
3402 // bottom
3403 add_sizes(step_sizes[1], pad_bottom, outer_radius);
3404 // center
3405 step_sizes[1].push_back(outer_radius);
3406 step_sizes[1].push_back(outer_radius);
3407 // top
3408 add_sizes(step_sizes[1], pad_top, outer_radius);
3409
3410 // now create bulk
3411 Triangulation<2> bulk_tria;
3413 bulk_tria, step_sizes, bl, tr, colorize);
3414
3415 // now remove cells reserved from the cylindrical hole
3416 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3417 for (const auto &cell : bulk_tria.active_cell_iterators())
3418 if (internal::point_in_2d_box(cell->center(), center, outer_radius))
3419 cells_to_remove.insert(cell);
3420
3421 Triangulation<2> tria_without_cylinder;
3423 bulk_tria, cells_to_remove, tria_without_cylinder);
3424
3425 const double tolerance =
3426 std::min(min_line_length(tria_without_cylinder),
3427 min_line_length(cylinder_tria)) /
3428 2.0;
3429
3430 GridGenerator::merge_triangulations(tria_without_cylinder,
3431 cylinder_tria,
3432 tria,
3433 tolerance);
3434 }
3435
3436 // now set manifold ids:
3437 for (const auto &cell : tria.active_cell_iterators())
3438 {
3439 // set all non-boundary manifold ids on the cells that came from the
3440 // grid around the cylinder to the new TFI manifold id.
3441 if (cell->manifold_id() == tfi_manifold_id)
3442 {
3443 for (const unsigned int face_n : GeometryInfo<2>::face_indices())
3444 {
3445 const auto &face = cell->face(face_n);
3446 if (face->at_boundary() &&
3447 internal::point_in_2d_box(face->center(),
3448 center,
3449 outer_radius))
3450 face->set_manifold_id(polar_manifold_id);
3451 else
3452 face->set_manifold_id(tfi_manifold_id);
3453 }
3454 }
3455 else
3456 {
3457 // ensure that all other manifold ids (including the faces
3458 // opposite the cylinder) are set to the flat id
3459 cell->set_all_manifold_ids(numbers::flat_manifold_id);
3460 }
3461 }
3462
3463 static constexpr double tol =
3465 if (colorize)
3466 for (const auto &cell : tria.active_cell_iterators())
3467 for (const unsigned int face_n : GeometryInfo<2>::face_indices())
3468 {
3469 const auto face = cell->face(face_n);
3470 if (face->at_boundary())
3471 {
3472 const Point<2> center = face->center();
3473 // left side
3474 if (std::abs(center[0] - bl[0]) < tol * std::abs(bl[0]))
3475 face->set_boundary_id(0);
3476 // right side
3477 else if (std::abs(center[0] - tr[0]) < tol * std::abs(tr[0]))
3478 face->set_boundary_id(1);
3479 // bottom
3480 else if (std::abs(center[1] - bl[1]) < tol * std::abs(bl[1]))
3481 face->set_boundary_id(2);
3482 // top
3483 else if (std::abs(center[1] - tr[1]) < tol * std::abs(tr[1]))
3484 face->set_boundary_id(3);
3485 // cylinder boundary
3486 else
3487 {
3488 Assert(cell->manifold_id() == tfi_manifold_id,
3490 face->set_boundary_id(4);
3491 }
3492 }
3493 }
3494
3495 // move to the new center
3496 GridTools::shift(new_center, tria);
3497
3498 PolarManifold<2> polar_manifold(new_center);
3499 tria.set_manifold(polar_manifold_id, polar_manifold);
3501 inner_manifold.initialize(tria);
3502 tria.set_manifold(tfi_manifold_id, inner_manifold);
3503 }
3504
3505
3506
3507 template <>
3508 void
3510 const double inner_radius,
3511 const double outer_radius,
3512 const double pad_bottom,
3513 const double pad_top,
3514 const double pad_left,
3515 const double pad_right,
3516 const Point<3> & new_center,
3517 const types::manifold_id polar_manifold_id,
3518 const types::manifold_id tfi_manifold_id,
3519 const double L,
3520 const unsigned int n_slices,
3521 const bool colorize)
3522 {
3523 Triangulation<2> tria_2;
3524 plate_with_a_hole(tria_2,
3525 inner_radius,
3526 outer_radius,
3527 pad_bottom,
3528 pad_top,
3529 pad_left,
3530 pad_right,
3531 Point<2>(new_center[0], new_center[1]),
3532 polar_manifold_id,
3533 tfi_manifold_id,
3534 L,
3535 n_slices,
3536 colorize);
3537
3538 // extrude to 3D
3539 extrude_triangulation(tria_2, n_slices, L, tria, true);
3540
3541 // shift in Z direction to match specified center
3542 GridTools::shift(Point<3>(0, 0, new_center[2] - L / 2.), tria);
3543
3544 // set up the new manifolds
3545 const Tensor<1, 3> direction{{0.0, 0.0, 1.0}};
3546 const CylindricalManifold<3> cylindrical_manifold(
3547 direction,
3548 /*axial_point*/ new_center);
3550 inner_manifold.initialize(tria);
3551 tria.set_manifold(polar_manifold_id, cylindrical_manifold);
3552 tria.set_manifold(tfi_manifold_id, inner_manifold);
3553 }
3554
3555
3556
3557 template <>
3558 void
3560 const double shell_region_width,
3561 const unsigned int n_shells,
3562 const double skewness,
3563 const bool colorize)
3564 {
3565 Assert(0.0 <= shell_region_width && shell_region_width < 0.05,
3566 ExcMessage("The width of the shell region must be less than 0.05 "
3567 "(and preferably close to 0.03)"));
3568 const types::manifold_id polar_manifold_id = 0;
3569 const types::manifold_id tfi_manifold_id = 1;
3570
3571 // We begin by setting up a grid that is 4 by 22 cells. While not
3572 // squares, these have pretty good aspect ratios.
3573 Triangulation<2> bulk_tria;
3575 {22u, 4u},
3576 Point<2>(0.0, 0.0),
3577 Point<2>(2.2, 0.41));
3578 // bulk_tria now looks like this:
3579 //
3580 // +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
3581 // | | | | | | | | | | | | | | | | | | | | | | |
3582 // +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
3583 // | |XX|XX| | | | | | | | | | | | | | | | | | | |
3584 // +--+--O--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
3585 // | |XX|XX| | | | | | | | | | | | | | | | | | | |
3586 // +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
3587 // | | | | | | | | | | | | | | | | | | | | | | |
3588 // +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
3589 //
3590 // Note that these cells are not quite squares: they are all 0.1 by
3591 // 0.1025.
3592 //
3593 // The next step is to remove the cells marked with XXs: we will place
3594 // the grid around the cylinder there later. The next loop does two
3595 // things:
3596 // 1. Determines which cells need to be removed from the Triangulation
3597 // (i.e., find the cells marked with XX in the picture).
3598 // 2. Finds the location of the vertex marked with 'O' and uses that to
3599 // calculate the shift vector for aligning cylinder_tria with
3600 // tria_without_cylinder.
3601 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3602 Tensor<1, 2> cylinder_triangulation_offset;
3603 for (const auto &cell : bulk_tria.active_cell_iterators())
3604 {
3605 if ((cell->center() - Point<2>(0.2, 0.2)).norm() < 0.15)
3606 cells_to_remove.insert(cell);
3607
3608 if (cylinder_triangulation_offset == Tensor<1, 2>())
3609 {
3610 for (const unsigned int vertex_n :
3612 if (cell->vertex(vertex_n) == Point<2>())
3613 {
3614 // cylinder_tria is centered at zero, so we need to
3615 // shift it up and to the right by two cells:
3616 cylinder_triangulation_offset =
3617 2.0 * (cell->vertex(3) - Point<2>());
3618 break;
3619 }
3620 }
3621 }
3622 Triangulation<2> tria_without_cylinder;
3624 bulk_tria, cells_to_remove, tria_without_cylinder);
3625
3626 // set up the cylinder triangulation. Note that this function sets the
3627 // manifold ids of the interior boundary cells to 0
3628 // (polar_manifold_id).
3629 Triangulation<2> cylinder_tria;
3631 0.05 + shell_region_width,
3632 0.41 / 4.0);
3633 // The bulk cells are not quite squares, so we need to move the left
3634 // and right sides of cylinder_tria inwards so that it fits in
3635 // bulk_tria:
3636 for (const auto &cell : cylinder_tria.active_cell_iterators())
3637 for (const unsigned int vertex_n : GeometryInfo<2>::vertex_indices())
3638 {
3639 if (std::abs(cell->vertex(vertex_n)[0] - -0.41 / 4.0) < 1e-10)
3640 cell->vertex(vertex_n)[0] = -0.1;
3641 else if (std::abs(cell->vertex(vertex_n)[0] - 0.41 / 4.0) < 1e-10)
3642 cell->vertex(vertex_n)[0] = 0.1;
3643 }
3644
3645 // Assign interior manifold ids to be the TFI id.
3646 for (const auto &cell : cylinder_tria.active_cell_iterators())
3647 {
3648 cell->set_manifold_id(tfi_manifold_id);
3649 for (const unsigned int face_n : GeometryInfo<2>::face_indices())
3650 if (!cell->face(face_n)->at_boundary())
3651 cell->face(face_n)->set_manifold_id(tfi_manifold_id);
3652 }
3653 if (0.0 < shell_region_width)
3654 {
3655 Assert(0 < n_shells,
3656 ExcMessage("If the shell region has positive width then "
3657 "there must be at least one shell."));
3658 Triangulation<2> shell_tria;
3660 Point<2>(),
3661 0.05,
3662 0.05 + shell_region_width,
3663 n_shells,
3664 skewness,
3665 8);
3666
3667 // Make the tolerance as large as possible since these cells can
3668 // be quite close together
3669 const double vertex_tolerance =
3670 std::min(internal::minimal_vertex_distance(shell_tria),
3671 internal::minimal_vertex_distance(cylinder_tria)) *
3672 0.5;
3673
3674 shell_tria.set_all_manifold_ids(polar_manifold_id);
3675 Triangulation<2> temp;
3677 shell_tria, cylinder_tria, temp, vertex_tolerance, true);
3678 cylinder_tria = std::move(temp);
3679 }
3680 GridTools::shift(cylinder_triangulation_offset, cylinder_tria);
3681
3682 // Compute the tolerance again, since the shells may be very close to
3683 // each-other:
3684 const double vertex_tolerance =
3685 std::min(internal::minimal_vertex_distance(tria_without_cylinder),
3686 internal::minimal_vertex_distance(cylinder_tria)) /
3687 10;
3689 tria_without_cylinder, cylinder_tria, tria, vertex_tolerance, true);
3690
3691 // Move the vertices in the middle of the faces of cylinder_tria slightly
3692 // to give a better mesh quality. We have to balance the quality of these
3693 // cells with the quality of the outer cells (initially rectangles). For
3694 // constant radial distance, we would place them at the distance 0.1 *
3695 // sqrt(2.) from the center. In case the shell region width is more than
3696 // 0.1/6., we choose to place them at 0.1 * 4./3. from the center, which
3697 // ensures that the shortest edge of the outer cells is 2./3. of the
3698 // original length. If the shell region width is less, we make the edge
3699 // length of the inner part and outer part (in the shorter x direction)
3700 // the same.
3701 {
3702 const double shift =
3703 std::min(0.125 + shell_region_width * 0.5, 0.1 * 4. / 3.);
3704 for (const auto &cell : tria.active_cell_iterators())
3705 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
3706 if (cell->vertex(v).distance(Point<2>(0.1, 0.205)) < 1e-10)
3707 cell->vertex(v) = Point<2>(0.2 - shift, 0.205);
3708 else if (cell->vertex(v).distance(Point<2>(0.3, 0.205)) < 1e-10)
3709 cell->vertex(v) = Point<2>(0.2 + shift, 0.205);
3710 else if (cell->vertex(v).distance(Point<2>(0.2, 0.1025)) < 1e-10)
3711 cell->vertex(v) = Point<2>(0.2, 0.2 - shift);
3712 else if (cell->vertex(v).distance(Point<2>(0.2, 0.3075)) < 1e-10)
3713 cell->vertex(v) = Point<2>(0.2, 0.2 + shift);
3714 }
3715
3716 // Ensure that all manifold ids on a polar cell really are set to the
3717 // polar manifold id:
3718 for (const auto &cell : tria.active_cell_iterators())
3719 if (cell->manifold_id() == polar_manifold_id)
3720 cell->set_all_manifold_ids(polar_manifold_id);
3721
3722 // Ensure that all other manifold ids (including the interior faces
3723 // opposite the cylinder) are set to the flat manifold id:
3724 for (const auto &cell : tria.active_cell_iterators())
3725 if (cell->manifold_id() != polar_manifold_id &&
3726 cell->manifold_id() != tfi_manifold_id)
3727 cell->set_all_manifold_ids(numbers::flat_manifold_id);
3728
3729 // We need to calculate the current center so that we can move it later:
3730 // to start get a unique list of (points to) vertices on the cylinder
3731 std::vector<Point<2> *> cylinder_pointers;
3732 for (const auto &face : tria.active_face_iterators())
3733 if (face->manifold_id() == polar_manifold_id)
3734 {
3735 cylinder_pointers.push_back(&face->vertex(0));
3736 cylinder_pointers.push_back(&face->vertex(1));
3737 }
3738 // de-duplicate
3739 std::sort(cylinder_pointers.begin(), cylinder_pointers.end());
3740 cylinder_pointers.erase(std::unique(cylinder_pointers.begin(),
3741 cylinder_pointers.end()),
3742 cylinder_pointers.end());
3743
3744 // find the current center...
3746 for (const Point<2> *const ptr : cylinder_pointers)
3747 center += *ptr / double(cylinder_pointers.size());
3748
3749 // and recenter at (0.2, 0.2)
3750 for (Point<2> *const ptr : cylinder_pointers)
3751 *ptr += Point<2>(0.2, 0.2) - center;
3752
3753 // attach manifolds
3754 PolarManifold<2> polar_manifold(Point<2>(0.2, 0.2));
3755 tria.set_manifold(polar_manifold_id, polar_manifold);
3757 inner_manifold.initialize(tria);
3758 tria.set_manifold(tfi_manifold_id, inner_manifold);
3759
3760 if (colorize)
3761 for (const auto &face : tria.active_face_iterators())
3762 if (face->at_boundary())
3763 {
3764 const Point<2> center = face->center();
3765 // left side
3766 if (std::abs(center[0] - 0.0) < 1e-10)
3767 face->set_boundary_id(0);
3768 // right side
3769 else if (std::abs(center[0] - 2.2) < 1e-10)
3770 face->set_boundary_id(1);
3771 // cylinder boundary
3772 else if (face->manifold_id() == polar_manifold_id)
3773 face->set_boundary_id(2);
3774 // sides of channel
3775 else
3776 {
3777 Assert(std::abs(center[1] - 0.00) < 1.0e-10 ||
3778 std::abs(center[1] - 0.41) < 1.0e-10,
3780 face->set_boundary_id(3);
3781 }
3782 }
3783 }
3784
3785
3786
3787 template <>
3788 void
3790 const double shell_region_width,
3791 const unsigned int n_shells,
3792 const double skewness,
3793 const bool colorize)
3794 {
3795 Triangulation<2> tria_2;
3797 tria_2, shell_region_width, n_shells, skewness, colorize);
3798 extrude_triangulation(tria_2, 5, 0.41, tria, true);
3799
3800 // set up the new 3D manifolds
3801 const types::manifold_id cylindrical_manifold_id = 0;
3802 const types::manifold_id tfi_manifold_id = 1;
3803 const PolarManifold<2> *const m_ptr =
3804 dynamic_cast<const PolarManifold<2> *>(
3805 &tria_2.get_manifold(cylindrical_manifold_id));
3806 Assert(m_ptr != nullptr, ExcInternalError());
3807 const Point<3> axial_point(m_ptr->center[0], m_ptr->center[1], 0.0);
3808 const Tensor<1, 3> direction{{0.0, 0.0, 1.0}};
3809
3810 const CylindricalManifold<3> cylindrical_manifold(direction, axial_point);
3812 inner_manifold.initialize(tria);
3813 tria.set_manifold(cylindrical_manifold_id, cylindrical_manifold);
3814 tria.set_manifold(tfi_manifold_id, inner_manifold);
3815
3816 // From extrude_triangulation: since the maximum boundary id of tria_2 was
3817 // 3, the bottom boundary id is 4 and the top is 5: both are walls, so set
3818 // them to 3
3819 if (colorize)
3820 for (const auto &face : tria.active_face_iterators())
3821 if (face->boundary_id() == 4 || face->boundary_id() == 5)
3822 face->set_boundary_id(3);
3823 }
3824
3825
3826
3827 template <int dim, int spacedim>
3828 void
3830 const std::vector<unsigned int> &sizes,
3831 const bool colorize)
3832 {
3834 Assert(dim > 1, ExcNotImplemented());
3835 Assert(dim < 4, ExcNotImplemented());
3836
3837 // If there is a desire at some point to change the geometry of
3838 // the cells, this tensor can be made an argument to the function.
3839 Tensor<1, dim> dimensions;
3840 for (unsigned int d = 0; d < dim; ++d)
3841 dimensions[d] = 1.;
3842
3843 std::vector<Point<spacedim>> points;
3844 unsigned int n_cells = 1;
3845 for (unsigned int i : GeometryInfo<dim>::face_indices())
3846 n_cells += sizes[i];
3847
3848 std::vector<CellData<dim>> cells(n_cells);
3849 // Vertices of the center cell
3850 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
3851 {
3853 for (unsigned int d = 0; d < dim; ++d)
3854 p(d) = 0.5 * dimensions[d] *
3857 points.push_back(p);
3858 cells[0].vertices[i] = i;
3859 }
3860 cells[0].material_id = 0;
3861
3862 // The index of the first cell of the leg.
3863 unsigned int cell_index = 1;
3864 // The legs of the cross
3865 for (const unsigned int face : GeometryInfo<dim>::face_indices())
3866 {
3867 const unsigned int oface = GeometryInfo<dim>::opposite_face[face];
3868 const unsigned int dir = GeometryInfo<dim>::unit_normal_direction[face];
3869
3870 // We are moving in the direction of face
3871 for (unsigned int j = 0; j < sizes[face]; ++j, ++cell_index)
3872 {
3873 const unsigned int last_cell = (j == 0) ? 0U : (cell_index - 1);
3874
3875 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
3876 ++v)
3877 {
3878 const unsigned int cellv =
3880 const unsigned int ocellv =
3882 // First the vertices which already exist
3883 cells[cell_index].vertices[ocellv] =
3884 cells[last_cell].vertices[cellv];
3885
3886 // Now the new vertices
3887 cells[cell_index].vertices[cellv] = points.size();
3888
3889 Point<spacedim> p = points[cells[cell_index].vertices[ocellv]];
3891 dimensions[dir];
3892 points.push_back(p);
3893 }
3894 cells[cell_index].material_id = (colorize) ? (face + 1U) : 0U;
3895 }
3896 }
3897 tria.create_triangulation(points, cells, SubCellData());
3898 }
3899
3900
3901 template <>
3902 void
3903 hyper_cube_slit(Triangulation<1> &, const double, const double, const bool)
3904 {
3905 Assert(false, ExcNotImplemented());
3906 }
3907
3908
3909
3910 template <>
3911 void
3913 const double,
3914 const double,
3915 const double,
3916 const bool)
3917 {
3918 Assert(false, ExcNotImplemented());
3919 }
3920
3921
3922
3923 template <>
3924 void
3925 hyper_L(Triangulation<1> &, const double, const double, const bool)
3926 {
3927 Assert(false, ExcNotImplemented());
3928 }
3929
3930
3931
3932 template <>
3933 void
3934 hyper_ball(Triangulation<1> &, const Point<1> &, const double, const bool)
3935 {
3936 Assert(false, ExcNotImplemented());
3937 }
3938
3939
3940
3941 template <>
3942 void
3943 hyper_ball_balanced(Triangulation<1> &, const Point<1> &, const double)
3944 {
3945 Assert(false, ExcNotImplemented());
3946 }
3947
3948
3949
3950 template <>
3951 void
3952 cylinder(Triangulation<1> &, const double, const double)
3953 {
3954 Assert(false, ExcNotImplemented());
3955 }
3956
3957
3958 template <>
3959 void
3961 const unsigned int,
3962 const double,
3963 const double)
3964 {
3965 Assert(false, ExcNotImplemented());
3966 }
3967
3968
3969
3970 template <>
3971 void
3972 truncated_cone(Triangulation<1> &, const double, const double, const double)
3973 {
3974 Assert(false, ExcNotImplemented());
3975 }
3976
3977
3978
3979 template <>
3980 void
3982 const Point<1> &,
3983 const double,
3984 const double,
3985 const unsigned int,
3986 const bool)
3987 {
3988 Assert(false, ExcNotImplemented());
3989 }
3990
3991 template <>
3992 void
3994 const double,
3995 const double,
3996 const double,
3997 const unsigned int,
3998 const unsigned int)
3999 {
4000 Assert(false, ExcNotImplemented());
4001 }
4002
4003
4004 template <>
4005 void
4006 quarter_hyper_ball(Triangulation<1> &, const Point<1> &, const double)
4007 {
4008 Assert(false, ExcNotImplemented());
4009 }
4010
4011
4012 template <>
4013 void
4014 half_hyper_ball(Triangulation<1> &, const Point<1> &, const double)
4015 {
4016 Assert(false, ExcNotImplemented());
4017 }
4018
4019
4020 template <>
4021 void
4023 const Point<1> &,
4024 const double,
4025 const double,
4026 const unsigned int,
4027 const bool)
4028 {
4029 Assert(false, ExcNotImplemented());
4030 }
4031
4032 template <>
4033 void
4035 const Point<1> &,
4036 const double,
4037 const double,
4038 const unsigned int,
4039 const bool)
4040 {
4041 Assert(false, ExcNotImplemented());
4042 }
4043
4044 template <>
4045 void
4047 const double left,
4048 const double right,
4049 const double thickness,
4050 const bool colorize)
4051 {
4052 Assert(left < right,
4053 ExcMessage("Invalid left-to-right bounds of enclosed hypercube"));
4054
4055 std::vector<Point<2>> vertices(16);
4056 double coords[4];
4057 coords[0] = left - thickness;
4058 coords[1] = left;
4059 coords[2] = right;
4060 coords[3] = right + thickness;
4061
4062 unsigned int k = 0;
4063 for (const double y : coords)
4064 for (const double x : coords)
4065 vertices[k++] = Point<2>(x, y);
4066
4067 const types::material_id materials[9] = {5, 4, 6, 1, 0, 2, 9, 8, 10};
4068
4069 std::vector<CellData<2>> cells(9);
4070 k = 0;
4071 for (unsigned int i0 = 0; i0 < 3; ++i0)
4072 for (unsigned int i1 = 0; i1 < 3; ++i1)
4073 {
4074 cells[k].vertices[0] = i1 + 4 * i0;
4075 cells[k].vertices[1] = i1 + 4 * i0 + 1;
4076 cells[k].vertices[2] = i1 + 4 * i0 + 4;
4077 cells[k].vertices[3] = i1 + 4 * i0 + 5;
4078 if (colorize)
4079 cells[k].material_id = materials[k];
4080 ++k;
4081 }
4083 cells,
4084 SubCellData()); // no boundary information
4085 }
4086
4087
4088
4089 // Implementation for 2D only
4090 template <>
4091 void
4093 const double left,
4094 const double right,
4095 const bool colorize)
4096 {
4097 const double rl2 = (right + left) / 2;
4098 const Point<2> vertices[10] = {Point<2>(left, left),
4099 Point<2>(rl2, left),
4100 Point<2>(rl2, rl2),
4101 Point<2>(left, rl2),
4102 Point<2>(right, left),
4103 Point<2>(right, rl2),
4104 Point<2>(rl2, right),
4105 Point<2>(left, right),
4106 Point<2>(right, right),
4107 Point<2>(rl2, left)};
4108 const int cell_vertices[4][4] = {{0, 1, 3, 2},
4109 {9, 4, 2, 5},
4110 {3, 2, 7, 6},
4111 {2, 5, 6, 8}};
4112 std::vector<CellData<2>> cells(4, CellData<2>());
4113 for (unsigned int i = 0; i < 4; ++i)
4114 {
4115 for (unsigned int j = 0; j < 4; ++j)
4116 cells[i].vertices[j] = cell_vertices[i][j];
4117 cells[i].material_id = 0;
4118 }
4121 cells,
4122 SubCellData()); // no boundary information
4123
4124 if (colorize)
4125 {
4127 cell->face(1)->set_boundary_id(1);
4128 ++cell;
4129 cell->face(0)->set_boundary_id(2);
4130 }
4131 }
4132
4133
4134
4135 template <>
4136 void
4138 const double radius_0,
4139 const double radius_1,
4140 const double half_length)
4141 {
4142 Point<2> vertices_tmp[4];
4143
4144 vertices_tmp[0] = Point<2>(-half_length, -radius_0);
4145 vertices_tmp[1] = Point<2>(half_length, -radius_1);
4146 vertices_tmp[2] = Point<2>(-half_length, radius_0);
4147 vertices_tmp[3] = Point<2>(half_length, radius_1);
4148
4149 const std::vector<Point<2>> vertices(std::begin(vertices_tmp),
4150 std::end(vertices_tmp));
4151 unsigned int cell_vertices[1][GeometryInfo<2>::vertices_per_cell];
4152
4153 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
4154 cell_vertices[0][i] = i;
4155
4156 std::vector<CellData<2>> cells(1, CellData<2>());
4157
4158 for (const unsigned int i : GeometryInfo<2>::vertex_indices())
4159 cells[0].vertices[i] = cell_vertices[0][i];
4160
4161 cells[0].material_id = 0;
4162 triangulation.create_triangulation(vertices, cells, SubCellData());
4163
4165
4166 cell->face(0)->set_boundary_id(1);
4167 cell->face(1)->set_boundary_id(2);
4168
4169 for (unsigned int i = 2; i < 4; ++i)
4170 cell->face(i)->set_boundary_id(0);
4171 }
4172
4173
4174
4175 // Implementation for 2D only
4176 template <>
4177 void
4179 const double a,
4180 const double b,
4181 const bool colorize)
4182 {
4183 const Point<2> vertices[8] = {Point<2>(a, a),
4184 Point<2>((a + b) / 2, a),
4185 Point<2>(b, a),
4186 Point<2>(a, (a + b) / 2),
4187 Point<2>((a + b) / 2, (a + b) / 2),
4188 Point<2>(b, (a + b) / 2),
4189 Point<2>(a, b),
4190 Point<2>((a + b) / 2, b)};
4191 const int cell_vertices[3][4] = {{0, 1, 3, 4}, {1, 2, 4, 5}, {3, 4, 6, 7}};
4192
4193 std::vector<CellData<2>> cells(3, CellData<2>());
4194
4195 for (unsigned int i = 0; i < 3; ++i)
4196 {
4197 for (unsigned int j = 0; j < 4; ++j)
4198 cells[i].vertices[j] = cell_vertices[i][j];
4199 cells[i].material_id = 0;
4200 }
4201
4204 cells,
4205 SubCellData());
4206
4207 if (colorize)
4208 {
4210
4211 cell->face(0)->set_boundary_id(0);
4212 cell->face(2)->set_boundary_id(1);
4213 cell++;
4214
4215 cell->face(1)->set_boundary_id(2);
4216 cell->face(2)->set_boundary_id(1);
4217 cell->face(3)->set_boundary_id(3);
4218 cell++;
4219
4220 cell->face(0)->set_boundary_id(0);
4221 cell->face(1)->set_boundary_id(4);
4222 cell->face(3)->set_boundary_id(5);
4223 }
4224 }
4225
4226
4227
4228 template <int dim, int spacedim>
4229 void
4231 const std::vector<unsigned int> &repetitions,
4232 const Point<dim> & bottom_left,
4233 const Point<dim> & top_right,
4234 const std::vector<int> & n_cells_to_remove)
4235 {
4236 Assert(dim > 1, ExcNotImplemented());
4237 // Check the consistency of the dimensions provided.
4238 AssertDimension(repetitions.size(), dim);
4239 AssertDimension(n_cells_to_remove.size(), dim);
4240 for (unsigned int d = 0; d < dim; ++d)
4241 {
4242 Assert(std::fabs(n_cells_to_remove[d]) <= repetitions[d],
4243 ExcMessage("Attempting to cut away too many cells."));
4244 }
4245 // Create the domain to be cut
4248 repetitions,
4249 bottom_left,
4250 top_right);
4251 // compute the vertex of the cut step, we will cut according to the
4252 // location of the cartesian coordinates of the cell centers
4253 std::array<double, dim> h;
4254 Point<dim> cut_step;
4255 for (unsigned int d = 0; d < dim; ++d)
4256 {
4257 // mesh spacing in each direction in cartesian coordinates
4258 h[d] = (top_right[d] - bottom_left[d]) / repetitions[d];
4259 // left to right, bottom to top, front to back
4260 if (n_cells_to_remove[d] >= 0)
4261 {
4262 // cartesian coordinates of vertex location
4263 cut_step[d] =
4264 h[d] * std::fabs(n_cells_to_remove[d]) + bottom_left[d];
4265 }
4266 // right to left, top to bottom, back to front
4267 else
4268 {
4269 cut_step[d] = top_right[d] - h[d] * std::fabs(n_cells_to_remove[d]);
4270 }
4271 }
4272
4273
4274 // compute cells to remove
4275 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>
4276 cells_to_remove;
4277 for (const auto &cell : rectangle.active_cell_iterators())
4278 {
4279 bool remove_cell = true;
4280 for (unsigned int d = 0; d < dim && remove_cell; ++d)
4281 if ((n_cells_to_remove[d] > 0 && cell->center()[d] >= cut_step[d]) ||
4282 (n_cells_to_remove[d] < 0 && cell->center()[d] <= cut_step[d]))
4283 remove_cell = false;
4284 if (remove_cell)
4285 cells_to_remove.insert(cell);
4286 }
4287
4289 cells_to_remove,
4290 tria);
4291 }
4292
4293
4294
4295 // Implementation for 2D only
4296 template <>
4297 void
4299 const Point<2> & p,
4300 const double radius,
4301 const bool internal_manifolds)
4302 {
4303 // equilibrate cell sizes at
4304 // transition from the inner part
4305 // to the radial cells
4306 const double a = 1. / (1 + std::sqrt(2.0));
4307 const Point<2> vertices[8] = {
4308 p + Point<2>(-1, -1) * (radius / std::sqrt(2.0)),
4309 p + Point<2>(+1, -1) * (radius / std::sqrt(2.0)),
4310 p + Point<2>(-1, -1) * (radius / std::sqrt(2.0) * a),
4311 p + Point<2>(+1, -1) * (radius / std::sqrt(2.0) * a),
4312 p + Point<2>(-1, +1) * (radius / std::sqrt(2.0) * a),
4313 p + Point<2>(+1, +1) * (radius / std::sqrt(2.0) * a),
4314 p + Point<2>(-1, +1) * (radius / std::sqrt(2.0)),
4315 p + Point<2>(+1, +1) * (radius / std::sqrt(2.0))};
4316
4317 const int cell_vertices[5][4] = {
4318 {0, 1, 2, 3}, {0, 2, 6, 4}, {2, 3, 4, 5}, {1, 7, 3, 5}, {6, 4, 7, 5}};
4319
4320 std::vector<CellData<2>> cells(5, CellData<2>());
4321
4322 for (unsigned int i = 0; i < 5; ++i)
4323 {
4324 for (unsigned int j = 0; j < 4; ++j)
4325 cells[i].vertices[j] = cell_vertices[i][j];
4326 cells[i].material_id = 0;
4327 cells[i].manifold_id = i == 2 ? numbers::flat_manifold_id : 1;
4328 }
4329
4332 cells,
4333 SubCellData()); // no boundary information
4336 if (internal_manifolds)
4338 }
4339
4340
4341
4342 template <>
4343 void
4345 const Point<2> & center,
4346 const double inner_radius,
4347 const double outer_radius,
4348 const unsigned int n_cells,
4349 const bool colorize)
4350 {
4351 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4352 ExcInvalidRadii());
4353
4354 const double pi = numbers::PI;
4355
4356 // determine the number of cells
4357 // for the grid. if not provided by
4358 // the user determine it such that
4359 // the length of each cell on the
4360 // median (in the middle between
4361 // the two circles) is equal to its
4362 // radial extent (which is the
4363 // difference between the two
4364 // radii)
4365 const unsigned int N =
4366 (n_cells == 0 ? static_cast<unsigned int>(
4367 std::ceil((2 * pi * (outer_radius + inner_radius) / 2) /
4368 (outer_radius - inner_radius))) :
4369 n_cells);
4370
4371 // set up N vertices on the
4372 // outer and N vertices on
4373 // the inner circle. the
4374 // first N ones are on the
4375 // outer one, and all are
4376 // numbered counter-clockwise
4377 std::vector<Point<2>> vertices(2 * N);
4378 for (unsigned int i = 0; i < N; ++i)
4379 {
4380 vertices[i] =
4381 Point<2>(std::cos(2 * pi * i / N), std::sin(2 * pi * i / N)) *
4382 outer_radius;
4383 vertices[i + N] = vertices[i] * (inner_radius / outer_radius);
4384
4385 vertices[i] += center;
4386 vertices[i + N] += center;
4387 }
4388
4389 std::vector<CellData<2>> cells(N, CellData<2>());
4390
4391 for (unsigned int i = 0; i < N; ++i)
4392 {
4393 cells[i].vertices[0] = i;
4394 cells[i].vertices[1] = (i + 1) % N;
4395 cells[i].vertices[2] = N + i;
4396 cells[i].vertices[3] = N + ((i + 1) % N);
4397
4398 cells[i].material_id = 0;
4399 }
4400
4402
4403 if (colorize)
4404 colorize_hyper_shell(tria, center, inner_radius, outer_radius);
4405
4408 }
4409
4410
4411
4412 template <int dim>
4413 void
4415 const Point<dim> & inner_center,
4416 const Point<dim> & outer_center,
4417 const double inner_radius,
4418 const double outer_radius,
4419 const unsigned int n_cells)
4420 {
4422 tria, outer_center, inner_radius, outer_radius, n_cells, true);
4423
4424 // check the consistency of the dimensions provided
4425 Assert(
4426 outer_radius - inner_radius > outer_center.distance(inner_center),
4428 "The inner radius is greater than or equal to the outer radius plus eccentricity."));
4429
4430 // shift nodes along the inner boundary according to the position of
4431 // inner_circle
4432 std::set<Point<dim> *> vertices_to_move;
4433
4434 for (const auto &face : tria.active_face_iterators())
4435 if (face->boundary_id() == 0)
4436 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
4437 vertices_to_move.insert(&face->vertex(v));
4438
4439 const auto shift = inner_center - outer_center;
4440 for (const auto &p : vertices_to_move)
4441 (*p) += shift;
4442
4443 // the original hyper_shell function assigns the same manifold id
4444 // to all cells and faces. Set all manifolds ids to a different
4445 // value (2), then use boundary ids to assign different manifolds to
4446 // the inner (0) and outer manifolds (1). Use a transfinite manifold
4447 // for all faces and cells aside from the boundaries.
4450
4451 SphericalManifold<dim> inner_manifold(inner_center);
4452 SphericalManifold<dim> outer_manifold(outer_center);
4453
4455 transfinite.initialize(tria);
4456
4457 tria.set_manifold(0, inner_manifold);
4458 tria.set_manifold(1, outer_manifold);
4459 tria.set_manifold(2, transfinite);
4460 }
4461
4462
4463
4464 // Implementation for 2D only
4465 template <>
4466 void
4468 const double radius,
4469 const double half_length)
4470 {
4471 Point<2> p1(-half_length, -radius);
4472 Point<2> p2(half_length, radius);
4473
4474 hyper_rectangle(tria, p1, p2, true);
4475
4478 while (f != end)
4479 {
4480 switch (f->boundary_id())
4481 {
4482 case 0:
4483 f->set_boundary_id(1);
4484 break;
4485 case 1:
4486 f->set_boundary_id(2);
4487 break;
4488 default:
4489 f->set_boundary_id(0);
4490 break;
4491 }
4492 ++f;
4493 }
4494 }
4495
4496 template <>
4497 void
4499 const unsigned int,
4500 const double,
4501 const double)
4502 {
4503 Assert(false, ExcNotImplemented());
4504 }
4505
4506
4507
4508 // Implementation for 2D only
4509 template <>
4510 void
4512 const double,
4513 const double,
4514 const double,
4515 const unsigned int,
4516 const unsigned int)
4517 {
4518 Assert(false, ExcNotImplemented());
4519 }
4520
4521
4522 template <>
4523 void
4525 const Point<2> & p,
4526 const double radius)
4527 {
4528 const unsigned int dim = 2;
4529
4530 // the numbers 0.55647 and 0.42883 have been found by a search for the
4531 // best aspect ratio (defined as the maximal between the minimal singular
4532 // value of the Jacobian)
4533 const Point<dim> vertices[7] = {p + Point<dim>(0, 0) * radius,
4534 p + Point<dim>(+1, 0) * radius,
4535 p + Point<dim>(+1, 0) * (radius * 0.55647),
4536 p + Point<dim>(0, +1) * (radius * 0.55647),
4537 p + Point<dim>(+1, +1) * (radius * 0.42883),
4538 p + Point<dim>(0, +1) * radius,
4539 p + Point<dim>(+1, +1) *
4540 (radius / std::sqrt(2.0))};
4541
4542 const int cell_vertices[3][4] = {{0, 2, 3, 4}, {1, 6, 2, 4}, {5, 3, 6, 4}};
4543
4544 std::vector<CellData<dim>> cells(3, CellData<dim>());
4545
4546 for (unsigned int i = 0; i < 3; ++i)
4547 {
4548 for (unsigned int j = 0; j < 4; ++j)
4549 cells[i].vertices[j] = cell_vertices[i][j];
4550 cells[i].material_id = 0;
4551 }
4552
4555 cells,
4556 SubCellData()); // no boundary information
4557
4560
4562
4563 while (cell != end)
4564 {
4565 for (unsigned int i : GeometryInfo<dim>::face_indices())
4566 {
4567 if (cell->face(i)->boundary_id() ==
4569 continue;
4570
4571 // If one the components is the same as the respective
4572 // component of the center, then this is part of the plane
4573 if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius ||
4574 cell->face(i)->center()(1) < p(1) + 1.e-5 * radius)
4575 {
4576 cell->face(i)->set_boundary_id(1);
4577 cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
4578 }
4579 }
4580 ++cell;
4581 }
4583 }
4584
4585
4586 template <>
4587 void
4589 const Point<2> & p,
4590 const double radius)
4591 {
4592 // equilibrate cell sizes at
4593 // transition from the inner part
4594 // to the radial cells
4595 const double a = 1. / (1 + std::sqrt(2.0));
4596 const Point<2> vertices[8] = {
4597 p + Point<2>(0, -1) * radius,
4598 p + Point<2>(+1, -1) * (radius / std::sqrt(2.0)),
4599 p + Point<2>(0, -1) * (radius / std::sqrt(2.0) * a),
4600 p + Point<2>(+1, -1) * (radius / std::sqrt(2.0) * a),
4601 p + Point<2>(0, +1) * (radius / std::sqrt(2.0) * a),
4602 p + Point<2>(+1, +1) * (radius / std::sqrt(2.0) * a),
4603 p + Point<2>(0, +1) * radius,
4604 p + Point<2>(+1, +1) * (radius / std::sqrt(2.0))};
4605
4606 const int cell_vertices[5][4] = {{0, 1, 2, 3},
4607 {2, 3, 4, 5},
4608 {1, 7, 3, 5},
4609 {6, 4, 7, 5}};
4610
4611 std::vector<CellData<2>> cells(4, CellData<2>());
4612
4613 for (unsigned int i = 0; i < 4; ++i)
4614 {
4615 for (unsigned int j = 0; j < 4; ++j)
4616 cells[i].vertices[j] = cell_vertices[i][j];
4617 cells[i].material_id = 0;
4618 }
4619
4622 cells,
4623 SubCellData()); // no boundary information
4624
4627
4629
4630 while (cell != end)
4631 {
4632 for (unsigned int i : GeometryInfo<2>::face_indices())
4633 {
4634 if (cell->face(i)->boundary_id() ==
4636 continue;
4637
4638 // If x is zero, then this is part of the plane
4639 if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius)
4640 {
4641 cell->face(i)->set_boundary_id(1);
4642 cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
4643 }
4644 }
4645 ++cell;
4646 }
4648 }
4649
4650
4651
4652 // Implementation for 2D only
4653 template <>
4654 void
4656 const Point<2> & center,
4657 const double inner_radius,
4658 const double outer_radius,
4659 const unsigned int n_cells,
4660 const bool colorize)
4661 {
4662 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4663 ExcInvalidRadii());
4664
4665 const double pi = numbers::PI;
4666 // determine the number of cells
4667 // for the grid. if not provided by
4668 // the user determine it such that
4669 // the length of each cell on the
4670 // median (in the middle between
4671 // the two circles) is equal to its
4672 // radial extent (which is the
4673 // difference between the two
4674 // radii)
4675 const unsigned int N =
4676 (n_cells == 0 ? static_cast<unsigned int>(
4677 std::ceil((pi * (outer_radius + inner_radius) / 2) /
4678 (outer_radius - inner_radius))) :
4679 n_cells);
4680
4681 // set up N+1 vertices on the
4682 // outer and N+1 vertices on
4683 // the inner circle. the
4684 // first N+1 ones are on the
4685 // outer one, and all are
4686 // numbered counter-clockwise
4687 std::vector<Point<2>> vertices(2 * (N + 1));
4688 for (unsigned int i = 0; i <= N; ++i)
4689 {
4690 // enforce that the x-coordinates
4691 // of the first and last point of
4692 // each half-circle are exactly
4693 // zero (contrary to what we may
4694 // compute using the imprecise
4695 // value of pi)
4696 vertices[i] =
4697 Point<2>(((i == 0) || (i == N) ? 0 : std::cos(pi * i / N - pi / 2)),
4698 std::sin(pi * i / N - pi / 2)) *
4699 outer_radius;
4700 vertices[i + N + 1] = vertices[i] * (inner_radius / outer_radius);
4701
4702 vertices[i] += center;
4703 vertices[i + N + 1] += center;
4704 }
4705
4706
4707 std::vector<CellData<2>> cells(N, CellData<2>());
4708
4709 for (unsigned int i = 0; i < N; ++i)
4710 {
4711 cells[i].vertices[0] = i;
4712 cells[i].vertices[1] = (i + 1) % (N + 1);
4713 cells[i].vertices[2] = N + 1 + i;
4714 cells[i].vertices[3] = N + 1 + ((i + 1) % (N + 1));
4715
4716 cells[i].material_id = 0;
4717 }
4718
4720
4721 if (colorize)
4722 {
4724 for (; cell != tria.end(); ++cell)
4725 {
4726 cell->face(2)->set_boundary_id(1);
4727 }
4728 tria.begin()->face(0)->set_boundary_id(3);
4729
4730 tria.last()->face(1)->set_boundary_id(2);
4731 }
4734 }
4735
4736
4737 template <>
4738 void
4740 const Point<2> & center,
4741 const double inner_radius,
4742 const double outer_radius,
4743 const unsigned int n_cells,
4744 const bool colorize)
4745 {
4746 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4747 ExcInvalidRadii());
4748
4749 const double pi = numbers::PI;
4750 // determine the number of cells
4751 // for the grid. if not provided by
4752 // the user determine it such that
4753 // the length of each cell on the
4754 // median (in the middle between
4755 // the two circles) is equal to its
4756 // radial extent (which is the
4757 // difference between the two
4758 // radii)
4759 const unsigned int N =
4760 (n_cells == 0 ? static_cast<unsigned int>(
4761 std::ceil((pi * (outer_radius + inner_radius) / 4) /
4762 (outer_radius - inner_radius))) :
4763 n_cells);
4764
4765 // set up N+1 vertices on the
4766 // outer and N+1 vertices on
4767 // the inner circle. the
4768 // first N+1 ones are on the
4769 // outer one, and all are
4770 // numbered counter-clockwise
4771 std::vector<Point<2>> vertices(2 * (N + 1));
4772 for (unsigned int i = 0; i <= N; ++i)
4773 {
4774 // enforce that the x-coordinates
4775 // of the last point is exactly
4776 // zero (contrary to what we may
4777 // compute using the imprecise
4778 // value of pi)
4779 vertices[i] = Point<2>(((i == N) ? 0 : std::cos(pi * i / N / 2)),
4780 std::sin(pi * i / N / 2)) *
4781 outer_radius;
4782 vertices[i + N + 1] = vertices[i] * (inner_radius / outer_radius);
4783
4784 vertices[i] += center;
4785 vertices[i + N + 1] += center;
4786 }
4787
4788
4789 std::vector<CellData<2>> cells(N, CellData<2>());
4790
4791 for (unsigned int i = 0; i < N; ++i)
4792 {
4793 cells[i].vertices[0] = i;
4794 cells[i].vertices[1] = (i + 1) % (N + 1);
4795 cells[i].vertices[2] = N + 1 + i;
4796 cells[i].vertices[3] = N + 1 + ((i + 1) % (N + 1));
4797
4798 cells[i].material_id = 0;
4799 }
4800
4802
4803 if (colorize)
4804 {
4806 for (; cell != tria.end(); ++cell)
4807 {
4808 cell->face(2)->set_boundary_id(1);
4809 }
4810 tria.begin()->face(0)->set_boundary_id(3);
4811
4812 tria.last()->face(1)->set_boundary_id(2);
4813 }
4814
4817 }
4818
4819
4820
4821 // Implementation for 3D only
4822 template <>
4823 void
4825 const double left,
4826 const double right,
4827 const bool colorize)
4828 {
4829 const double rl2 = (right + left) / 2;
4830 const double len = (right - left) / 2.;
4831
4832 const Point<3> vertices[20] = {
4833 Point<3>(left, left, -len / 2.), Point<3>(rl2, left, -len / 2.),
4834 Point<3>(rl2, rl2, -len / 2.), Point<3>(left, rl2, -len / 2.),
4835 Point<3>(right, left, -len / 2.), Point<3>(right, rl2, -len / 2.),
4836 Point<3>(rl2, right, -len / 2.), Point<3>(left, right, -len / 2.),
4837 Point<3>(right, right, -len / 2.), Point<3>(rl2, left, -len / 2.),
4838 Point<3>(left, left, len / 2.), Point<3>(rl2, left, len / 2.),
4839 Point<3>(rl2, rl2, len / 2.), Point<3>(left, rl2, len / 2.),
4840 Point<3>(right, left, len / 2.), Point<3>(right, rl2, len / 2.),
4841 Point<3>(rl2, right, len / 2.), Point<3>(left, right, len / 2.),
4842 Point<3>(right, right, len / 2.), Point<3>(rl2, left, len / 2.)};
4843 const int cell_vertices[4][8] = {{0, 1, 3, 2, 10, 11, 13, 12},
4844 {9, 4, 2, 5, 19, 14, 12, 15},
4845 {3, 2, 7, 6, 13, 12, 17, 16},
4846 {2, 5, 6, 8, 12, 15, 16, 18}};
4847 std::vector<CellData<3>> cells(4, CellData<3>());
4848 for (unsigned int i = 0; i < 4; ++i)
4849 {
4850 for (unsigned int j = 0; j < 8; ++j)
4851 cells[i].vertices[j] = cell_vertices[i][j];
4852 cells[i].material_id = 0;
4853 }
4856 cells,
4857 SubCellData()); // no boundary information
4858
4859 if (colorize)
4860 {
4862 cell->face(1)->set_boundary_id(1);
4863 ++cell;
4864 cell->face(0)->set_boundary_id(2);
4865 }
4866 }
4867
4868
4869
4870 // Implementation for 3D only
4871 template <>
4872 void
4874 const double left,
4875 const double right,
4876 const double thickness,
4877 const bool colorize)
4878 {
4879 Assert(left < right,
4880 ExcMessage("Invalid left-to-right bounds of enclosed hypercube"));
4881
4882 std::vector<Point<3>> vertices(64);
4883 double coords[4];
4884 coords[0] = left - thickness;
4885 coords[1] = left;
4886 coords[2] = right;
4887 coords[3] = right + thickness;
4888
4889 unsigned int k = 0;
4890 for (const double z : coords)
4891 for (const double y : coords)
4892 for (const double x : coords)
4893 vertices[k++] = Point<3>(x, y, z);
4894
4895 const types::material_id materials[27] = {21, 20, 22, 17, 16, 18, 25,
4896 24, 26, 5, 4, 6, 1, 0,
4897 2, 9, 8, 10, 37, 36, 38,
4898 33, 32, 34, 41, 40, 42};
4899
4900 std::vector<CellData<3>> cells(27);
4901 k = 0;
4902 for (unsigned int z = 0; z < 3; ++z)
4903 for (unsigned int y = 0; y < 3; ++y)
4904 for (unsigned int x = 0; x < 3; ++x)
4905 {
4906 cells[k].vertices[0] = x + 4 * y + 16 * z;
4907 cells[k].vertices[1] = x + 4 * y + 16 * z + 1;
4908 cells[k].vertices[2] = x + 4 * y + 16 * z + 4;
4909 cells[k].vertices[3] = x + 4 * y + 16 * z + 5;
4910 cells[k].vertices[4] = x + 4 * y + 16 * z + 16;
4911 cells[k].vertices[5] = x + 4 * y + 16 * z + 17;
4912 cells[k].vertices[6] = x + 4 * y + 16 * z + 20;
4913 cells[k].vertices[7] = x + 4 * y + 16 * z + 21;
4914 if (colorize)
4915 cells[k].material_id = materials[k];
4916 ++k;
4917 }
4919 cells,
4920 SubCellData()); // no boundary information
4921 }
4922
4923
4924
4925 template <>
4926 void
4928 const double radius_0,
4929 const double radius_1,
4930 const double half_length)
4931 {
4932 Assert(triangulation.n_cells() == 0,
4933 ExcMessage("The output triangulation object needs to be empty."));
4934 Assert(0 < radius_0, ExcMessage("The radii must be positive."));
4935 Assert(0 < radius_1, ExcMessage("The radii must be positive."));
4936 Assert(0 < half_length, ExcMessage("The half length must be positive."));
4937
4938 const auto n_slices = 1 + static_cast<unsigned int>(std::ceil(
4939 half_length / std::max(radius_0, radius_1)));
4940
4941 Triangulation<2> triangulation_2;
4942 GridGenerator::hyper_ball(triangulation_2, Point<2>(), radius_0);
4944 n_slices,
4945 2 * half_length,
4948 GridTools::shift(Tensor<1, 3>({-half_length, 0.0, 0.0}), triangulation);
4949 // At this point we have a cylinder. Multiply the y and z coordinates by a
4950 // factor that scales (with x) linearly between radius_0 and radius_1 to fix
4951 // the circle radii and interior points:
4952 auto shift_radii = [=](const Point<3> &p) {
4953 const double slope = (radius_1 / radius_0 - 1.0) / (2.0 * half_length);
4954 const double factor = slope * (p[0] - -half_length) + 1.0;
4955 return Point<3>(p[0], factor * p[1], factor * p[2]);
4956 };
4957 GridTools::transform(shift_radii, triangulation);
4958
4959 // Set boundary ids at -half_length to 1 and at half_length to 2. Set the
4960 // manifold id on hull faces (i.e., faces not on either end) to 0.
4961 for (const auto &face : triangulation.active_face_iterators())
4962 if (face->at_boundary())
4963 {
4964 if (std::abs(face->center()[0] - -half_length) < 1e-8 * half_length)
4965 face->set_boundary_id(1);
4966 else if (std::abs(face->center()[0] - half_length) <
4967 1e-8 * half_length)
4968 face->set_boundary_id(2);
4969 else
4970 face->set_all_manifold_ids(0);
4971 }
4972
4973 triangulation.set_manifold(0, CylindricalManifold<3>());
4974 }
4975
4976
4977 // Implementation for 3D only
4978 template <>
4979 void
4981 const double a,
4982 const double b,
4983 const bool colorize)
4984 {
4985 // we slice out the top back right
4986 // part of the cube
4987 const Point<3> vertices[26] = {
4988 // front face of the big cube
4989 Point<3>(a, a, a),
4990 Point<3>((a + b) / 2, a, a),
4991 Point<3>(b, a, a),
4992 Point<3>(a, a, (a + b) / 2),
4993 Point<3>((a + b) / 2, a, (a + b) / 2),
4994 Point<3>(b, a, (a + b) / 2),
4995 Point<3>(a, a, b),
4996 Point<3>((a + b) / 2, a, b),
4997 Point<3>(b, a, b),
4998 // middle face of the big cube
4999 Point<3>(a, (a + b) / 2, a),
5000 Point<3>((a + b) / 2, (a + b) / 2, a),
5001 Point<3>(b, (a + b) / 2, a),
5002 Point<3>(a, (a + b) / 2, (a + b) / 2),
5003 Point<3>((a + b) / 2, (a + b) / 2, (a + b) / 2),
5004 Point<3>(b, (a + b) / 2, (a + b) / 2),
5005 Point<3>(a, (a + b) / 2, b),
5006 Point<3>((a + b) / 2, (a + b) / 2, b),
5007 Point<3>(b, (a + b) / 2, b),
5008 // back face of the big cube
5009 // last (top right) point is missing
5010 Point<3>(a, b, a),
5011 Point<3>((a + b) / 2, b, a),
5012 Point<3>(b, b, a),
5013 Point<3>(a, b, (a + b) / 2),
5014 Point<3>((a + b) / 2, b, (a + b) / 2),
5015 Point<3>(b, b, (a + b) / 2),
5016 Point<3>(a, b, b),
5017 Point<3>((a + b) / 2, b, b)};
5018 const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
5019 {1, 2, 10, 11, 4, 5, 13, 14},
5020 {3, 4, 12, 13, 6, 7, 15, 16},
5021 {4, 5, 13, 14, 7, 8, 16, 17},
5022 {9, 10, 18, 19, 12, 13, 21, 22},
5023 {10, 11, 19, 20, 13, 14, 22, 23},
5024 {12, 13, 21, 22, 15, 16, 24, 25}};
5025
5026 std::vector<CellData<3>> cells(7, CellData<3>());
5027
5028 for (unsigned int i = 0; i < 7; ++i)
5029 {
5030 for (unsigned int j = 0; j < 8; ++j)
5031 cells[i].vertices[j] = cell_vertices[i][j];
5032 cells[i].material_id = 0;
5033 }
5034
5037 cells,
5038 SubCellData()); // no boundary information
5039
5040 if (colorize)
5041 {
5042 Assert(false, ExcNotImplemented());
5043 }
5044 }
5045
5046
5047
5048 // Implementation for 3D only
5049 template <>
5050 void
5052 const Point<3> & p,
5053 const double radius,
5054 const bool internal_manifold)
5055 {
5056 const double a =
5057 1. / (1 + std::sqrt(3.0)); // equilibrate cell sizes at transition
5058 // from the inner part to the radial
5059 // cells
5060 const unsigned int n_vertices = 16;
5061 const Point<3> vertices[n_vertices] = {
5062 // first the vertices of the inner
5063 // cell
5064 p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0) * a),
5065 p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0) * a),
5066 p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0) * a),
5067 p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0) * a),
5068 p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0) * a),
5069 p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0) * a),
5070 p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0) * a),
5071 p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0) * a),
5072 // now the eight vertices at
5073 // the outer sphere
5074 p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0)),
5075 p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0)),
5076 p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0)),
5077 p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0)),
5078 p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0)),
5079 p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0)),
5080 p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0)),
5081 p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0)),
5082 };
5083
5084 // one needs to draw the seven cubes to
5085 // understand what's going on here
5086 const unsigned int n_cells = 7;
5087 const int cell_vertices[n_cells][8] = {
5088 {0, 1, 4, 5, 3, 2, 7, 6}, // center
5089 {8, 9, 12, 13, 0, 1, 4, 5}, // bottom
5090 {9, 13, 1, 5, 10, 14, 2, 6}, // right
5091 {11, 10, 3, 2, 15, 14, 7, 6}, // top
5092 {8, 0, 12, 4, 11, 3, 15, 7}, // left
5093 {8, 9, 0, 1, 11, 10, 3, 2}, // front
5094 {12, 4, 13, 5, 15, 7, 14, 6}}; // back
5095
5096 std::vector<CellData<3>> cells(n_cells, CellData<3>());
5097
5098 for (unsigned int i = 0; i < n_cells; ++i)
5099 {
5100 for (const unsigned int j : GeometryInfo<3>::vertex_indices())
5101 cells[i].vertices[j] = cell_vertices[i][j];
5102 cells[i].material_id = 0;
5103 cells[i].manifold_id = i == 0 ? numbers::flat_manifold_id : 1;
5104 }
5105
5108 cells,
5109 SubCellData()); // no boundary information
5112 if (internal_manifold)
5114 }
5115
5116
5117
5118 void
5120 const unsigned int n_rotate_middle_square)
5121 {
5122 AssertThrow(n_rotate_middle_square < 4,
5123 ExcMessage("The number of rotation by pi/2 of the right square "
5124 "must be in the half-open range [0,4)."))
5125
5126 constexpr unsigned int dim = 2;
5127
5128 const unsigned int n_cells = 5;
5129 std::vector<CellData<dim>> cells(n_cells);
5130
5131 // Corner points of the cube [0,1]^2
5132 const std::vector<Point<dim>> vertices = {Point<dim>(0, 0), // 0
5133 Point<dim>(1, 0), // 1
5134 Point<dim>(0, 1), // 2
5135 Point<dim>(1, 1), // 3
5136 Point<dim>(2, 0), // 4
5137 Point<dim>(2, 1), // 5
5138 Point<dim>(3, 0), // 6
5139 Point<dim>(3, 1), // 7
5140 Point<dim>(1, -1), // 8
5141 Point<dim>(2, -1), // 9
5142 Point<dim>(1, 2), // 10
5143 Point<dim>(2, 2)}; // 11
5144
5145
5146 // consistent orientation
5147 unsigned int cell_vertices[n_cells][4] = {{0, 1, 2, 3},
5148 {1, 4, 3, 5}, // rotating cube
5149 {8, 9, 1, 4},
5150 {4, 6, 5, 7},
5151 {3, 5, 10, 11}};
5152
5153 switch (n_rotate_middle_square)
5154 {
5155 case /* rotate right square */ 1:
5156 {
5157 cell_vertices[1][0] = 4;
5158 cell_vertices[1][1] = 5;
5159 cell_vertices[1][2] = 1;
5160 cell_vertices[1][3] = 3;
5161 break;
5162 }
5163
5164 case /* rotate right square */ 2:
5165 {
5166 cell_vertices[1][0] = 5;
5167 cell_vertices[1][1] = 3;
5168 cell_vertices[1][2] = 4;
5169 cell_vertices[1][3] = 1;
5170 break;
5171 }
5172
5173 case /* rotate right square */ 3:
5174 {
5175 cell_vertices[1][0] = 3;
5176 cell_vertices[1][1] = 1;
5177 cell_vertices[1][2] = 5;
5178 cell_vertices[1][3] = 4;
5179 break;
5180 }
5181
5182 default /* 0 */:
5183 break;
5184 } // switch
5185
5186 cells.resize(n_cells, CellData<dim>());
5187
5188 for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
5189 {
5190 for (const unsigned int vertex_index :
5192 {
5193 cells[cell_index].vertices[vertex_index] =
5194 cell_vertices[cell_index][vertex_index];
5195 cells[cell_index].material_id = 0;
5196 }
5197 }
5198
5200 }
5201
5202
5203 void
5205 const bool face_orientation,
5206 const bool face_flip,
5207 const bool face_rotation,
5208 const bool manipulate_left_cube)
5209 {
5210 constexpr unsigned int dim = 3;
5211
5212 const unsigned int n_cells = 2;
5213 std::vector<CellData<dim>> cells(n_cells);
5214
5215 // Corner points of the cube [0,1]^3
5216 const std::vector<Point<dim>> vertices = {Point<dim>(0, 0, 0), // 0
5217 Point<dim>(1, 0, 0), // 1
5218 Point<dim>(0, 1, 0), // 2
5219 Point<dim>(1, 1, 0), // 3
5220 Point<dim>(0, 0, 1), // 4
5221 Point<dim>(1, 0, 1), // 5
5222 Point<dim>(0, 1, 1), // 6
5223 Point<dim>(1, 1, 1), // 7
5224 Point<dim>(2, 0, 0), // 8
5225 Point<dim>(2, 1, 0), // 9
5226 Point<dim>(2, 0, 1), // 10
5227 Point<dim>(2, 1, 1)}; // 11
5228
5229 unsigned int cell_vertices[n_cells][8] = {
5230 {0, 1, 2, 3, 4, 5, 6, 7}, // unit cube
5231 {1, 8, 3, 9, 5, 10, 7, 11}}; // shifted cube
5232
5233 // binary to case number
5234 const unsigned int this_case = 4 * static_cast<int>(face_orientation) +
5235 2 * static_cast<int>(face_flip) +
5236 static_cast<int>(face_rotation);
5237
5238 if (manipulate_left_cube)
5239 {
5240 switch (this_case)
5241 {
5242 case 0:
5243 {
5244 cell_vertices[0][0] = 1;
5245 cell_vertices[0][1] = 0;
5246 cell_vertices[0][2] = 5;
5247 cell_vertices[0][3] = 4;
5248 cell_vertices[0][4] = 3;
5249 cell_vertices[0][5] = 2;
5250 cell_vertices[0][6] = 7;
5251 cell_vertices[0][7] = 6;
5252 break;
5253 }
5254
5255 case 1:
5256 {
5257 cell_vertices[0][0] = 5;
5258 cell_vertices[0][1] = 4;
5259 cell_vertices[0][2] = 7;
5260 cell_vertices[0][3] = 6;
5261 cell_vertices[0][4] = 1;
5262 cell_vertices[0][5] = 0;
5263 cell_vertices[0][6] = 3;
5264 cell_vertices[0][7] = 2;
5265 break;
5266 }
5267
5268 case 2:
5269 {
5270 cell_vertices[0][0] = 7;
5271 cell_vertices[0][1] = 6;
5272 cell_vertices[0][2] = 3;
5273 cell_vertices[0][3] = 2;
5274 cell_vertices[0][4] = 5;
5275 cell_vertices[0][5] = 4;
5276 cell_vertices[0][6] = 1;
5277 cell_vertices[0][7] = 0;
5278 break;
5279 }
5280 case 3:
5281 {
5282 cell_vertices[0][0] = 3;
5283 cell_vertices[0][1] = 2;
5284 cell_vertices[0][2] = 1;
5285 cell_vertices[0][3] = 0;
5286 cell_vertices[0][4] = 7;
5287 cell_vertices[0][5] = 6;
5288 cell_vertices[0][6] = 5;
5289 cell_vertices[0][7] = 4;
5290 break;
5291 }
5292
5293 case 4:
5294 {
5295 cell_vertices[0][0] = 0;
5296 cell_vertices[0][1] = 1;
5297 cell_vertices[0][2] = 2;
5298 cell_vertices[0][3] = 3;
5299 cell_vertices[0][4] = 4;
5300 cell_vertices[0][5] = 5;
5301 cell_vertices[0][6] = 6;
5302 cell_vertices[0][7] = 7;
5303 break;
5304 }
5305
5306 case 5:
5307 {
5308 cell_vertices[0][0] = 2;
5309 cell_vertices[0][1] = 3;
5310 cell_vertices[0][2] = 6;
5311 cell_vertices[0][3] = 7;
5312 cell_vertices[0][4] = 0;
5313 cell_vertices[0][5] = 1;
5314 cell_vertices[0][6] = 4;
5315 cell_vertices[0][7] = 5;
5316 break;
5317 }
5318
5319 case 6:
5320 {
5321 cell_vertices[0][0] = 6;
5322 cell_vertices[0][1] = 7;
5323 cell_vertices[0][2] = 4;
5324 cell_vertices[0][3] = 5;
5325 cell_vertices[0][4] = 2;
5326 cell_vertices[0][5] = 3;
5327 cell_vertices[0][6] = 0;
5328 cell_vertices[0][7] = 1;
5329 break;
5330 }
5331
5332 case 7:
5333 {
5334 cell_vertices[0][0] = 4;
5335 cell_vertices[0][1] = 5;
5336 cell_vertices[0][2] = 0;
5337 cell_vertices[0][3] = 1;
5338 cell_vertices[0][4] = 6;
5339 cell_vertices[0][5] = 7;
5340 cell_vertices[0][6] = 2;
5341 cell_vertices[0][7] = 3;
5342 break;
5343 }
5344 } // switch
5345 }
5346 else
5347 {
5348 switch (this_case)
5349 {
5350 case 0:
5351 {
5352 cell_vertices[1][0] = 8;
5353 cell_vertices[1][1] = 1;
5354 cell_vertices[1][2] = 10;
5355 cell_vertices[1][3] = 5;
5356 cell_vertices[1][4] = 9;
5357 cell_vertices[1][5] = 3;
5358 cell_vertices[1][6] = 11;
5359 cell_vertices[1][7] = 7;
5360 break;
5361 }
5362
5363 case 1:
5364 {
5365 cell_vertices[1][0] = 10;
5366 cell_vertices[1][1] = 5;
5367 cell_vertices[1][2] = 11;
5368 cell_vertices[1][3] = 7;
5369 cell_vertices[1][4] = 8;
5370 cell_vertices[1][5] = 1;
5371 cell_vertices[1][6] = 9;
5372 cell_vertices[1][7] = 3;
5373 break;
5374 }
5375
5376 case 2:
5377 {
5378 cell_vertices[1][0] = 11;
5379 cell_vertices[1][1] = 7;
5380 cell_vertices[1][2] = 9;
5381 cell_vertices[1][3] = 3;
5382 cell_vertices[1][4] = 10;
5383 cell_vertices[1][5] = 5;
5384 cell_vertices[1][6] = 8;
5385 cell_vertices[1][7] = 1;
5386 break;
5387 }
5388
5389 case 3:
5390 {
5391 cell_vertices[1][0] = 9;
5392 cell_vertices[1][1] = 3;
5393 cell_vertices[1][2] = 8;
5394 cell_vertices[1][3] = 1;
5395 cell_vertices[1][4] = 11;
5396 cell_vertices[1][5] = 7;
5397 cell_vertices[1][6] = 10;
5398 cell_vertices[1][7] = 5;
5399 break;
5400 }
5401
5402 case 4:
5403 {
5404 cell_vertices[1][0] = 1;
5405 cell_vertices[1][1] = 8;
5406 cell_vertices[1][2] = 3;
5407 cell_vertices[1][3] = 9;
5408 cell_vertices[1][4] = 5;
5409 cell_vertices[1][5] = 10;
5410 cell_vertices[1][6] = 7;
5411 cell_vertices[1][7] = 11;
5412 break;
5413 }
5414
5415 case 5:
5416 {
5417 cell_vertices[1][0] = 5;
5418 cell_vertices[1][1] = 10;
5419 cell_vertices[1][2] = 1;
5420 cell_vertices[1][3] = 8;
5421 cell_vertices[1][4] = 7;
5422 cell_vertices[1][5] = 11;
5423 cell_vertices[1][6] = 3;
5424 cell_vertices[1][7] = 9;
5425 break;
5426 }
5427
5428 case 6:
5429 {
5430 cell_vertices[1][0] = 7;
5431 cell_vertices[1][1] = 11;
5432 cell_vertices[1][2] = 5;
5433 cell_vertices[1][3] = 10;
5434 cell_vertices[1][4] = 3;
5435 cell_vertices[1][5] = 9;
5436 cell_vertices[1][6] = 1;
5437 cell_vertices[1][7] = 8;
5438 break;
5439 }
5440
5441 case 7:
5442 {
5443 cell_vertices[1][0] = 3;
5444 cell_vertices[1][1] = 9;
5445 cell_vertices[1][2] = 7;
5446 cell_vertices[1][3] = 11;
5447 cell_vertices[1][4] = 1;
5448 cell_vertices[1][5] = 8;
5449 cell_vertices[1][6] = 5;
5450 cell_vertices[1][7] = 10;
5451 break;
5452 }
5453 } // switch
5454 }
5455
5456 cells.resize(n_cells, CellData<dim>());
5457
5458 for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
5459 {
5460 for (const unsigned int vertex_index :
5462 {
5463 cells[cell_index].vertices[vertex_index] =
5464 cell_vertices[cell_index][vertex_index];
5465 cells[cell_index].material_id = 0;
5466 }
5467 }
5468
5470 }
5471
5472
5473
5474 template <int spacedim>
5475 void
5477 const Point<spacedim> & p,
5478 const double radius)
5479 {
5480 Triangulation<spacedim> volume_mesh;
5481 GridGenerator::hyper_ball(volume_mesh, p, radius);
5482 const std::set<types::boundary_id> boundary_ids = {0};
5483 GridGenerator::extract_boundary_mesh(volume_mesh, tria, boundary_ids);
5486 }
5487
5488
5489
5490 // Implementation for 3D only
5491 template <>
5492 void
5494 const unsigned int x_subdivisions,
5495 const double radius,
5496 const double half_length)
5497 {
5498 // Copy the base from hyper_ball<3>
5499 // and transform it to yz
5500 const double d = radius / std::sqrt(2.0);
5501 const double a = d / (1 + std::sqrt(2.0));
5502
5503 std::vector<Point<3>> vertices;
5504 const double initial_height = -half_length;
5505 const double height_increment = 2. * half_length / x_subdivisions;
5506
5507 for (unsigned int rep = 0; rep < (x_subdivisions + 1); ++rep)
5508 {
5509 const double height = initial_height + height_increment * rep;
5510
5511 vertices.emplace_back(Point<3>(-d, height, -d));
5512 vertices.emplace_back(Point<3>(d, height, -d));
5513 vertices.emplace_back(Point<3>(-a, height, -a));
5514 vertices.emplace_back(Point<3>(a, height, -a));
5515 vertices.emplace_back(Point<3>(-a, height, a));
5516 vertices.emplace_back(Point<3>(a, height, a));
5517 vertices.emplace_back(Point<3>(-d, height, d));
5518 vertices.emplace_back(Point<3>(d, height, d));
5519 }
5520
5521 // Turn cylinder such that y->x
5522 for (auto &vertex : vertices)
5523 {
5524 const double h = vertex(1);
5525 vertex(1) = -vertex(0);
5526 vertex(0) = h;
5527 }
5528
5529 std::vector<std::vector<int>> cell_vertices;
5530 cell_vertices.push_back({0, 1, 8, 9, 2, 3, 10, 11});
5531 cell_vertices.push_back({0, 2, 8, 10, 6, 4, 14, 12});
5532 cell_vertices.push_back({2, 3, 10, 11, 4, 5, 12, 13});
5533 cell_vertices.push_back({1, 7, 9, 15, 3, 5, 11, 13});
5534 cell_vertices.push_back({6, 4, 14, 12, 7, 5, 15, 13});
5535
5536 for (unsigned int rep = 1; rep < x_subdivisions; ++rep)
5537 {
5538 for (unsigned int i = 0; i < 5; ++i)
5539 {
5540 std::vector<int> new_cell_vertices(8);
5541 for (unsigned int j = 0; j < 8; ++j)
5542 new_cell_vertices[j] = cell_vertices[i][j] + 8 * rep;
5543 cell_vertices.push_back(new_cell_vertices);
5544 }
5545 }
5546
5547 unsigned int n_cells = x_subdivisions * 5;
5548
5549 std::vector<CellData<3>> cells(n_cells, CellData<3>());
5550
5551 for (unsigned int i = 0; i < n_cells; ++i)
5552 {
5553 for (unsigned int j = 0; j < 8; ++j)
5554 cells[i].vertices[j] = cell_vertices[i][j];
5555 cells[i].material_id = 0;
5556 }
5557
5560 cells,
5561 SubCellData()); // no boundary information
5562
5563 // set boundary indicators for the
5564 // faces at the ends to 1 and 2,
5565 // respectively. note that we also
5566 // have to deal with those lines
5567 // that are purely in the interior
5568 // of the ends. we determine whether
5569 // an edge is purely in the
5570 // interior if one of its vertices
5571 // is at coordinates '+-a' as set
5572 // above
5574
5575 // Tolerance is calculated using the minimal length defining
5576 // the cylinder
5577 const double tolerance = 1e-5 * std::min(radius, half_length);
5578
5579 for (const auto &cell : tria.cell_iterators())
5580 for (unsigned int i : GeometryInfo<3>::face_indices())
5581 if (cell->at_boundary(i))
5582 {
5583 if (cell->face(i)->center()(0) > half_length - tolerance)
5584 {
5585 cell->face(i)->set_boundary_id(2);
5586 cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
5587
5588 for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5589 ++e)
5590 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5591 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5592 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5593 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5594 {
5595 cell->face(i)->line(e)->set_boundary_id(2);
5596 cell->face(i)->line(e)->set_manifold_id(
5598 }
5599 }
5600 else if (cell->face(i)->center()(0) < -half_length + tolerance)
5601 {
5602 cell->face(i)->set_boundary_id(1);
5603 cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
5604
5605 for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5606 ++e)
5607 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5608 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5609 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5610 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5611 {
5612 cell->face(i)->line(e)->set_boundary_id(1);
5613 cell->face(i)->line(e)->set_manifold_id(
5615 }
5616 }
5617 }
5619 }
5620
5621 // Implementation for 3D only
5622 template <>
5623 void
5625 const double radius,
5626 const double half_length)
5627 {
5628 subdivided_cylinder(tria, 2, radius, half_length);
5629 }
5630
5631 template <>
5632 void
5634 const Point<3> & center,
5635 const double radius)
5636 {
5637 const unsigned int dim = 3;
5638
5639 // the parameters a (intersection on the octant lines from center), b
5640 // (intersection within the octant faces) and c (position inside the
5641 // octant) have been derived by equilibrating the minimal singular value
5642 // of the Jacobian of the four cells around the center point c and, as a
5643 // secondary measure, to minimize the aspect ratios defined as the maximal
5644 // divided by the minimal singular values throughout cells
5645 const double a = 0.528;
5646 const double b = 0.4533;
5647 const double c = 0.3752;
5648 const Point<dim> vertices[15] = {
5649 center + Point<dim>(0, 0, 0) * radius,
5650 center + Point<dim>(+1, 0, 0) * radius,
5651 center + Point<dim>(+1, 0, 0) * (radius * a),
5652 center + Point<dim>(0, +1, 0) * (radius * a),
5653 center + Point<dim>(+1, +1, 0) * (radius * b),
5654 center + Point<dim>(0, +1, 0) * radius,
5655 center + Point<dim>(+1, +1, 0) * radius / std::sqrt(2.0),
5656 center + Point<dim>(0, 0, 1) * radius * a,
5657 center + Point<dim>(+1, 0, 1) * radius / std::sqrt(2.0),
5658 center + Point<dim>(+1, 0, 1) * (radius * b),
5659 center + Point<dim>(0, +1, 1) * (radius * b),
5660 center + Point<dim>(+1, +1, 1) * (radius * c),
5661 center + Point<dim>(0, +1, 1) * radius / std::sqrt(2.0),
5662 center + Point<dim>(+1, +1, 1) * (radius / (std::sqrt(3.0))),
5663 center + Point<dim>(0, 0, 1) * radius};
5664 const int cell_vertices[4][8] = {{0, 2, 3, 4, 7, 9, 10, 11},
5665 {1, 6, 2, 4, 8, 13, 9, 11},
5666 {5, 3, 6, 4, 12, 10, 13, 11},
5667 {7, 9, 10, 11, 14, 8, 12, 13}};
5668
5669 std::vector<CellData<dim>> cells(4, CellData<dim>());
5670
5671 for (unsigned int i = 0; i < 4; ++i)
5672 {
5673 for (unsigned int j = 0; j < 8; ++j)
5674 cells[i].vertices[j] = cell_vertices[i][j];
5675 cells[i].material_id = 0;
5676 }
5677
5680 cells,
5681 SubCellData()); // no boundary information
5682
5685
5687 while (cell != end)
5688 {
5689 for (unsigned int i : GeometryInfo<dim>::face_indices())
5690 {
5691 if (cell->face(i)->boundary_id() ==
5693 continue;
5694
5695 // If x,y or z is zero, then this is part of the plane
5696 if (cell->face(i)->center()(0) < center(0) + 1.e-5 * radius ||
5697 cell->face(i)->center()(1) < center(1) + 1.e-5 * radius ||
5698 cell->face(i)->center()(2) < center(2) + 1.e-5 * radius)
5699 {
5700 cell->face(i)->set_boundary_id(1);
5701 cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
5702 // also set the boundary indicators of the bounding lines,
5703 // unless both vertices are on the perimeter
5704 for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5705 ++j)
5706 {
5707 const Point<3> line_vertices[2] = {
5708 cell->face(i)->line(j)->vertex(0),
5709 cell->face(i)->line(j)->vertex(1)};
5710 if ((std::fabs(line_vertices[0].distance(center) - radius) >
5711 1e-5 * radius) ||
5712 (std::fabs(line_vertices[1].distance(center) - radius) >
5713 1e-5 * radius))
5714 {
5715 cell->face(i)->line(j)->set_boundary_id(1);
5716 cell->face(i)->line(j)->set_manifold_id(
5718 }
5719 }
5720 }
5721 }
5722 ++cell;
5723 }
5725 }
5726
5727
5728
5729 // Implementation for 3D only
5730 template <>
5731 void
5733 const Point<3> & center,
5734 const double radius)
5735 {
5736 // These are for the two lower squares
5737 const double d = radius / std::sqrt(2.0);
5738 const double a = d / (1 + std::sqrt(2.0));
5739 // These are for the two upper square
5740 const double b = a / 2.0;
5741 const double c =