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-2022 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  //=================================================================
40  //=================================================================
41  class EighthSphereDomain : public Domain
42  {
43  public:
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 
61  EighthSphereDomain(const EighthSphereDomain&) = delete;
62 
64  void operator=(const EighthSphereDomain&) = delete;
65 
67  ~EighthSphereDomain() {}
68 
69 
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 
276  void r_centr_L(const unsigned& t,
277  const Vector<double>& zeta,
278  Vector<double>& f);
279 
282  void r_centr_R(const unsigned& t,
283  const Vector<double>& zeta,
284  Vector<double>& f);
285 
288  void r_centr_D(const unsigned& t,
289  const Vector<double>& zeta,
290  Vector<double>& f);
291 
294  void r_centr_U(const unsigned& t,
295  const Vector<double>& zeta,
296  Vector<double>& f);
297 
300  void r_centr_B(const unsigned& t,
301  const Vector<double>& zeta,
302  Vector<double>& f);
303 
306  void r_centr_F(const unsigned& t,
307  const Vector<double>& zeta,
308  Vector<double>& f);
309 
312  void r_right_L(const unsigned& t,
313  const Vector<double>& zeta,
314  Vector<double>& f);
315 
318  void r_right_R(const unsigned& t,
319  const Vector<double>& zeta,
320  Vector<double>& f);
321 
324  void r_right_D(const unsigned& t,
325  const Vector<double>& zeta,
326  Vector<double>& f);
327 
330  void r_right_U(const unsigned& t,
331  const Vector<double>& zeta,
332  Vector<double>& f);
333 
336  void r_right_B(const unsigned& t,
337  const Vector<double>& zeta,
338  Vector<double>& f);
339 
342  void r_right_F(const unsigned& t,
343  const Vector<double>& zeta,
344  Vector<double>& f);
345 
348  void r_up_L(const unsigned& t,
349  const Vector<double>& zeta,
350  Vector<double>& f);
351 
354  void r_up_R(const unsigned& t,
355  const Vector<double>& zeta,
356  Vector<double>& f);
357 
360  void r_up_D(const unsigned& t,
361  const Vector<double>& zeta,
362  Vector<double>& f);
363 
366  void r_up_U(const unsigned& t,
367  const Vector<double>& zeta,
368  Vector<double>& f);
369 
372  void r_up_B(const unsigned& t,
373  const Vector<double>& zeta,
374  Vector<double>& f);
375 
378  void r_up_F(const unsigned& t,
379  const Vector<double>& zeta,
380  Vector<double>& f);
381 
384  void r_front_L(const unsigned& t,
385  const Vector<double>& zeta,
386  Vector<double>& f);
387 
390  void r_front_R(const unsigned& t,
391  const Vector<double>& zeta,
392  Vector<double>& f);
393 
396  void r_front_D(const unsigned& t,
397  const Vector<double>& zeta,
398  Vector<double>& f);
399 
402  void r_front_U(const unsigned& t,
403  const Vector<double>& zeta,
404  Vector<double>& f);
405 
408  void r_front_B(const unsigned& t,
409  const Vector<double>& zeta,
410  Vector<double>& f);
411 
414  void r_front_F(const unsigned& t,
415  const Vector<double>& zeta,
416  Vector<double>& f);
417  };
418 
422 
423 
424  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
D
dominoes D
Definition: Domino.cpp:55
L
MatrixXd L
Definition: LLT_example.cpp:6
R
@ R
Definition: StatisticsVector.h:21
p
float * p
Definition: Tutorial_Map_using.cpp:9
oomph::Domain
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
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
Definition: eighth_sphere_domain.h:784
oomph::EighthSphereDomain::Radius
double Radius
Definition: eighth_sphere_domain.h:272
oomph::EighthSphereDomain::r_right_F
void r_right_F(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
Definition: eighth_sphere_domain.h:470
oomph::Matrix
Definition: matrices.h:74
oomph::OcTree::Direct_string
static Vector< std::string > Direct_string
Translate (enumerated) directions into strings.
Definition: octree.h:329
oomph::OomphLibError
Definition: oomph_definitions.h:222
oomph::OomphLibWarning
Definition: oomph_definitions.h:267
oomph::QMacroElement< 3 >
Definition: macro_element.h:373
f
static int f(const TensorMap< Tensor< int, 3 > > &tensor)
Definition: cxx11_tensor_map.cpp:237
s
RealScalar s
Definition: level1_cplx_impl.h:130
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:151
Eigen::zeta
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
RachelsAdvectionDiffusion::U
double U
Swimming speed.
Definition: two_d_variable_diff_adapt.cc:53
UniformPSDSelfTest.radius
radius
Definition: UniformPSDSelfTest.py:15
oomph::OcTreeNames::F
@ F
Definition: octree.h:74
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition: AnisotropicHookean.h:10
plotPSD.t
t
Definition: plotPSD.py:36
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86