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