eighth_sphere_domain.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Include guards
27#ifndef OOMPH_EIGHTH_SPHERE_DOMAIN_HEADER
28#define OOMPH_EIGHTH_SPHERE_DOMAIN_HEADER
29
30// Generic oomph-lib includes
31#include "../generic/quadtree.h"
32#include "../generic/domain.h"
33#include "../generic/geom_objects.h"
34
35namespace oomph
36{
37 //=================================================================
38 /// Eighth sphere as domain. Domain is
39 /// parametrised by four macro elements
40 //=================================================================
41 class EighthSphereDomain : public Domain
42 {
43 public:
44 /// Constructor: Pass the radius of the sphere.
45 EighthSphereDomain(const double& radius) : Radius(radius)
46 {
47 // There are four macro elements
48 unsigned nmacro = 4;
49
50 // Resize
51 Macro_element_pt.resize(nmacro);
52
53 // Create the 3D Q macro elements
54 for (unsigned i = 0; i < nmacro; i++)
55 {
56 Macro_element_pt[i] = new QMacroElement<3>(this, i);
57 }
58 }
59
60 /// Broken copy constructor
61 EighthSphereDomain(const EighthSphereDomain&) = delete;
62
63 /// Broken assignment operator
64 void operator=(const EighthSphereDomain&) = delete;
65
66 /// Destructor: Empty; cleanup done in base class
67 ~EighthSphereDomain() {}
68
69
70 /// Vector representation of the imacro-th macro element
71 /// boundary idirect (L/R/D/U/B/F) at time level t
72 /// (t=0: present; t>0: previous):
73 /// f(s).
74 void macro_element_boundary(const unsigned& t,
75 const unsigned& imacro,
76 const unsigned& idirect,
77 const Vector<double>& s,
78 Vector<double>& f)
79 {
80 using namespace OcTreeNames;
81
82#ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
83 // Warn about time argument being moved to the front
84 OomphLibWarning(
85 "Order of function arguments has changed between versions 0.8 and 0.85",
86 "EighthSphereDomain::macro_element_boundary(...)",
87 OOMPH_EXCEPTION_LOCATION);
88#endif
89
90 // Which macro element?
91 // --------------------
92 switch (imacro)
93 {
94 // Macro element 0: Central box
95 case 0:
96
97 if (idirect == L)
98 {
99 r_centr_L(t, s, f);
100 }
101 else if (idirect == R)
102 {
103 r_centr_R(t, s, f);
104 }
105 else if (idirect == D)
106 {
107 r_centr_D(t, s, f);
108 }
109 else if (idirect == U)
110 {
111 r_centr_U(t, s, f);
112 }
113 else if (idirect == B)
114 {
115 r_centr_B(t, s, f);
116 }
117 else if (idirect == F)
118 {
119 r_centr_F(t, s, f);
120 }
121 else
122 {
123 std::ostringstream error_message;
124 error_message << "idirect is " << OcTree::Direct_string[idirect]
125 << "not one of L, R, U, D, B, F" << std::endl;
126
127 throw OomphLibError(error_message.str(),
128 OOMPH_CURRENT_FUNCTION,
129 OOMPH_EXCEPTION_LOCATION);
130 }
131
132 break;
133
134
135 // Macro element 1:right
136 case 1:
137
138 // Which direction?
139 if (idirect == L)
140 {
141 r_right_L(t, s, f);
142 }
143 else if (idirect == R)
144 {
145 r_right_R(t, s, f);
146 }
147 else if (idirect == D)
148 {
149 r_right_D(t, s, f);
150 }
151 else if (idirect == U)
152 {
153 r_right_U(t, s, f);
154 }
155 else if (idirect == B)
156 {
157 r_right_B(t, s, f);
158 }
159 else if (idirect == F)
160 {
161 r_right_F(t, s, f);
162 }
163 else
164 {
165 std::ostringstream error_message;
166 error_message << "idirect is " << OcTree::Direct_string[idirect]
167 << "not one of L, R, U, D, B, F" << std::endl;
168
169 throw OomphLibError(error_message.str(),
170 OOMPH_CURRENT_FUNCTION,
171 OOMPH_EXCEPTION_LOCATION);
172 }
173
174 break;
175
176 // Macro element 2: Up
177 case 2:
178
179 // Which direction?
180 if (idirect == L)
181 {
182 r_up_L(t, s, f);
183 }
184 else if (idirect == R)
185 {
186 r_up_R(t, s, f);
187 }
188 else if (idirect == D)
189 {
190 r_up_D(t, s, f);
191 }
192 else if (idirect == U)
193 {
194 r_up_U(t, s, f);
195 }
196 else if (idirect == B)
197 {
198 r_up_B(t, s, f);
199 }
200 else if (idirect == F)
201 {
202 r_up_F(t, s, f);
203 }
204 else
205 {
206 std::ostringstream error_message;
207 error_message << "idirect is " << OcTree::Direct_string[idirect]
208 << "not one of L, R, U, D, B, F" << std::endl;
209
210 throw OomphLibError(error_message.str(),
211 OOMPH_CURRENT_FUNCTION,
212 OOMPH_EXCEPTION_LOCATION);
213 }
214
215 break;
216
217 // Macro element 3: Front
218 case 3:
219 // Which direction?
220 if (idirect == L)
221 {
222 r_front_L(t, s, f);
223 }
224 else if (idirect == R)
225 {
226 r_front_R(t, s, f);
227 }
228 else if (idirect == D)
229 {
230 r_front_D(t, s, f);
231 }
232 else if (idirect == U)
233 {
234 r_front_U(t, s, f);
235 }
236 else if (idirect == B)
237 {
238 r_front_B(t, s, f);
239 }
240 else if (idirect == F)
241 {
242 r_front_F(t, s, f);
243 }
244 else
245 {
246 std::ostringstream error_message;
247 error_message << "idirect is " << OcTree::Direct_string[idirect]
248 << "not one of L, R, U, D, B, F" << std::endl;
249
250 throw OomphLibError(error_message.str(),
251 OOMPH_CURRENT_FUNCTION,
252 OOMPH_EXCEPTION_LOCATION);
253 }
254
255 break;
256
257 default:
258
259 std::ostringstream error_message;
260 error_message << "imacro is " << OcTree::Direct_string[idirect]
261 << ", but should be in the range 0-3" << std::endl;
262
263 throw OomphLibError(error_message.str(),
264 OOMPH_CURRENT_FUNCTION,
265 OOMPH_EXCEPTION_LOCATION);
266 }
267 }
268
269
270 private:
271 // Radius of the sphere
272 double Radius;
273
274 /// Boundary of central box macro element
275 /// zeta \f$ \in [-1,1]^2 \f$
276 void r_centr_L(const unsigned& t,
277 const Vector<double>& zeta,
278 Vector<double>& f);
279
280 /// Boundary of central box macro element
281 /// zeta \f$ \in [-1,1]^2 \f$
282 void r_centr_R(const unsigned& t,
283 const Vector<double>& zeta,
284 Vector<double>& f);
285
286 /// Boundary of central box macro element
287 /// zeta \f$ \in [-1,1]^2 \f$
288 void r_centr_D(const unsigned& t,
289 const Vector<double>& zeta,
290 Vector<double>& f);
291
292 /// Boundary of central box macro element
293 /// zeta \f$ \in [-1,1]^2 \f$
294 void r_centr_U(const unsigned& t,
295 const Vector<double>& zeta,
296 Vector<double>& f);
297
298 /// Boundary of central box macro element
299 /// zeta \f$ \in [-1,1]^2 \f$
300 void r_centr_B(const unsigned& t,
301 const Vector<double>& zeta,
302 Vector<double>& f);
303
304 /// Boundary of central box macro element
305 /// zeta \f$ \in [-1,1]^2 \f$
306 void r_centr_F(const unsigned& t,
307 const Vector<double>& zeta,
308 Vector<double>& f);
309
310 /// Boundary of right box macro element
311 /// zeta \f$ \in [-1,1]^2 \f$
312 void r_right_L(const unsigned& t,
313 const Vector<double>& zeta,
314 Vector<double>& f);
315
316 /// Boundary of right box macro element
317 /// zeta \f$ \in [-1,1]^2 \f$
318 void r_right_R(const unsigned& t,
319 const Vector<double>& zeta,
320 Vector<double>& f);
321
322 /// Boundary of right box macro element
323 /// zeta \f$ \in [-1,1]^2 \f$
324 void r_right_D(const unsigned& t,
325 const Vector<double>& zeta,
326 Vector<double>& f);
327
328 /// Boundary of right box macro element
329 /// zeta \f$ \in [-1,1]^2 \f$
330 void r_right_U(const unsigned& t,
331 const Vector<double>& zeta,
332 Vector<double>& f);
333
334 /// Boundary of right box macro element
335 /// zeta \f$ \in [-1,1]^2 \f$
336 void r_right_B(const unsigned& t,
337 const Vector<double>& zeta,
338 Vector<double>& f);
339
340 /// Boundary of right box macro element
341 /// zeta \f$ \in [-1,1]^2 \f$
342 void r_right_F(const unsigned& t,
343 const Vector<double>& zeta,
344 Vector<double>& f);
345
346 /// Boundary of top left box macro element
347 /// zeta \f$ \in [-1,1]^2 \f$
348 void r_up_L(const unsigned& t,
349 const Vector<double>& zeta,
350 Vector<double>& f);
351
352 /// Boundary of top left box macro element
353 /// zeta \f$ \in [-1,1]^2 \f$
354 void r_up_R(const unsigned& t,
355 const Vector<double>& zeta,
356 Vector<double>& f);
357
358 /// Boundary of top left box macro element
359 /// zeta \f$ \in [-1,1]^2 \f$
360 void r_up_D(const unsigned& t,
361 const Vector<double>& zeta,
362 Vector<double>& f);
363
364 /// Boundary of top left box macro element
365 /// zeta \f$ \in [-1,1]^2 \f$
366 void r_up_U(const unsigned& t,
367 const Vector<double>& zeta,
368 Vector<double>& f);
369
370 /// Boundary of top left box macro element
371 /// zeta \f$ \in [-1,1]^2 \f$
372 void r_up_B(const unsigned& t,
373 const Vector<double>& zeta,
374 Vector<double>& f);
375
376 /// Boundary of top left box macro element
377 /// zeta \f$ \in [-1,1]^2 \f$
378 void r_up_F(const unsigned& t,
379 const Vector<double>& zeta,
380 Vector<double>& f);
381
382 /// Boundary of top left box macro element
383 /// zeta \f$ \in [-1,1]^2 \f$
384 void r_front_L(const unsigned& t,
385 const Vector<double>& zeta,
386 Vector<double>& f);
387
388 /// Boundary of top left box macro element
389 /// zeta \f$ \in [-1,1]^2 \f$
390 void r_front_R(const unsigned& t,
391 const Vector<double>& zeta,
392 Vector<double>& f);
393
394 /// Boundary of top left box macro element
395 /// zeta \f$ \in [-1,1]^2 \f$
396 void r_front_D(const unsigned& t,
397 const Vector<double>& zeta,
398 Vector<double>& f);
399
400 /// Boundary of top left box macro element
401 /// zeta \f$ \in [-1,1]^2 \f$
402 void r_front_U(const unsigned& t,
403 const Vector<double>& zeta,
404 Vector<double>& f);
405
406 /// Boundary of top left box macro element
407 /// zeta \f$ \in [-1,1]^2 \f$
408 void r_front_B(const unsigned& t,
409 const Vector<double>& zeta,
410 Vector<double>& f);
411
412 /// Boundary of top left box macro element
413 /// zeta \f$ \in [-1,1]^2 \f$
414 void r_front_F(const unsigned& t,
415 const Vector<double>& zeta,
416 Vector<double>& f);
417 };
418
419 /// /////////////////////////////////////////////////////////////////////
420 /// /////////////////////////////////////////////////////////////////////
421 /// /////////////////////////////////////////////////////////////////////
422
423
424 //=================================================================
425 /// Boundary of central box macro element
426 /// zeta \f$ \in [-1,1]^2 \f$
427 //=================================================================
428 void EighthSphereDomain::r_centr_L(const unsigned& t,
429 const Vector<double>& zeta,
430 Vector<double>& f)
431 {
432 f[0] = 0;
433 f[1] = Radius * 0.25 * (1.0 + zeta[0]);
434 f[2] = Radius * 0.25 * (1.0 + zeta[1]);
435 }
436
437
438 //=================================================================
439 /// Boundary of central box macro element
440 /// zeta \f$ \in [-1,1]^2 \f$
441 //=================================================================
442 void EighthSphereDomain::r_centr_R(const unsigned& t,
443 const Vector<double>& zeta,
444 Vector<double>& f)
445 {
446 f[0] = Radius * 0.5;
447 f[1] = Radius * 0.25 * (1.0 + zeta[0]);
448 f[2] = Radius * 0.25 * (1.0 + zeta[1]);
449 }
450
451
452 //=================================================================
453 /// Boundary of central box macro element
454 /// zeta \f$ \in [-1,1]^2 \f$
455 //=================================================================
456 void EighthSphereDomain::r_centr_D(const unsigned& t,
457 const Vector<double>& zeta,
458 Vector<double>& f)
459 {
460 f[0] = Radius * 0.25 * (1.0 + zeta[0]);
461 f[1] = 0;
462 f[2] = Radius * 0.25 * (1.0 + zeta[1]);
463 }
464
465
466 //=================================================================
467 /// Boundary of central box macro element
468 /// zeta \f$ \in [-1,1]^2 \f$
469 //=================================================================
470 void EighthSphereDomain::r_centr_U(const unsigned& t,
471 const Vector<double>& zeta,
472 Vector<double>& f)
473 {
474 f[0] = Radius * 0.25 * (1.0 + zeta[0]);
475 f[1] = Radius * 0.5;
476 f[2] = Radius * 0.25 * (1.0 + zeta[1]);
477 }
478
479
480 //=================================================================
481 /// Boundary of central box macro element
482 /// zeta \f$ \in [-1,1]^2 \f$
483 //=================================================================
484 void EighthSphereDomain::r_centr_B(const unsigned& t,
485 const Vector<double>& zeta,
486 Vector<double>& f)
487 {
488 f[0] = Radius * 0.25 * (1.0 + zeta[0]);
489 f[1] = Radius * 0.25 * (1.0 + zeta[1]);
490 f[2] = 0.0;
491 }
492
493 //=================================================================
494 /// Boundary of central box macro element
495 /// zeta \f$ \in [-1,1]^2 \f$
496 //=================================================================
497 void EighthSphereDomain::r_centr_F(const unsigned& t,
498 const Vector<double>& zeta,
499 Vector<double>& f)
500 {
501 f[0] = Radius * 0.25 * (1.0 + zeta[0]);
502 f[1] = Radius * 0.25 * (1.0 + zeta[1]);
503 f[2] = Radius * 0.5;
504 }
505
506
507 //=================================================================
508 /// Boundary of right box macro element
509 /// zeta \f$ \in [-1,1]^2 \f$
510 //=================================================================
511 void EighthSphereDomain::r_right_L(const unsigned& t,
512 const Vector<double>& zeta,
513 Vector<double>& f)
514 {
515 r_centr_R(t, zeta, f);
516 }
517
518
519 //=================================================================
520 /// Boundary of right box macro element
521 /// zeta \f$ \in [-1,1]^2 \f$
522 //=================================================================
523 void EighthSphereDomain::r_right_R(const unsigned& t,
524 const Vector<double>& zeta,
525 Vector<double>& f)
526 {
527 double k0 = 0.5 * (1.0 + zeta[0]);
528 double k1 = 0.5 * (1.0 + zeta[1]);
529 Vector<double> p(3);
530 Vector<double> point1(3);
531 Vector<double> point2(3);
532
533 point1[0] = Radius - 0.5 * Radius * k0;
534 point1[1] = 0.5 * Radius * k0;
535 point1[2] = 0.0;
536 point2[0] = 0.5 * Radius - k0 * Radius / 6.0;
537 point2[1] = k0 * Radius / 3.0;
538 point2[2] = 0.5 * Radius - k0 * Radius / 6.0;
539
540 for (unsigned i = 0; i < 3; i++)
541 {
542 p[i] = point1[i] + k1 * (point2[i] - point1[i]);
543 }
544 double alpha = Radius / std::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
545
546 for (unsigned i = 0; i < 3; i++)
547 {
548 f[i] = alpha * p[i];
549 }
550 }
551
552 //=================================================================
553 /// Boundary of right box macro element
554 /// zeta \f$ \in [-1,1]^2 \f$
555 //=================================================================
556 void EighthSphereDomain::r_right_D(const unsigned& t,
557 const Vector<double>& zeta,
558 Vector<double>& f)
559 {
560 // position vector on sphere
561 Vector<double> on_sphere(3);
562 Vector<double> temp_zeta(2);
563 temp_zeta[0] = -1.0;
564 temp_zeta[1] = zeta[1];
565 r_right_R(t, temp_zeta, on_sphere);
566
567 // position vector on center box
568 Vector<double> on_center(3);
569 on_center[0] = 0.5 * Radius;
570 on_center[1] = 0.0;
571 on_center[2] = 0.5 * Radius * 0.5 * (zeta[1] + 1.0);
572 // strait line across
573 for (unsigned i = 0; i < 3; i++)
574 {
575 f[i] =
576 on_center[i] + 0.5 * (zeta[0] + 1.0) * (on_sphere[i] - on_center[i]);
577 }
578 }
579
580
581 //=================================================================
582 /// Boundary of right box macro element
583 /// zeta \f$ \in [-1,1]^2 \f$
584 //=================================================================
585 void EighthSphereDomain::r_right_U(const unsigned& t,
586 const Vector<double>& zeta,
587 Vector<double>& f)
588 {
589 // position vector on sphere
590 Vector<double> on_sphere(3);
591 Vector<double> temp_zeta(2);
592 temp_zeta[0] = 1.0;
593 temp_zeta[1] = zeta[1];
594 r_right_R(t, temp_zeta, on_sphere);
595
596 // position vector on center box
597 Vector<double> on_center(3);
598 on_center[0] = 0.5 * Radius;
599 on_center[1] = 0.5 * Radius;
600 on_center[2] = 0.5 * Radius * 0.5 * (zeta[1] + 1.0);
601 // strait line across
602 for (unsigned i = 0; i < 3; i++)
603 {
604 f[i] =
605 on_center[i] + 0.5 * (zeta[0] + 1.0) * (on_sphere[i] - on_center[i]);
606 }
607 }
608
609 //=================================================================
610 /// Boundary of right box macro element
611 /// zeta \f$ \in [-1,1]^2 \f$
612 //=================================================================
613 void EighthSphereDomain::r_right_B(const unsigned& t,
614 const Vector<double>& zeta,
615 Vector<double>& f)
616 {
617 // position vector on sphere
618 Vector<double> on_sphere(3);
619 Vector<double> temp_zeta(2);
620 temp_zeta[0] = zeta[1];
621 temp_zeta[1] = -1.0;
622 r_right_R(t, temp_zeta, on_sphere);
623
624 // position vector on center box
625 Vector<double> on_center(3);
626 on_center[0] = 0.5 * Radius;
627 on_center[1] = 0.5 * Radius * 0.5 * (zeta[1] + 1.0);
628 on_center[2] = 0.0;
629 // strait line across
630 for (unsigned i = 0; i < 3; i++)
631 {
632 f[i] =
633 on_center[i] + 0.5 * (zeta[0] + 1.0) * (on_sphere[i] - on_center[i]);
634 }
635 }
636
637
638 //=================================================================
639 /// Boundary of right box macro element
640 /// zeta \f$ \in [-1,1]^2 \f$
641 //=================================================================
642 void EighthSphereDomain::r_right_F(const unsigned& t,
643 const Vector<double>& zeta,
644 Vector<double>& f)
645 {
646 // position vector on sphere
647 Vector<double> on_sphere(3);
648 Vector<double> temp_zeta(2);
649 temp_zeta[0] = zeta[1];
650 temp_zeta[1] = 1.0;
651 r_right_R(t, temp_zeta, on_sphere);
652
653 // position vector on center box
654 Vector<double> on_center(3);
655 on_center[0] = 0.5 * Radius;
656 on_center[1] = 0.5 * Radius * 0.5 * (zeta[1] + 1.0);
657 on_center[2] = 0.5 * Radius;
658 // strait line across
659 for (unsigned i = 0; i < 3; i++)
660 {
661 f[i] =
662 on_center[i] + 0.5 * (zeta[0] + 1.0) * (on_sphere[i] - on_center[i]);
663 }
664 }
665
666 //=================================================================
667 /// Boundary of top left box macro element
668 /// zeta \f$ \in [-1,1]^2 \f$
669 //=================================================================
670 void EighthSphereDomain::r_up_L(const unsigned& t,
671 const Vector<double>& zeta,
672 Vector<double>& f)
673 {
674 // position vector on sphere
675 Vector<double> on_sphere(3);
676 Vector<double> temp_zeta(2);
677 temp_zeta[0] = -1.0;
678 temp_zeta[1] = zeta[1];
679 r_up_U(t, temp_zeta, on_sphere);
680
681 // position vector on center box
682 Vector<double> on_center(3);
683 on_center[0] = 0.0;
684 on_center[1] = 0.5 * Radius;
685 on_center[2] = 0.5 * Radius * 0.5 * (zeta[1] + 1.0);
686 // strait line across
687 for (unsigned i = 0; i < 3; i++)
688 {
689 f[i] =
690 on_center[i] + 0.5 * (zeta[0] + 1.0) * (on_sphere[i] - on_center[i]);
691 }
692 }
693
694
695 //=================================================================
696 /// Boundary of top left box macro element
697 /// zeta \f$ \in [-1,1]^2 \f$
698 //=================================================================
699 void EighthSphereDomain::r_up_R(const unsigned& t,
700 const Vector<double>& zeta,
701 Vector<double>& f)
702 {
703 r_right_U(t, zeta, f);
704 }
705
706
707 //=================================================================
708 /// Boundary of top left box macro element
709 /// zeta \f$ \in [-1,1]^2 \f$
710 //=================================================================
711 void EighthSphereDomain::r_up_D(const unsigned& t,
712 const Vector<double>& zeta,
713 Vector<double>& f)
714 {
715 r_centr_U(t, zeta, f);
716 }
717
718 //=================================================================
719 /// Boundary of top left box macro element
720 /// zeta \f$ \in [-1,1]^2 \f$
721 //=================================================================
722 void EighthSphereDomain::r_up_U(const unsigned& t,
723 const Vector<double>& zeta,
724 Vector<double>& f)
725 {
726 double k0 = 0.5 * (1.0 + zeta[0]);
727 double k1 = 0.5 * (1.0 + zeta[1]);
728 Vector<double> p(3);
729 Vector<double> point1(3);
730 Vector<double> point2(3);
731
732 point1[0] = 0.5 * Radius * k0;
733 point1[1] = Radius - 0.5 * Radius * k0;
734 point1[2] = 0;
735 point2[0] = k0 * Radius / 3.0;
736 point2[1] = 0.5 * Radius - k0 * Radius / 6.0;
737 point2[2] = 0.5 * Radius - k0 * Radius / 6.0;
738
739 for (unsigned i = 0; i < 3; i++)
740 {
741 p[i] = point1[i] + k1 * (point2[i] - point1[i]);
742 }
743 double alpha = Radius / std::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
744
745 for (unsigned i = 0; i < 3; i++)
746 {
747 f[i] = alpha * p[i];
748 }
749 }
750
751 //=================================================================
752 /// Boundary of top left box macro element
753 /// zeta \f$ \in [-1,1]^2 \f$
754 //=================================================================
755 void EighthSphereDomain::r_up_B(const unsigned& t,
756 const Vector<double>& zeta,
757 Vector<double>& f)
758 {
759 // position vector on sphere
760 Vector<double> on_sphere(3);
761 Vector<double> temp_zeta(2);
762 temp_zeta[0] = zeta[0];
763 temp_zeta[1] = -1.0;
764 r_up_U(t, temp_zeta, on_sphere);
765
766 // position vector on center box
767 Vector<double> on_center(3);
768 on_center[0] = 0.5 * Radius * 0.5 * (zeta[0] + 1.0);
769 on_center[1] = 0.5 * Radius;
770 on_center[2] = 0.0;
771 // strait line across
772 for (unsigned i = 0; i < 3; i++)
773 {
774 f[i] =
775 on_center[i] + 0.5 * (zeta[1] + 1.0) * (on_sphere[i] - on_center[i]);
776 }
777 }
778
779
780 //=================================================================
781 /// Boundary of top left box macro element
782 /// zeta \f$ \in [-1,1]^2 \f$
783 //=================================================================
784 void EighthSphereDomain::r_up_F(const unsigned& t,
785 const Vector<double>& zeta,
786 Vector<double>& f)
787 {
788 // position vector on sphere
789 Vector<double> on_sphere(3);
790 Vector<double> temp_zeta(2);
791 temp_zeta[0] = zeta[0];
792 temp_zeta[1] = 1.0;
793 r_up_U(t, temp_zeta, on_sphere);
794
795 // position vector on center box
796 Vector<double> on_center(3);
797 on_center[0] = 0.5 * Radius * 0.5 * (zeta[0] + 1.0);
798 on_center[1] = 0.5 * Radius;
799 on_center[2] = 0.5 * Radius;
800 // straight line across
801 for (unsigned i = 0; i < 3; i++)
802 {
803 f[i] =
804 on_center[i] + 0.5 * (zeta[1] + 1.0) * (on_sphere[i] - on_center[i]);
805 }
806 }
807
808
809 //=================================================================
810 /// Boundary of top left box macro element
811 /// zeta \f$ \in [-1,1]^2 \f$
812 //=================================================================
813 void EighthSphereDomain::r_front_L(const unsigned& t,
814 const Vector<double>& zeta,
815 Vector<double>& f)
816 {
817 // position vector on sphere
818 Vector<double> on_sphere(3);
819 Vector<double> temp_zeta(2);
820 temp_zeta[0] = -1.0;
821 temp_zeta[1] = zeta[0];
822 r_front_F(t, temp_zeta, on_sphere);
823
824 // position vector on center box
825 Vector<double> on_center(3);
826 on_center[0] = 0.0;
827 on_center[1] = 0.5 * Radius * 0.5 * (zeta[0] + 1.0);
828 on_center[2] = 0.5 * Radius;
829 // straight line across
830 for (unsigned i = 0; i < 3; i++)
831 {
832 f[i] =
833 on_center[i] + 0.5 * (zeta[1] + 1.0) * (on_sphere[i] - on_center[i]);
834 }
835 }
836
837
838 //=================================================================
839 /// Boundary of top left box macro element
840 /// zeta \f$ \in [-1,1]^2 \f$
841 //=================================================================
842 void EighthSphereDomain::r_front_R(const unsigned& t,
843 const Vector<double>& zeta,
844 Vector<double>& f)
845 {
846 Vector<double> zeta2(2);
847 zeta2[0] = zeta[1];
848 zeta2[1] = zeta[0];
849 r_right_F(t, zeta2, f);
850 }
851
852
853 //=================================================================
854 /// Boundary of top left box macro element
855 /// zeta \f$ \in [-1,1]^2 \f$
856 //=================================================================
857 void EighthSphereDomain::r_front_D(const unsigned& t,
858 const Vector<double>& zeta,
859 Vector<double>& f)
860 {
861 // position vector on sphere
862 Vector<double> on_sphere(3);
863 Vector<double> temp_zeta(2);
864 temp_zeta[0] = zeta[0];
865 temp_zeta[1] = -1.0;
866 r_front_F(t, temp_zeta, on_sphere);
867
868 // position vector on center box
869 Vector<double> on_center(3);
870 on_center[0] = 0.5 * Radius * 0.5 * (zeta[0] + 1.0);
871 on_center[1] = 0.0;
872 on_center[2] = 0.5 * Radius;
873 // straight line across
874 for (unsigned i = 0; i < 3; i++)
875 {
876 f[i] =
877 on_center[i] + 0.5 * (zeta[1] + 1.0) * (on_sphere[i] - on_center[i]);
878 }
879 }
880
881
882 //=================================================================
883 /// Boundary of top left box macro element
884 /// zeta \f$ \in [-1,1]^2 \f$
885 //=================================================================
886 void EighthSphereDomain::r_front_U(const unsigned& t,
887 const Vector<double>& zeta,
888 Vector<double>& f)
889 {
890 r_up_F(t, zeta, f);
891 }
892
893
894 //=================================================================
895 /// Boundary of top left box macro element
896 /// zeta \f$ \in [-1,1]^2 \f$
897 //=================================================================
898 void EighthSphereDomain::r_front_B(const unsigned& t,
899 const Vector<double>& zeta,
900 Vector<double>& f)
901 {
902 r_centr_F(t, zeta, f);
903 }
904
905
906 //=================================================================
907 /// Boundary of top left box macro element
908 /// zeta \f$ \in [-1,1]^2 \f$
909 //=================================================================
910 void EighthSphereDomain::r_front_F(const unsigned& t,
911 const Vector<double>& zeta,
912 Vector<double>& f)
913 {
914 double k0 = 0.5 * (1.0 + zeta[0]);
915 double k1 = 0.5 * (1.0 + zeta[1]);
916 Vector<double> p(3);
917 Vector<double> point1(3);
918 Vector<double> point2(3);
919
920 point1[0] = 0.5 * Radius * k0;
921 point1[1] = 0.0;
922 point1[2] = Radius - k0 * 0.5 * Radius;
923 point2[0] = k0 * Radius / 3.0;
924 point2[1] = 0.5 * Radius - k0 * Radius / 6.0;
925 point2[2] = 0.5 * Radius - k0 * Radius / 6.0;
926
927 for (unsigned i = 0; i < 3; i++)
928 {
929 p[i] = point1[i] + k1 * (point2[i] - point1[i]);
930 }
931 double alpha = Radius / std::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
932
933 for (unsigned i = 0; i < 3; i++)
934 {
935 f[i] = alpha * p[i];
936 }
937 }
938
939
940} // namespace oomph
941
942#endif
oomph::EighthSphereDomain
Eighth sphere as domain. Domain is parametrised by four macro elements.
Definition eighth_sphere_domain.h:42
oomph::EighthSphereDomain::r_right_U
void r_right_U(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
Definition eighth_sphere_domain.h:585
oomph::EighthSphereDomain::r_centr_F
void r_centr_F(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
Definition eighth_sphere_domain.h:497
oomph::EighthSphereDomain::r_front_L
void r_front_L(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:813
oomph::EighthSphereDomain::r_centr_R
void r_centr_R(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
Definition eighth_sphere_domain.h:442
oomph::EighthSphereDomain::macro_element_boundary
void macro_element_boundary(const unsigned &t, const unsigned &imacro, const unsigned &idirect, const Vector< double > &s, Vector< double > &f)
Vector representation of the imacro-th macro element boundary idirect (L/R/D/U/B/F) at time level t (...
Definition eighth_sphere_domain.h:74
oomph::EighthSphereDomain::r_up_U
void r_up_U(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:722
oomph::EighthSphereDomain::r_centr_B
void r_centr_B(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
Definition eighth_sphere_domain.h:484
oomph::EighthSphereDomain::r_front_D
void r_front_D(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:857
oomph::EighthSphereDomain::r_right_R
void r_right_R(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
Definition eighth_sphere_domain.h:523
oomph::EighthSphereDomain::r_front_R
void r_front_R(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:842
oomph::EighthSphereDomain::EighthSphereDomain
EighthSphereDomain(const EighthSphereDomain &)=delete
Broken copy constructor.
oomph::EighthSphereDomain::r_front_F
void r_front_F(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:910
oomph::EighthSphereDomain::r_right_D
void r_right_D(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
Definition eighth_sphere_domain.h:556
oomph::EighthSphereDomain::EighthSphereDomain
EighthSphereDomain(const double &radius)
Constructor: Pass the radius of the sphere.
Definition eighth_sphere_domain.h:45
oomph::EighthSphereDomain::operator=
void operator=(const EighthSphereDomain &)=delete
Broken assignment operator.
oomph::EighthSphereDomain::r_up_F
void r_up_F(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:784
oomph::EighthSphereDomain::r_right_F
void r_right_F(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
Definition eighth_sphere_domain.h:642
oomph::EighthSphereDomain::r_front_B
void r_front_B(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:898
oomph::EighthSphereDomain::r_centr_D
void r_centr_D(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
Definition eighth_sphere_domain.h:456
oomph::EighthSphereDomain::r_up_B
void r_up_B(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:755
oomph::EighthSphereDomain::~EighthSphereDomain
~EighthSphereDomain()
Destructor: Empty; cleanup done in base class.
Definition eighth_sphere_domain.h:67
oomph::EighthSphereDomain::r_right_B
void r_right_B(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
Definition eighth_sphere_domain.h:613
oomph::EighthSphereDomain::r_right_L
void r_right_L(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
Definition eighth_sphere_domain.h:511
oomph::EighthSphereDomain::r_up_R
void r_up_R(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:699
oomph::EighthSphereDomain::r_up_D
void r_up_D(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:711
oomph::EighthSphereDomain::r_centr_L
void r_centr_L(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
Definition eighth_sphere_domain.h:428
oomph::EighthSphereDomain::r_front_U
void r_front_U(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:886
oomph::EighthSphereDomain::r_up_L
void r_up_L(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
Definition eighth_sphere_domain.h:670
oomph::EighthSphereDomain::r_centr_U
void r_centr_U(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
Definition eighth_sphere_domain.h:470
oomph::RefineableTriangleMesh
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Definition triangle_mesh.template.h:2249
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition annular_domain.h:35