ExaDG
Loading...
Searching...
No Matches
one_sided_cylindrical_manifold.h
1/* ______________________________________________________________________
2 *
3 * ExaDG - High-Order Discontinuous Galerkin for the Exa-Scale
4 *
5 * Copyright (C) 2021 by the ExaDG authors
6 *
7 * This program is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see <https://www.gnu.org/licenses/>.
19 * ______________________________________________________________________
20 */
21
22#ifndef EXADG_GRID_ONE_SIDED_CYLINDRICAL_MANIFOLD_H_
23#define EXADG_GRID_ONE_SIDED_CYLINDRICAL_MANIFOLD_H_
24
25// deal.II
26#include <deal.II/grid/manifold_lib.h>
27
28namespace ExaDG
29{
38template<int dim>
39class OneSidedCylindricalManifold : public dealii::ChartManifold<dim, dim, dim>
40{
41public:
42 OneSidedCylindricalManifold(dealii::Triangulation<dim> const & tria_in,
43 typename dealii::Triangulation<dim>::cell_iterator const & cell_in,
44 unsigned int const face_in,
45 dealii::Point<dim> const & center_in)
46 : alpha(1.0), radius(1.0), tria(tria_in), cell(cell_in), face(face_in), center(center_in)
47 {
48 AssertThrow(tria.all_reference_cells_are_hyper_cube(),
49 dealii::ExcMessage("This class is only implemented for hypercube elements."));
50
51 AssertThrow(face <= 3,
52 dealii::ExcMessage(
53 "One sided spherical manifold can only be applied to face f=0,1,2,3."));
54
55 // get center coordinates in x1-x2 plane
56 x_C[0] = center[0];
57 x_C[1] = center[1];
58
59 // determine x_1 and x_2 which denote the end points of the face that is
60 // subject to the spherical manifold.
61 dealii::Point<dim> x_1, x_2;
62 x_1 = cell->vertex(get_vertex_id(0));
63 x_2 = cell->vertex(get_vertex_id(1));
64
65 dealii::Point<2> x_1_2d = dealii::Point<2>(x_1[0], x_1[1]);
66 dealii::Point<2> x_2_2d = dealii::Point<2>(x_2[0], x_2[1]);
67
68 initialize(x_1_2d, x_2_2d);
69 }
70
71 void
72 initialize(dealii::Point<2> const & x_1, dealii::Point<2> const & x_2)
73 {
74 double const tol = 1.e-12;
75
76 v_1 = x_1 - x_C;
77 v_2 = x_2 - x_C,
78
79 // calculate radius of spherical manifold
80 radius = v_1.norm();
81
82 // check correctness of geometry and parameters
83 double radius_check = v_2.norm();
84 AssertThrow(std::abs(radius - radius_check) < tol * radius,
85 dealii::ExcMessage(
86 "Invalid geometry parameters. To apply a spherical manifold both "
87 "end points of the face must have the same distance from the center."));
88
89 // normalize v_1 and v_2
90 v_1 /= v_1.norm();
91 v_2 /= v_2.norm();
92
93 // calculate angle between v_1 and v_2
94 alpha = std::acos(v_1 * v_2);
95
96 // calculate vector that is perpendicular to v_1 in plane that is spanned by v_1 and v_2
97 normal = v_2 - (v_2 * v_1) * v_1;
98
99 AssertThrow(normal.norm() > tol, dealii::ExcMessage("Vector must not have length 0."));
100
101 normal /= normal.norm();
102 }
103
104 /*
105 * push_forward operation that maps point xi in reference coordinates [0,1]^d to
106 * point x in physical coordinates
107 */
108 dealii::Point<dim>
109 push_forward(dealii::Point<dim> const & xi) const override
110 {
111 dealii::Point<dim> x;
112
113 // standard mapping from reference space to physical space using d-linear shape functions
114 for(unsigned int const v : cell->vertex_indices())
115 {
116 double shape_function_value = cell->reference_cell().d_linear_shape_function(xi, v);
117 x += shape_function_value * cell->vertex(v);
118 }
119
120 // Add contribution of spherical manifold.
121 // Here, we only operate in the xi1-xi2 plane.
122
123 // set xi_face, xi_other to xi[0],xi[1] depending on the face that is subject to the manifold
124 unsigned int index_face = get_index_face();
125 unsigned int index_other = get_index_other();
126 double const xi_face = xi[index_face];
127 double const xi_other = xi[index_other];
128
129 // calculate deformation related to the spherical manifold
130 double beta = xi_face * alpha;
131
132 dealii::Tensor<1, 2> direction;
133 direction = std::cos(beta) * v_1 + std::sin(beta) * normal;
134
135 Assert(std::abs(direction.norm() - 1.0) < 1.e-12,
136 dealii::ExcMessage("Vector must have length 1."));
137
138 // calculate point x_S on spherical manifold
139 dealii::Tensor<1, 2> x_S;
140 x_S = x_C + radius * direction;
141
142 // calculate displacement as compared to straight sided quadrilateral element
143 // on the face that is subject to the manifold
144 dealii::Tensor<1, 2> displ, x_lin;
145 for(unsigned int v : dealii::ReferenceCells::template get_hypercube<1>().vertex_indices())
146 {
147 double shape_function_value =
148 dealii::ReferenceCells::template get_hypercube<1>().d_linear_shape_function(
149 dealii::Point<1>(xi_face), v);
150
151 unsigned int vertex_id = get_vertex_id(v);
152 dealii::Point<dim> vertex = cell->vertex(vertex_id);
153
154 x_lin[0] += shape_function_value * vertex[0];
155 x_lin[1] += shape_function_value * vertex[1];
156 }
157
158 displ = x_S - x_lin;
159
160 // deformation decreases linearly in the second (other) direction
161 dealii::Point<1> xi_other_1d = dealii::Point<1>(xi_other);
162 unsigned int index_1d = get_index_1d();
163 double fading_value =
164 dealii::ReferenceCells::template get_hypercube<1>().d_linear_shape_function(xi_other_1d,
165 index_1d);
166 x[0] += fading_value * displ[0];
167 x[1] += fading_value * displ[1];
168
169 Assert(dealii::numbers::is_finite(x.norm_square()), dealii::ExcMessage("Invalid point found"));
170
171 return x;
172 }
173
174 /*
175 * Calculate vertex_id of 2d object (cell in 2d, face4 in 3d)
176 * given the vertex_id of the 1d object (vertex_id_1d = 0,1).
177 */
178 unsigned int
179 get_vertex_id(unsigned int vertex_id_1d) const
180 {
181 unsigned int vertex_id = 0;
182
183 if(face == 0)
184 vertex_id = 2 * vertex_id_1d;
185 else if(face == 1)
186 vertex_id = 1 + 2 * vertex_id_1d;
187 else if(face == 2)
188 vertex_id = vertex_id_1d;
189 else if(face == 3)
190 vertex_id = 2 + vertex_id_1d;
191
192 return vertex_id;
193 }
194
195 /*
196 * Calculate index of 1d linear shape function (0 or 1)
197 * that takes a value of 1 on the specified face.
198 */
199 unsigned int
200 get_index_1d() const
201 {
202 unsigned int index_1d = 0;
203
204 if(face == 0 or face == 2)
205 index_1d = 0;
206 else if(face == 1 or face == 3)
207 index_1d = 1;
208 else
209 Assert(false, dealii::ExcMessage("Face ID is invalid."));
210
211 return index_1d;
212 }
213
214 /*
215 * Calculate which xi-coordinate corresponds to the
216 * tangent direction of the respective face
217 */
218 unsigned int
219 get_index_face() const
220 {
221 unsigned int index_face = 0;
222
223 if(face == 0 or face == 1)
224 index_face = 1;
225 else if(face == 2 or face == 3)
226 index_face = 0;
227 else
228 Assert(false, dealii::ExcMessage("Face ID is invalid."));
229
230 return index_face;
231 }
232
233 /*
234 * Calculate which xi-coordinate corresponds to
235 * the normal direction of the respective face
236 * in xi1-xi2-plane.
237 */
238 unsigned int
239 get_index_other() const
240 {
241 return 1 - get_index_face();
242 }
243
244 /*
245 * This function calculates the inverse Jacobi matrix dx/d(xi) = d(phi(xi))/d(xi) of
246 * the push-forward operation phi: [0,1]^d -> R^d: xi -> x = phi(xi)
247 * We assume that the gradient of the standard bilinear shape functions is sufficient
248 * to find the solution.
249 */
250 dealii::Tensor<2, dim>
251 get_inverse_jacobian(dealii::Point<dim> const & xi) const
252 {
253 dealii::Tensor<2, dim> jacobian;
254
255 // standard mapping from reference space to physical space using d-linear shape functions
256 for(unsigned int const v : cell->vertex_indices())
257 {
258 dealii::Tensor<1, dim> shape_function_gradient =
259 cell->reference_cell().d_linear_shape_function_gradient(xi, v);
260 jacobian += outer_product(cell->vertex(v), shape_function_gradient);
261 }
262
263 return invert(jacobian);
264 }
265
266 /*
267 * pull_back operation that maps point x in physical coordinates
268 * to point xi in reference coordinates [0,1]^d using the
269 * push_forward operation and Newton's method
270 */
271 dealii::Point<dim>
272 pull_back(dealii::Point<dim> const & x) const override
273 {
274 dealii::Point<dim> xi;
275 dealii::Tensor<1, dim> residual = push_forward(xi) - x;
276 dealii::Tensor<1, dim> delta_xi;
277
278 // Newton method to solve nonlinear pull_back operation
279 unsigned int n_iter = 0, MAX_ITER = 100;
280 double const TOL = 1.e-12;
281 while(residual.norm() > TOL and n_iter < MAX_ITER)
282 {
283 // multiply by -1.0, i.e., shift residual to the rhs
284 residual *= -1.0;
285
286 // solve linear problem
287 delta_xi = get_inverse_jacobian(xi) * residual;
288
289 // add increment
290 xi += delta_xi;
291
292 // make sure that xi is in the valid range [0,1]^d
293 if(xi[0] < 0.0)
294 xi[0] = 0.0;
295 else if(xi[0] > 1.0)
296 xi[0] = 1.0;
297
298 if(xi[1] < 0.0)
299 xi[1] = 0.0;
300 else if(xi[1] > 1.0)
301 xi[1] = 1.0;
302
303 // evaluate residual
304 residual = push_forward(xi) - x;
305
306 // increment counter
307 ++n_iter;
308 }
309
310 Assert(n_iter < MAX_ITER,
311 dealii::ExcMessage("Newton solver did not converge to given tolerance. "
312 "Maximum number of iterations exceeded."));
313
314 Assert(xi[0] >= 0.0 and xi[0] <= 1.0,
315 dealii::ExcMessage("Pull back operation generated invalid xi[0] values."));
316
317 Assert(xi[1] >= 0.0 and xi[1] <= 1.0,
318 dealii::ExcMessage("Pull back operation generated invalid xi[1] values."));
319
320 return xi;
321 }
322
323 std::unique_ptr<dealii::Manifold<dim>>
324 clone() const override
325 {
326 return std::make_unique<OneSidedCylindricalManifold<dim>>(tria, cell, face, center);
327 }
328
329private:
330 dealii::Point<2> x_C;
331 dealii::Tensor<1, 2> v_1;
332 dealii::Tensor<1, 2> v_2;
333 dealii::Tensor<1, 2> normal;
334 double alpha;
335 double radius;
336
337 dealii::Triangulation<dim> const & tria;
338 typename dealii::Triangulation<dim>::cell_iterator cell;
339 unsigned int face;
340 dealii::Point<dim> center;
341};
342
351template<int dim>
352class OneSidedConicalManifold : public dealii::ChartManifold<dim, dim, dim>
353{
354public:
355 OneSidedConicalManifold(dealii::Triangulation<dim> const & tria_in,
356 typename dealii::Triangulation<dim>::cell_iterator const & cell_in,
357 unsigned int const face_in,
358 dealii::Point<dim> const & center_in,
359 double const r_0_in,
360 double const r_1_in)
361 : alpha(1.0),
362 tria(tria_in),
363 cell(cell_in),
364 face(face_in),
365 center(center_in),
366 r_0(r_0_in),
367 r_1(r_1_in)
368 {
369 AssertThrow(tria.all_reference_cells_are_hyper_cube(),
370 dealii::ExcMessage("This class is only implemented for hypercube elements."));
371
372 AssertThrow(dim == 3,
373 dealii::ExcMessage("OneSidedConicalManifold can only be used for 3D problems."));
374
375 AssertThrow(face <= 3,
376 dealii::ExcMessage(
377 "One sided spherical manifold can only be applied to face f=0,1,2,3."));
378
379 // get center coordinates in x1-x2 plane
380 x_C[0] = center[0];
381 x_C[1] = center[1];
382
383 // determine x_1 and x_2 which denote the end points of the face that is
384 // subject to the spherical manifold.
385 dealii::Point<dim> x_1, x_2;
386 x_1 = cell->vertex(get_vertex_id(0));
387 x_2 = cell->vertex(get_vertex_id(1));
388
389 dealii::Point<2> x_1_2d = dealii::Point<2>(x_1[0], x_1[1]);
390 dealii::Point<2> x_2_2d = dealii::Point<2>(x_2[0], x_2[1]);
391
392 initialize(x_1_2d, x_2_2d);
393 }
394
395 void
396 initialize(dealii::Point<2> const & x_1, dealii::Point<2> const & x_2)
397 {
398 double const tol = 1.e-12;
399
400 v_1 = x_1 - x_C;
401 v_2 = x_2 - x_C,
402
403 // calculate radius of spherical manifold
404 r_0 = v_1.norm();
405
406 // check correctness of geometry and parameters
407 double radius_check = v_2.norm();
408
409 AssertThrow(std::abs(r_0 - radius_check) < tol * r_0,
410 dealii::ExcMessage(
411 "Invalid geometry parameters. To apply a spherical manifold both "
412 "end points of the face must have the same distance from the center."));
413
414 // normalize v_1 and v_2
415 v_1 /= v_1.norm();
416 v_2 /= v_2.norm();
417
418 // calculate angle between v_1 and v_2
419 alpha = std::acos(v_1 * v_2);
420
421 // calculate vector that is perpendicular to v_1 in plane that is spanned by v_1 and v_2
422 normal = v_2 - (v_2 * v_1) * v_1;
423
424 AssertThrow(normal.norm() > tol, dealii::ExcMessage("Vector must not have length 0."));
425
426 normal /= normal.norm();
427 }
428
429 /*
430 * push_forward operation that maps point xi in reference coordinates [0,1]^d to
431 * point x in physical coordinates
432 */
433 dealii::Point<dim>
434 push_forward(dealii::Point<dim> const & xi) const override
435 {
436 dealii::Point<dim> x;
437
438 // standard mapping from reference space to physical space using d-linear shape functions
439 for(unsigned int const v : cell->vertex_indices())
440 {
441 double shape_function_value = cell->reference_cell().d_linear_shape_function(xi, v);
442 x += shape_function_value * cell->vertex(v);
443 }
444
445 // Add contribution of conical manifold.
446 // Here, we only operate in the xi1-xi2 plane.
447
448 // set xi_face, xi_other to xi[0],xi[1] depending on the face that is subject to the manifold
449 unsigned int index_face = get_index_face();
450 unsigned int index_other = get_index_other();
451 double const xi_face = xi[index_face];
452 double const xi_other = xi[index_other];
453
454 // calculate deformation related to the conical manifold
455 double beta = xi_face * alpha;
456
457 dealii::Tensor<1, 2> direction;
458 direction = std::cos(beta) * v_1 + std::sin(beta) * normal;
459
460 Assert(std::abs(direction.norm() - 1.0) < 1.e-12,
461 dealii::ExcMessage("Vector must have length 1."));
462
463 // calculate point x_S on spherical manifold
464 dealii::Tensor<1, 2> x_S;
465 x_S = x_C + r_0 * direction;
466
467 // calculate displacement as compared to straight sided quadrilateral element
468 // on the face that is subject to the manifold
469 dealii::Tensor<1, 2> displ, x_lin;
470 for(unsigned int const v : dealii::ReferenceCells::template get_hypercube<1>().vertex_indices())
471 {
472 double shape_function_value =
473 dealii::ReferenceCells::template get_hypercube<1>().d_linear_shape_function(
474 dealii::Point<1>(xi_face), v);
475
476 unsigned int vertex_id = get_vertex_id(v);
477 dealii::Point<dim> vertex = cell->vertex(vertex_id);
478
479 x_lin[0] += shape_function_value * vertex[0];
480 x_lin[1] += shape_function_value * vertex[1];
481 }
482
483 displ = x_S - x_lin;
484
485 // conical manifold
486 displ *= (1 - xi[2] * (r_0 - r_1) / r_0);
487
488 // deformation decreases linearly in the second (other) direction
489 dealii::Point<1> xi_other_1d = dealii::Point<1>(xi_other);
490 unsigned int index_1d = get_index_1d();
491 double fading_value =
492 dealii::ReferenceCells::template get_hypercube<1>().d_linear_shape_function(xi_other_1d,
493 index_1d);
494 x[0] += fading_value * displ[0];
495 x[1] += fading_value * displ[1];
496
497 Assert(dealii::numbers::is_finite(x.norm_square()), dealii::ExcMessage("Invalid point found"));
498
499 return x;
500 }
501
502 /*
503 * Calculate vertex_id of 2d object (cell in 2d, face4 in 3d)
504 * given the vertex_id of the 1d object (vertex_id_1d = 0,1).
505 */
506 unsigned int
507 get_vertex_id(unsigned int vertex_id_1d) const
508 {
509 unsigned int vertex_id = 0;
510
511 if(face == 0)
512 vertex_id = 2 * vertex_id_1d;
513 else if(face == 1)
514 vertex_id = 1 + 2 * vertex_id_1d;
515 else if(face == 2)
516 vertex_id = vertex_id_1d;
517 else if(face == 3)
518 vertex_id = 2 + vertex_id_1d;
519
520 return vertex_id;
521 }
522
523 /*
524 * Calculate index of 1d linear shape function (0 or 1)
525 * that takes a value of 1 on the specified face.
526 */
527 unsigned int
528 get_index_1d() const
529 {
530 unsigned int index_1d = 0;
531
532 if(face == 0 or face == 2)
533 index_1d = 0;
534 else if(face == 1 or face == 3)
535 index_1d = 1;
536 else
537 Assert(false, dealii::ExcMessage("Face ID is invalid."));
538
539 return index_1d;
540 }
541
542 /*
543 * Calculate which xi-coordinate corresponds to the
544 * tangent direction of the respective face
545 */
546 unsigned int
547 get_index_face() const
548 {
549 unsigned int index_face = 0;
550
551 if(face == 0 or face == 1)
552 index_face = 1;
553 else if(face == 2 or face == 3)
554 index_face = 0;
555 else
556 Assert(false, dealii::ExcMessage("Face ID is invalid."));
557
558 return index_face;
559 }
560
561 /*
562 * Calculate which xi-coordinate corresponds to
563 * the normal direction of the respective face
564 * in xi1-xi2-plane.
565 */
566 unsigned int
567 get_index_other() const
568 {
569 return 1 - get_index_face();
570 }
571
572 /*
573 * This function calculates the inverse Jacobi matrix dx/d(xi) = d(phi(xi))/d(xi) of
574 * the push-forward operation phi: [0,1]^d -> R^d: xi -> x = phi(xi)
575 * We assume that the gradient of the standard bilinear shape functions is sufficient
576 * to find the solution.
577 */
578 dealii::Tensor<2, dim>
579 get_inverse_jacobian(dealii::Point<dim> const & xi) const
580 {
581 dealii::Tensor<2, dim> jacobian;
582
583 // standard mapping from reference space to physical space using d-linear shape functions
584 for(unsigned int const v : cell->vertex_indices())
585 {
586 dealii::Tensor<1, dim> shape_function_gradient =
587 cell->reference_cell().d_linear_shape_function_gradient(xi, v);
588 jacobian += outer_product(cell->vertex(v), shape_function_gradient);
589 }
590
591 return invert(jacobian);
592 }
593
594 /*
595 * pull_back operation that maps point x in physical coordinates
596 * to point xi in reference coordinates [0,1]^d using the
597 * push_forward operation and Newton's method
598 */
599 dealii::Point<dim>
600 pull_back(dealii::Point<dim> const & x) const override
601 {
602 dealii::Point<dim> xi;
603 dealii::Tensor<1, dim> residual = push_forward(xi) - x;
604 dealii::Tensor<1, dim> delta_xi;
605
606 // Newton method to solve nonlinear pull_back operation
607 unsigned int n_iter = 0, MAX_ITER = 100;
608 double const TOL = 1.e-12;
609 while(residual.norm() > TOL and n_iter < MAX_ITER)
610 {
611 // multiply by -1.0, i.e., shift residual to the rhs
612 residual *= -1.0;
613
614 // solve linear problem
615 delta_xi = get_inverse_jacobian(xi) * residual;
616
617 // add increment
618 xi += delta_xi;
619
620 // make sure that xi is in the valid range [0,1]^d
621 if(xi[0] < 0.0)
622 xi[0] = 0.0;
623 else if(xi[0] > 1.0)
624 xi[0] = 1.0;
625
626 if(xi[1] < 0.0)
627 xi[1] = 0.0;
628 else if(xi[1] > 1.0)
629 xi[1] = 1.0;
630
631 if(xi[2] < 0.0)
632 xi[2] = 0.0;
633 else if(xi[2] > 1.0)
634 xi[2] = 1.0;
635
636 // evaluate residual
637 residual = push_forward(xi) - x;
638
639 // increment counter
640 ++n_iter;
641 }
642
643 Assert(n_iter < MAX_ITER,
644 dealii::ExcMessage("Newton solver did not converge to given tolerance. "
645 "Maximum number of iterations exceeded."));
646
647 Assert(xi[0] >= 0.0 and xi[0] <= 1.0,
648 dealii::ExcMessage("Pull back operation generated invalid xi[0] values."));
649
650 Assert(xi[1] >= 0.0 and xi[1] <= 1.0,
651 dealii::ExcMessage("Pull back operation generated invalid xi[1] values."));
652
653 Assert(xi[2] >= 0.0 and xi[2] <= 1.0,
654 dealii::ExcMessage("Pull back operation generated invalid xi[2] values."));
655
656 return xi;
657 }
658
659 std::unique_ptr<dealii::Manifold<dim>>
660 clone() const override
661 {
662 return std::make_unique<OneSidedConicalManifold<dim>>(tria, cell, face, center, r_0, r_1);
663 }
664
665
666private:
667 dealii::Point<2> x_C;
668 dealii::Tensor<1, 2> v_1;
669 dealii::Tensor<1, 2> v_2;
670 dealii::Tensor<1, 2> normal;
671 double alpha;
672
673 dealii::Triangulation<dim> const & tria;
674 typename dealii::Triangulation<dim>::cell_iterator cell;
675 unsigned int face;
676
677 dealii::Point<dim> center;
678
679 // radius of cone at xi_3 = 0 (-> r_0) and at xi_3 = 1 (-> r_1)
680 double r_0, r_1;
681};
682
683
688template<int dim, int spacedim = dim>
689class MyCylindricalManifold : public dealii::ChartManifold<dim, spacedim, spacedim>
690{
691public:
692 MyCylindricalManifold(dealii::Point<spacedim> const center_in)
693 : dealii::ChartManifold<dim, spacedim, spacedim>(
694 MyCylindricalManifold<dim, spacedim>::get_periodicity()),
695 center(center_in)
696 {
697 }
698
699 dealii::Tensor<1, spacedim>
700 get_periodicity()
701 {
702 dealii::Tensor<1, spacedim> periodicity;
703
704 // angle theta is 2*pi periodic
705 periodicity[1] = 2 * dealii::numbers::PI;
706 return periodicity;
707 }
708
709 dealii::Point<spacedim>
710 push_forward(dealii::Point<spacedim> const & ref_point) const override
711 {
712 double const radius = ref_point[0];
713 double const theta = ref_point[1];
714
715 Assert(ref_point[0] >= 0.0, dealii::ExcMessage("Radius must be positive."));
716
717 dealii::Point<spacedim> space_point;
718 if(radius > 1e-10)
719 {
720 AssertThrow(spacedim == 2 or spacedim == 3,
721 dealii::ExcMessage("Only implemented for 2D and 3D case."));
722
723 space_point[0] = radius * cos(theta);
724 space_point[1] = radius * sin(theta);
725
726 if(spacedim == 3)
727 space_point[2] = ref_point[2];
728 }
729
730 return space_point + center;
731 }
732
733 dealii::Point<spacedim>
734 pull_back(dealii::Point<spacedim> const & space_point) const override
735 {
736 dealii::Tensor<1, spacedim> vector;
737 vector[0] = space_point[0] - center[0];
738 vector[1] = space_point[1] - center[1];
739 // for the 3D case: vector[2] will always be 0.
740
741 double const radius = vector.norm();
742
743 dealii::Point<spacedim> ref_point;
744 ref_point[0] = radius;
745 ref_point[1] = atan2(vector[1], vector[0]);
746 if(ref_point[1] < 0)
747 ref_point[1] += 2.0 * dealii::numbers::PI;
748 if(spacedim == 3)
749 ref_point[2] = space_point[2];
750
751 return ref_point;
752 }
753
754 std::unique_ptr<dealii::Manifold<dim>>
755 clone() const override
756 {
757 return std::make_unique<MyCylindricalManifold<dim, spacedim>>(center);
758 }
759
760
761private:
762 dealii::Point<dim> center;
763};
764
765} // namespace ExaDG
766
767#endif /* EXADG_GRID_ONE_SIDED_CYLINDRICAL_MANIFOLD_H_ */
ExaDG
Definition driver.cpp:33
Generated by doxygen 1.13.2