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-2025 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 // Macro element 1:right
135 case 1:
136
137 // Which direction?
138 if (idirect == L)
139 {
140 r_right_L(t, s, f);
141 }
142 else if (idirect == R)
143 {
144 r_right_R(t, s, f);
145 }
146 else if (idirect == D)
147 {
148 r_right_D(t, s, f);
149 }
150 else if (idirect == U)
151 {
152 r_right_U(t, s, f);
153 }
154 else if (idirect == B)
155 {
156 r_right_B(t, s, f);
157 }
158 else if (idirect == F)
159 {
160 r_right_F(t, s, f);
161 }
162 else
163 {
164 std::ostringstream error_message;
165 error_message << "idirect is " << OcTree::Direct_string[idirect]
166 << "not one of L, R, U, D, B, F" << std::endl;
167
168 throw OomphLibError(error_message.str(),
169 OOMPH_CURRENT_FUNCTION,
170 OOMPH_EXCEPTION_LOCATION);
171 }
172
173 break;
174
175 // Macro element 2: Up
176 case 2:
177
178 // Which direction?
179 if (idirect == L)
180 {
181 r_up_L(t, s, f);
182 }
183 else if (idirect == R)
184 {
185 r_up_R(t, s, f);
186 }
187 else if (idirect == D)
188 {
189 r_up_D(t, s, f);
190 }
191 else if (idirect == U)
192 {
193 r_up_U(t, s, f);
194 }
195 else if (idirect == B)
196 {
197 r_up_B(t, s, f);
198 }
199 else if (idirect == F)
200 {
201 r_up_F(t, s, f);
202 }
203 else
204 {
205 std::ostringstream error_message;
206 error_message << "idirect is " << OcTree::Direct_string[idirect]
207 << "not one of L, R, U, D, B, F" << std::endl;
208
209 throw OomphLibError(error_message.str(),
210 OOMPH_CURRENT_FUNCTION,
211 OOMPH_EXCEPTION_LOCATION);
212 }
213
214 break;
215
216 // Macro element 3: Front
217 case 3:
218 // Which direction?
219 if (idirect == L)
220 {
221 r_front_L(t, s, f);
222 }
223 else if (idirect == R)
224 {
225 r_front_R(t, s, f);
226 }
227 else if (idirect == D)
228 {
229 r_front_D(t, s, f);
230 }
231 else if (idirect == U)
232 {
233 r_front_U(t, s, f);
234 }
235 else if (idirect == B)
236 {
237 r_front_B(t, s, f);
238 }
239 else if (idirect == F)
240 {
241 r_front_F(t, s, f);
242 }
243 else
244 {
245 std::ostringstream error_message;
246 error_message << "idirect is " << OcTree::Direct_string[idirect]
247 << "not one of L, R, U, D, B, F" << std::endl;
248
249 throw OomphLibError(error_message.str(),
250 OOMPH_CURRENT_FUNCTION,
251 OOMPH_EXCEPTION_LOCATION);
252 }
253
254 break;
255
256 default:
257
258 std::ostringstream error_message;
259 error_message << "imacro is " << OcTree::Direct_string[idirect]
260 << ", but should be in the range 0-3" << std::endl;
261
262 throw OomphLibError(error_message.str(),
263 OOMPH_CURRENT_FUNCTION,
264 OOMPH_EXCEPTION_LOCATION);
265 }
266 }
267
268 private:
269 // Radius of the sphere
270 double Radius;
271
272 /// Boundary of central box macro element
273 /// zeta \f$ \in [-1,1]^2 \f$
274 void r_centr_L(const unsigned& t,
275 const Vector<double>& zeta,
276 Vector<double>& f);
277
278 /// Boundary of central box macro element
279 /// zeta \f$ \in [-1,1]^2 \f$
280 void r_centr_R(const unsigned& t,
281 const Vector<double>& zeta,
282 Vector<double>& f);
283
284 /// Boundary of central box macro element
285 /// zeta \f$ \in [-1,1]^2 \f$
286 void r_centr_D(const unsigned& t,
287 const Vector<double>& zeta,
288 Vector<double>& f);
289
290 /// Boundary of central box macro element
291 /// zeta \f$ \in [-1,1]^2 \f$
292 void r_centr_U(const unsigned& t,
293 const Vector<double>& zeta,
294 Vector<double>& f);
295
296 /// Boundary of central box macro element
297 /// zeta \f$ \in [-1,1]^2 \f$
298 void r_centr_B(const unsigned& t,
299 const Vector<double>& zeta,
300 Vector<double>& f);
301
302 /// Boundary of central box macro element
303 /// zeta \f$ \in [-1,1]^2 \f$
304 void r_centr_F(const unsigned& t,
305 const Vector<double>& zeta,
306 Vector<double>& f);
307
308 /// Boundary of right box macro element
309 /// zeta \f$ \in [-1,1]^2 \f$
310 void r_right_L(const unsigned& t,
311 const Vector<double>& zeta,
312 Vector<double>& f);
313
314 /// Boundary of right box macro element
315 /// zeta \f$ \in [-1,1]^2 \f$
316 void r_right_R(const unsigned& t,
317 const Vector<double>& zeta,
318 Vector<double>& f);
319
320 /// Boundary of right box macro element
321 /// zeta \f$ \in [-1,1]^2 \f$
322 void r_right_D(const unsigned& t,
323 const Vector<double>& zeta,
324 Vector<double>& f);
325
326 /// Boundary of right box macro element
327 /// zeta \f$ \in [-1,1]^2 \f$
328 void r_right_U(const unsigned& t,
329 const Vector<double>& zeta,
330 Vector<double>& f);
331
332 /// Boundary of right box macro element
333 /// zeta \f$ \in [-1,1]^2 \f$
334 void r_right_B(const unsigned& t,
335 const Vector<double>& zeta,
336 Vector<double>& f);
337
338 /// Boundary of right box macro element
339 /// zeta \f$ \in [-1,1]^2 \f$
340 void r_right_F(const unsigned& t,
341 const Vector<double>& zeta,
342 Vector<double>& f);
343
344 /// Boundary of top left box macro element
345 /// zeta \f$ \in [-1,1]^2 \f$
346 void r_up_L(const unsigned& t,
347 const Vector<double>& zeta,
348 Vector<double>& f);
349
350 /// Boundary of top left box macro element
351 /// zeta \f$ \in [-1,1]^2 \f$
352 void r_up_R(const unsigned& t,
353 const Vector<double>& zeta,
354 Vector<double>& f);
355
356 /// Boundary of top left box macro element
357 /// zeta \f$ \in [-1,1]^2 \f$
358 void r_up_D(const unsigned& t,
359 const Vector<double>& zeta,
360 Vector<double>& f);
361
362 /// Boundary of top left box macro element
363 /// zeta \f$ \in [-1,1]^2 \f$
364 void r_up_U(const unsigned& t,
365 const Vector<double>& zeta,
366 Vector<double>& f);
367
368 /// Boundary of top left box macro element
369 /// zeta \f$ \in [-1,1]^2 \f$
370 void r_up_B(const unsigned& t,
371 const Vector<double>& zeta,
372 Vector<double>& f);
373
374 /// Boundary of top left box macro element
375 /// zeta \f$ \in [-1,1]^2 \f$
376 void r_up_F(const unsigned& t,
377 const Vector<double>& zeta,
378 Vector<double>& f);
379
380 /// Boundary of top left box macro element
381 /// zeta \f$ \in [-1,1]^2 \f$
382 void r_front_L(const unsigned& t,
383 const Vector<double>& zeta,
384 Vector<double>& f);
385
386 /// Boundary of top left box macro element
387 /// zeta \f$ \in [-1,1]^2 \f$
388 void r_front_R(const unsigned& t,
389 const Vector<double>& zeta,
390 Vector<double>& f);
391
392 /// Boundary of top left box macro element
393 /// zeta \f$ \in [-1,1]^2 \f$
394 void r_front_D(const unsigned& t,
395 const Vector<double>& zeta,
396 Vector<double>& f);
397
398 /// Boundary of top left box macro element
399 /// zeta \f$ \in [-1,1]^2 \f$
400 void r_front_U(const unsigned& t,
401 const Vector<double>& zeta,
402 Vector<double>& f);
403
404 /// Boundary of top left box macro element
405 /// zeta \f$ \in [-1,1]^2 \f$
406 void r_front_B(const unsigned& t,
407 const Vector<double>& zeta,
408 Vector<double>& f);
409
410 /// Boundary of top left box macro element
411 /// zeta \f$ \in [-1,1]^2 \f$
412 void r_front_F(const unsigned& t,
413 const Vector<double>& zeta,
414 Vector<double>& f);
415 };
416
417 ////////////////////////////////////////////////////////////////////////
418 ////////////////////////////////////////////////////////////////////////
419 ////////////////////////////////////////////////////////////////////////
420
421 //=================================================================
422 /// Boundary of central box macro element
423 /// zeta \f$ \in [-1,1]^2 \f$
424 //=================================================================
425 void EighthSphereDomain::r_centr_L(const unsigned& t,
426 const Vector<double>& zeta,
427 Vector<double>& f)
428 {
429 f[0] = 0;
430 f[1] = Radius * 0.25 * (1.0 + zeta[0]);
431 f[2] = Radius * 0.25 * (1.0 + zeta[1]);
432 }
433
434 //=================================================================
435 /// Boundary of central box macro element
436 /// zeta \f$ \in [-1,1]^2 \f$
437 //=================================================================
438 void EighthSphereDomain::r_centr_R(const unsigned& t,
439 const Vector<double>& zeta,
440 Vector<double>& f)
441 {
442 f[0] = Radius * 0.5;
443 f[1] = Radius * 0.25 * (1.0 + zeta[0]);
444 f[2] = Radius * 0.25 * (1.0 + zeta[1]);
445 }
446
447 //=================================================================
448 /// Boundary of central box macro element
449 /// zeta \f$ \in [-1,1]^2 \f$
450 //=================================================================
451 void EighthSphereDomain::r_centr_D(const unsigned& t,
452 const Vector<double>& zeta,
453 Vector<double>& f)
454 {
455 f[0] = Radius * 0.25 * (1.0 + zeta[0]);
456 f[1] = 0;
457 f[2] = Radius * 0.25 * (1.0 + zeta[1]);
458 }
459
460 //=================================================================
461 /// Boundary of central box macro element
462 /// zeta \f$ \in [-1,1]^2 \f$
463 //=================================================================
464 void EighthSphereDomain::r_centr_U(const unsigned& t,
465 const Vector<double>& zeta,
466 Vector<double>& f)
467 {
468 f[0] = Radius * 0.25 * (1.0 + zeta[0]);
469 f[1] = Radius * 0.5;
470 f[2] = Radius * 0.25 * (1.0 + zeta[1]);
471 }
472
473 //=================================================================
474 /// Boundary of central box macro element
475 /// zeta \f$ \in [-1,1]^2 \f$
476 //=================================================================
477 void EighthSphereDomain::r_centr_B(const unsigned& t,
478 const Vector<double>& zeta,
479 Vector<double>& f)
480 {
481 f[0] = Radius * 0.25 * (1.0 + zeta[0]);
482 f[1] = Radius * 0.25 * (1.0 + zeta[1]);
483 f[2] = 0.0;
484 }
485
486 //=================================================================
487 /// Boundary of central box macro element
488 /// zeta \f$ \in [-1,1]^2 \f$
489 //=================================================================
490 void EighthSphereDomain::r_centr_F(const unsigned& t,
491 const Vector<double>& zeta,
492 Vector<double>& f)
493 {
494 f[0] = Radius * 0.25 * (1.0 + zeta[0]);
495 f[1] = Radius * 0.25 * (1.0 + zeta[1]);
496 f[2] = Radius * 0.5;
497 }
498
499 //=================================================================
500 /// Boundary of right box macro element
501 /// zeta \f$ \in [-1,1]^2 \f$
502 //=================================================================
503 void EighthSphereDomain::r_right_L(const unsigned& t,
504 const Vector<double>& zeta,
505 Vector<double>& f)
506 {
507 r_centr_R(t, zeta, f);
508 }
509
510 //=================================================================
511 /// Boundary of right box macro element
512 /// zeta \f$ \in [-1,1]^2 \f$
513 //=================================================================
514 void EighthSphereDomain::r_right_R(const unsigned& t,
515 const Vector<double>& zeta,
516 Vector<double>& f)
517 {
518 double k0 = 0.5 * (1.0 + zeta[0]);
519 double k1 = 0.5 * (1.0 + zeta[1]);
520 Vector<double> p(3);
521 Vector<double> point1(3);
522 Vector<double> point2(3);
523
524 point1[0] = Radius - 0.5 * Radius * k0;
525 point1[1] = 0.5 * Radius * k0;
526 point1[2] = 0.0;
527 point2[0] = 0.5 * Radius - k0 * Radius / 6.0;
528 point2[1] = k0 * Radius / 3.0;
529 point2[2] = 0.5 * Radius - k0 * Radius / 6.0;
530
531 for (unsigned i = 0; i < 3; i++)
532 {
533 p[i] = point1[i] + k1 * (point2[i] - point1[i]);
534 }
535 double alpha = Radius / std::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
536
537 for (unsigned i = 0; i < 3; i++)
538 {
539 f[i] = alpha * p[i];
540 }
541 }
542
543 //=================================================================
544 /// Boundary of right box macro element
545 /// zeta \f$ \in [-1,1]^2 \f$
546 //=================================================================
547 void EighthSphereDomain::r_right_D(const unsigned& t,
548 const Vector<double>& zeta,
549 Vector<double>& f)
550 {
551 // position vector on sphere
552 Vector<double> on_sphere(3);
553 Vector<double> temp_zeta(2);
554 temp_zeta[0] = -1.0;
555 temp_zeta[1] = zeta[1];
556 r_right_R(t, temp_zeta, on_sphere);
557
558 // position vector on center box
559 Vector<double> on_center(3);
560 on_center[0] = 0.5 * Radius;
561 on_center[1] = 0.0;
562 on_center[2] = 0.5 * Radius * 0.5 * (zeta[1] + 1.0);
563 // strait line across
564 for (unsigned i = 0; i < 3; i++)
565 {
566 f[i] =
567 on_center[i] + 0.5 * (zeta[0] + 1.0) * (on_sphere[i] - on_center[i]);
568 }
569 }
570
571 //=================================================================
572 /// Boundary of right box macro element
573 /// zeta \f$ \in [-1,1]^2 \f$
574 //=================================================================
575 void EighthSphereDomain::r_right_U(const unsigned& t,
576 const Vector<double>& zeta,
577 Vector<double>& f)
578 {
579 // position vector on sphere
580 Vector<double> on_sphere(3);
581 Vector<double> temp_zeta(2);
582 temp_zeta[0] = 1.0;
583 temp_zeta[1] = zeta[1];
584 r_right_R(t, temp_zeta, on_sphere);
585
586 // position vector on center box
587 Vector<double> on_center(3);
588 on_center[0] = 0.5 * Radius;
589 on_center[1] = 0.5 * Radius;
590 on_center[2] = 0.5 * Radius * 0.5 * (zeta[1] + 1.0);
591 // strait line across
592 for (unsigned i = 0; i < 3; i++)
593 {
594 f[i] =
595 on_center[i] + 0.5 * (zeta[0] + 1.0) * (on_sphere[i] - on_center[i]);
596 }
597 }
598
599 //=================================================================
600 /// Boundary of right box macro element
601 /// zeta \f$ \in [-1,1]^2 \f$
602 //=================================================================
603 void EighthSphereDomain::r_right_B(const unsigned& t,
604 const Vector<double>& zeta,
605 Vector<double>& f)
606 {
607 // position vector on sphere
608 Vector<double> on_sphere(3);
609 Vector<double> temp_zeta(2);
610 temp_zeta[0] = zeta[1];
611 temp_zeta[1] = -1.0;
612 r_right_R(t, temp_zeta, on_sphere);
613
614 // position vector on center box
615 Vector<double> on_center(3);
616 on_center[0] = 0.5 * Radius;
617 on_center[1] = 0.5 * Radius * 0.5 * (zeta[1] + 1.0);
618 on_center[2] = 0.0;
619 // strait line across
620 for (unsigned i = 0; i < 3; i++)
621 {
622 f[i] =
623 on_center[i] + 0.5 * (zeta[0] + 1.0) * (on_sphere[i] - on_center[i]);
624 }
625 }
626
627 //=================================================================
628 /// Boundary of right box macro element
629 /// zeta \f$ \in [-1,1]^2 \f$
630 //=================================================================
631 void EighthSphereDomain::r_right_F(const unsigned& t,
632 const Vector<double>& zeta,
633 Vector<double>& f)
634 {
635 // position vector on sphere
636 Vector<double> on_sphere(3);
637 Vector<double> temp_zeta(2);
638 temp_zeta[0] = zeta[1];
639 temp_zeta[1] = 1.0;
640 r_right_R(t, temp_zeta, on_sphere);
641
642 // position vector on center box
643 Vector<double> on_center(3);
644 on_center[0] = 0.5 * Radius;
645 on_center[1] = 0.5 * Radius * 0.5 * (zeta[1] + 1.0);
646 on_center[2] = 0.5 * Radius;
647 // strait line across
648 for (unsigned i = 0; i < 3; i++)
649 {
650 f[i] =
651 on_center[i] + 0.5 * (zeta[0] + 1.0) * (on_sphere[i] - on_center[i]);
652 }
653 }
654
655 //=================================================================
656 /// Boundary of top left box macro element
657 /// zeta \f$ \in [-1,1]^2 \f$
658 //=================================================================
659 void EighthSphereDomain::r_up_L(const unsigned& t,
660 const Vector<double>& zeta,
661 Vector<double>& f)
662 {
663 // position vector on sphere
664 Vector<double> on_sphere(3);
665 Vector<double> temp_zeta(2);
666 temp_zeta[0] = -1.0;
667 temp_zeta[1] = zeta[1];
668 r_up_U(t, temp_zeta, on_sphere);
669
670 // position vector on center box
671 Vector<double> on_center(3);
672 on_center[0] = 0.0;
673 on_center[1] = 0.5 * Radius;
674 on_center[2] = 0.5 * Radius * 0.5 * (zeta[1] + 1.0);
675 // strait line across
676 for (unsigned i = 0; i < 3; i++)
677 {
678 f[i] =
679 on_center[i] + 0.5 * (zeta[0] + 1.0) * (on_sphere[i] - on_center[i]);
680 }
681 }
682
683 //=================================================================
684 /// Boundary of top left box macro element
685 /// zeta \f$ \in [-1,1]^2 \f$
686 //=================================================================
687 void EighthSphereDomain::r_up_R(const unsigned& t,
688 const Vector<double>& zeta,
689 Vector<double>& f)
690 {
691 r_right_U(t, zeta, f);
692 }
693
694 //=================================================================
695 /// Boundary of top left box macro element
696 /// zeta \f$ \in [-1,1]^2 \f$
697 //=================================================================
698 void EighthSphereDomain::r_up_D(const unsigned& t,
699 const Vector<double>& zeta,
700 Vector<double>& f)
701 {
702 r_centr_U(t, zeta, f);
703 }
704
705 //=================================================================
706 /// Boundary of top left box macro element
707 /// zeta \f$ \in [-1,1]^2 \f$
708 //=================================================================
709 void EighthSphereDomain::r_up_U(const unsigned& t,
710 const Vector<double>& zeta,
711 Vector<double>& f)
712 {
713 double k0 = 0.5 * (1.0 + zeta[0]);
714 double k1 = 0.5 * (1.0 + zeta[1]);
715 Vector<double> p(3);
716 Vector<double> point1(3);
717 Vector<double> point2(3);
718
719 point1[0] = 0.5 * Radius * k0;
720 point1[1] = Radius - 0.5 * Radius * k0;
721 point1[2] = 0;
722 point2[0] = k0 * Radius / 3.0;
723 point2[1] = 0.5 * Radius - k0 * Radius / 6.0;
724 point2[2] = 0.5 * Radius - k0 * Radius / 6.0;
725
726 for (unsigned i = 0; i < 3; i++)
727 {
728 p[i] = point1[i] + k1 * (point2[i] - point1[i]);
729 }
730 double alpha = Radius / std::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
731
732 for (unsigned i = 0; i < 3; i++)
733 {
734 f[i] = alpha * p[i];
735 }
736 }
737
738 //=================================================================
739 /// Boundary of top left box macro element
740 /// zeta \f$ \in [-1,1]^2 \f$
741 //=================================================================
742 void EighthSphereDomain::r_up_B(const unsigned& t,
743 const Vector<double>& zeta,
744 Vector<double>& f)
745 {
746 // position vector on sphere
747 Vector<double> on_sphere(3);
748 Vector<double> temp_zeta(2);
749 temp_zeta[0] = zeta[0];
750 temp_zeta[1] = -1.0;
751 r_up_U(t, temp_zeta, on_sphere);
752
753 // position vector on center box
754 Vector<double> on_center(3);
755 on_center[0] = 0.5 * Radius * 0.5 * (zeta[0] + 1.0);
756 on_center[1] = 0.5 * Radius;
757 on_center[2] = 0.0;
758 // strait line across
759 for (unsigned i = 0; i < 3; i++)
760 {
761 f[i] =
762 on_center[i] + 0.5 * (zeta[1] + 1.0) * (on_sphere[i] - on_center[i]);
763 }
764 }
765
766 //=================================================================
767 /// Boundary of top left box macro element
768 /// zeta \f$ \in [-1,1]^2 \f$
769 //=================================================================
770 void EighthSphereDomain::r_up_F(const unsigned& t,
771 const Vector<double>& zeta,
772 Vector<double>& f)
773 {
774 // position vector on sphere
775 Vector<double> on_sphere(3);
776 Vector<double> temp_zeta(2);
777 temp_zeta[0] = zeta[0];
778 temp_zeta[1] = 1.0;
779 r_up_U(t, temp_zeta, on_sphere);
780
781 // position vector on center box
782 Vector<double> on_center(3);
783 on_center[0] = 0.5 * Radius * 0.5 * (zeta[0] + 1.0);
784 on_center[1] = 0.5 * Radius;
785 on_center[2] = 0.5 * Radius;
786 // straight line across
787 for (unsigned i = 0; i < 3; i++)
788 {
789 f[i] =
790 on_center[i] + 0.5 * (zeta[1] + 1.0) * (on_sphere[i] - on_center[i]);
791 }
792 }
793
794 //=================================================================
795 /// Boundary of top left box macro element
796 /// zeta \f$ \in [-1,1]^2 \f$
797 //=================================================================
798 void EighthSphereDomain::r_front_L(const unsigned& t,
799 const Vector<double>& zeta,
800 Vector<double>& f)
801 {
802 // position vector on sphere
803 Vector<double> on_sphere(3);
804 Vector<double> temp_zeta(2);
805 temp_zeta[0] = -1.0;
806 temp_zeta[1] = zeta[0];
807 r_front_F(t, temp_zeta, on_sphere);
808
809 // position vector on center box
810 Vector<double> on_center(3);
811 on_center[0] = 0.0;
812 on_center[1] = 0.5 * Radius * 0.5 * (zeta[0] + 1.0);
813 on_center[2] = 0.5 * Radius;
814 // straight line across
815 for (unsigned i = 0; i < 3; i++)
816 {
817 f[i] =
818 on_center[i] + 0.5 * (zeta[1] + 1.0) * (on_sphere[i] - on_center[i]);
819 }
820 }
821
822 //=================================================================
823 /// Boundary of top left box macro element
824 /// zeta \f$ \in [-1,1]^2 \f$
825 //=================================================================
826 void EighthSphereDomain::r_front_R(const unsigned& t,
827 const Vector<double>& zeta,
828 Vector<double>& f)
829 {
830 Vector<double> zeta2(2);
831 zeta2[0] = zeta[1];
832 zeta2[1] = zeta[0];
833 r_right_F(t, zeta2, f);
834 }
835
836 //=================================================================
837 /// Boundary of top left box macro element
838 /// zeta \f$ \in [-1,1]^2 \f$
839 //=================================================================
840 void EighthSphereDomain::r_front_D(const unsigned& t,
841 const Vector<double>& zeta,
842 Vector<double>& f)
843 {
844 // position vector on sphere
845 Vector<double> on_sphere(3);
846 Vector<double> temp_zeta(2);
847 temp_zeta[0] = zeta[0];
848 temp_zeta[1] = -1.0;
849 r_front_F(t, temp_zeta, on_sphere);
850
851 // position vector on center box
852 Vector<double> on_center(3);
853 on_center[0] = 0.5 * Radius * 0.5 * (zeta[0] + 1.0);
854 on_center[1] = 0.0;
855 on_center[2] = 0.5 * Radius;
856 // straight line across
857 for (unsigned i = 0; i < 3; i++)
858 {
859 f[i] =
860 on_center[i] + 0.5 * (zeta[1] + 1.0) * (on_sphere[i] - on_center[i]);
861 }
862 }
863
864 //=================================================================
865 /// Boundary of top left box macro element
866 /// zeta \f$ \in [-1,1]^2 \f$
867 //=================================================================
868 void EighthSphereDomain::r_front_U(const unsigned& t,
869 const Vector<double>& zeta,
870 Vector<double>& f)
871 {
872 r_up_F(t, zeta, f);
873 }
874
875 //=================================================================
876 /// Boundary of top left box macro element
877 /// zeta \f$ \in [-1,1]^2 \f$
878 //=================================================================
879 void EighthSphereDomain::r_front_B(const unsigned& t,
880 const Vector<double>& zeta,
881 Vector<double>& f)
882 {
883 r_centr_F(t, zeta, f);
884 }
885
886 //=================================================================
887 /// Boundary of top left box macro element
888 /// zeta \f$ \in [-1,1]^2 \f$
889 //=================================================================
890 void EighthSphereDomain::r_front_F(const unsigned& t,
891 const Vector<double>& zeta,
892 Vector<double>& f)
893 {
894 double k0 = 0.5 * (1.0 + zeta[0]);
895 double k1 = 0.5 * (1.0 + zeta[1]);
896 Vector<double> p(3);
897 Vector<double> point1(3);
898 Vector<double> point2(3);
899
900 point1[0] = 0.5 * Radius * k0;
901 point1[1] = 0.0;
902 point1[2] = Radius - k0 * 0.5 * Radius;
903 point2[0] = k0 * Radius / 3.0;
904 point2[1] = 0.5 * Radius - k0 * Radius / 6.0;
905 point2[2] = 0.5 * Radius - k0 * Radius / 6.0;
906
907 for (unsigned i = 0; i < 3; i++)
908 {
909 p[i] = point1[i] + k1 * (point2[i] - point1[i]);
910 }
911 double alpha = Radius / std::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
912
913 for (unsigned i = 0; i < 3; i++)
914 {
915 f[i] = alpha * p[i];
916 }
917 }
918
919} // namespace oomph
920
921#endif
s
static char t char * s
Definition cfortran.h:568
i
cstr elem_len * i
Definition cfortran.h:603
t
char t
Definition cfortran.h:568
oomph::Domain
Base class for Domains with curvilinear and/or time-dependent boundaries. Domain boundaries are typic...
Definition domain.h:67
oomph::Domain::Macro_element_pt
Vector< MacroElement * > Macro_element_pt
Vector of pointers to macro elements.
Definition domain.h:301
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:575
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:490
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:798
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:438
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:709
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:477
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:840
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:514
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:826
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:890
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:547
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:770
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:631
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:879
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:451
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:742
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:603
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:503
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:687
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:698
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:425
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:868
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:659
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:464
oomph::OcTree::Direct_string
static Vector< std::string > Direct_string
Translate (enumerated) directions into strings.
Definition octree.h:329
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition oomph_definitions.h:229
oomph::OomphLibWarning
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Definition oomph_definitions.h:274
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition Tadvection_diffusion_reaction_elements.h:66
geom_objects.h
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition advection_diffusion_elements.cc:30