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 
35 namespace 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
62 
63  /// Broken assignment operator
64  void operator=(const EighthSphereDomain&) = delete;
65 
66  /// Destructor: Empty; cleanup done in base class
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
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
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
Base class for Domains with curvilinear and/or time-dependent boundaries. Domain boundaries are typic...
Definition: domain.h:67
Vector< MacroElement * > Macro_element_pt
Vector of pointers to macro elements.
Definition: domain.h:301
Eighth sphere as domain. Domain is parametrised by four macro elements.
void r_right_U(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
void r_centr_F(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
void r_front_L(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
void r_centr_R(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
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 (...
void r_up_U(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
void r_centr_B(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
void r_front_D(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
void r_right_R(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
void r_front_R(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
EighthSphereDomain(const EighthSphereDomain &)=delete
Broken copy constructor.
void r_front_F(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
void r_right_D(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
EighthSphereDomain(const double &radius)
Constructor: Pass the radius of the sphere.
void operator=(const EighthSphereDomain &)=delete
Broken assignment operator.
void r_up_F(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
void r_right_F(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
void r_front_B(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
void r_centr_D(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
void r_up_B(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
~EighthSphereDomain()
Destructor: Empty; cleanup done in base class.
void r_right_B(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
void r_right_L(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of right box macro element zeta .
void r_up_R(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
void r_up_D(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
void r_centr_L(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
void r_front_U(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
void r_up_L(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left box macro element zeta .
void r_centr_U(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
static Vector< std::string > Direct_string
Translate (enumerated) directions into strings.
Definition: octree.h:329
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
QMacroElement specialised to 3 spatial dimensions.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...