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