quarter_circle_sector_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_QUARTER_CIRCLE_SECTOR_DOMAIN_HEADER
28 #define OOMPH_QUARTER_CIRCLE_SECTOR_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  /// Circular sector as domain. Domain is bounded by
39  /// curved boundary which is represented by a GeomObject. Domain is
40  /// parametrised by three macro elements.
41  //=================================================================
42  class QuarterCircleSectorDomain : public Domain
43  {
44  public:
45  /// Constructor: Pass boundary object and start and end coordinates
46  /// and fraction along boundary object where outer ring is divided.
47  QuarterCircleSectorDomain(GeomObject* boundary_geom_object_pt,
48  const double& xi_lo,
49  const double& fract_mid,
50  const double& xi_hi)
51  : Xi_lo(xi_lo),
52  Fract_mid(fract_mid),
53  Xi_hi(xi_hi),
54  Wall_pt(boundary_geom_object_pt),
56  {
57  // There are three macro elements
58  unsigned nmacro = 3;
59 
60  // Resize
61  Macro_element_pt.resize(nmacro);
62 
63  // Create macro elements
64  for (unsigned i = 0; i < nmacro; i++)
65  {
66  Macro_element_pt[i] = new QMacroElement<2>(this, i);
67  }
68  }
69 
70 
71  /// Broken copy constructor
73 
74  /// Broken assignment operator
75  void operator=(const QuarterCircleSectorDomain&) = delete;
76 
77  /// Destructor: empty; cleanup done in base class
79 
80  /// Typedef for function pointer for function that squashes
81  /// the outer two macro elements towards
82  /// the wall by mapping the input value of the "radial" macro element
83  /// coordinate to the return value
84  typedef double (*BLSquashFctPt)(const double& s);
85 
86 
87  /// Function pointer for function that squashes
88  /// the outer two macro elements towards
89  /// the wall by mapping the input value of the "radial" macro element
90  /// coordinate to the return value
92  {
93  return BL_squash_fct_pt;
94  }
95 
96 
97  /// Function that squashes the outer two macro elements towards
98  /// the wall by mapping the input value of the "radial" macro element
99  /// coordinate to the return value
100  double s_squashed(const double& s)
101  {
102  return BL_squash_fct_pt(s);
103  }
104 
105  /// Vector representation of the i_macro-th macro element
106  /// boundary i_direct (N/S/W/E) at time level t
107  /// (t=0: present; t>0: previous):
108  /// f(s). Note that the local coordinate \b s is a 1D
109  /// Vector rather than a scalar -- this is unavoidable because
110  /// this function implements the pure virtual function in the
111  /// Domain base class.
112  void macro_element_boundary(const unsigned& t,
113  const unsigned& i_macro,
114  const unsigned& i_direct,
115  const Vector<double>& s,
116  Vector<double>& f);
117 
118 
119  private:
120  /// Lower limit for the (1D) coordinates along the wall
121  double Xi_lo;
122 
123  /// Fraction along wall where outer ring is to be divided
124  double Fract_mid;
125 
126  /// Upper limit for the (1D) coordinates along the wall
127  double Xi_hi;
128 
129  /// Pointer to geometric object that represents the curved wall
130  GeomObject* Wall_pt;
131 
132  /// Function pointer for function that squashes
133  /// the outer two macro elements towards
134  /// the wall by mapping the input value of the "radial" macro element
135  /// coordinate to the return value
137 
138  /// Default for function that squashes
139  /// the outer two macro elements towards
140  /// the wall by mapping the input value of the "radial" macro element
141  /// coordinate to the return value: Identity.
142  static double default_BL_squash_fct(const double& s)
143  {
144  return s;
145  }
146 
147  /// Boundary of top left macro element zeta \f$ \in [-1,1] \f$
148  void r_top_left_N(const unsigned& t,
149  const Vector<double>& zeta,
150  Vector<double>& f);
151 
152  /// Boundary of top left macro element zeta \f$ \in [-1,1] \f$
153  void r_top_left_W(const unsigned& t,
154  const Vector<double>& zeta,
155  Vector<double>& f);
156 
157  /// Boundary of top left macro element zeta \f$ \in [-1,1] \f$
158  void r_top_left_S(const unsigned& t,
159  const Vector<double>& zeta,
160  Vector<double>& f);
161 
162  /// Boundary of top left macro element zeta \f$ \in [-1,1] \f$
163  void r_top_left_E(const unsigned& t,
164  const Vector<double>& zeta,
165  Vector<double>& f);
166 
167  /// Boundary of bottom right macro element zeta \f$ \in [-1,1] \f$
168  void r_bot_right_N(const unsigned& t,
169  const Vector<double>& zeta,
170  Vector<double>& f);
171 
172  /// Boundary of bottom right macro element zeta \f$ \in [-1,1] \f$
173  void r_bot_right_W(const unsigned& t,
174  const Vector<double>& zeta,
175  Vector<double>& f);
176 
177  /// Boundary of bottom right macro element zeta \f$ \in [-1,1] \f$
178  void r_bot_right_S(const unsigned& t,
179  const Vector<double>& zeta,
180  Vector<double>& f);
181 
182  /// Boundary of bottom right macro element zeta \f$ \in [-1,1] \f$
183  void r_bot_right_E(const unsigned& t,
184  const Vector<double>& zeta,
185  Vector<double>& f);
186 
187  /// Boundary of central box macro element zeta \f$ \in [-1,1] \f$
188  void r_centr_N(const unsigned& t,
189  const Vector<double>& zeta,
190  Vector<double>& f);
191 
192  /// Boundary of central box macro element zeta \f$ \in [-1,1] \f$
193  void r_centr_E(const unsigned& t,
194  const Vector<double>& zeta,
195  Vector<double>& f);
196 
197  /// Boundary of central box macro element zeta \f$ \in [-1,1] \f$
198  void r_centr_S(const unsigned& t,
199  const Vector<double>& zeta,
200  Vector<double>& f);
201 
202  /// Boundary of central box macro element zeta \f$ \in [-1,1] \f$
203  void r_centr_W(const unsigned& t,
204  const Vector<double>& zeta,
205  Vector<double>& f);
206  };
207 
208 
209  /// //////////////////////////////////////////////////////////////////////
210  /// //////////////////////////////////////////////////////////////////////
211  /// //////////////////////////////////////////////////////////////////////
212 
213 
214  //=================================================================
215  /// Vector representation of the imacro-th macro element
216  /// boundary idirect (N/S/W/E) at time level t (t=0: present; t>0: previous):
217  /// f(s)
218  //=================================================================
220  const unsigned& t,
221  const unsigned& imacro,
222  const unsigned& idirect,
223  const Vector<double>& s,
224  Vector<double>& f)
225  {
226 #ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
227  // Warn about time argument being moved to the front
228  OomphLibWarning(
229  "Order of function arguments has changed between versions 0.8 and 0.85",
230  "QuarterCircleSectorDomain::macro_element_boundary(...)",
231  OOMPH_EXCEPTION_LOCATION);
232 #endif
233 
234  // Which macro element?
235  // --------------------
236  switch (imacro)
237  {
238  using namespace QuadTreeNames;
239 
240  // Macro element 0: Central box
241  case 0:
242 
243  // Which direction?
244  if (idirect == N)
245  {
246  r_centr_N(t, s, f);
247  }
248  else if (idirect == S)
249  {
250  r_centr_S(t, s, f);
251  }
252  else if (idirect == W)
253  {
254  r_centr_W(t, s, f);
255  }
256  else if (idirect == E)
257  {
258  r_centr_E(t, s, f);
259  }
260  else
261  {
262  std::ostringstream error_stream;
263  error_stream << "idirect is " << idirect << " not one of N, S, E, W"
264  << std::endl;
265 
266  throw OomphLibError(error_stream.str(),
267  OOMPH_CURRENT_FUNCTION,
268  OOMPH_EXCEPTION_LOCATION);
269  }
270 
271  break;
272 
273  // Macro element 1: Bottom right
274  case 1:
275 
276  // Which direction?
277  if (idirect == N)
278  {
279  r_bot_right_N(t, s, f);
280  }
281  else if (idirect == S)
282  {
283  r_bot_right_S(t, s, f);
284  }
285  else if (idirect == W)
286  {
287  r_bot_right_W(t, s, f);
288  }
289  else if (idirect == E)
290  {
291  r_bot_right_E(t, s, f);
292  }
293  else
294  {
295  std::ostringstream error_stream;
296  error_stream << "idirect is " << idirect << " not one of N, S, E, W"
297  << std::endl;
298 
299  throw OomphLibError(error_stream.str(),
300  OOMPH_CURRENT_FUNCTION,
301  OOMPH_EXCEPTION_LOCATION);
302  }
303 
304  break;
305 
306  // Macro element 2:Top left
307  case 2:
308 
309  // Which direction?
310  if (idirect == N)
311  {
312  r_top_left_N(t, s, f);
313  }
314  else if (idirect == S)
315  {
316  r_top_left_S(t, s, f);
317  }
318  else if (idirect == W)
319  {
320  r_top_left_W(t, s, f);
321  }
322  else if (idirect == E)
323  {
324  r_top_left_E(t, s, f);
325  }
326  else
327  {
328  std::ostringstream error_stream;
329  error_stream << "idirect is " << idirect << " not one of N, S, E, W"
330  << std::endl;
331 
332  throw OomphLibError(error_stream.str(),
333  OOMPH_CURRENT_FUNCTION,
334  OOMPH_EXCEPTION_LOCATION);
335  }
336 
337  break;
338 
339  default:
340 
341  // Error
342  std::ostringstream error_stream;
343  error_stream << "Wrong imacro " << imacro << std::endl;
344 
345  throw OomphLibError(
346  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
347  }
348  }
349 
350 
351  //=================================================================
352  /// Northern edge of top left macro element \f$ s \in [-1,1] \f$
353  //=================================================================
355  const Vector<double>& s,
356  Vector<double>& f)
357  {
358  Vector<double> x(1);
359 
360  // Coordinate along wall
361  x[0] = Xi_hi + (Fract_mid * (Xi_hi - Xi_lo) - Xi_hi) * 0.5 * (s[0] + 1.0);
362 
363  Wall_pt->position(t, x, f);
364  }
365 
366 
367  //=================================================================
368  /// Western edge of top left macro element \f$s \in [-1,1] \f$
369  //=================================================================
371  const Vector<double>& s,
372  Vector<double>& f)
373  {
374  Vector<double> x(1);
375 
376  // Top left corner
377  Vector<double> r_top(2);
378  x[0] = Xi_hi;
379 
380  Wall_pt->position(t, x, r_top);
381 
382  f[0] = 0.0;
383  f[1] = 0.5 * r_top[1] * (1.0 + s_squashed(0.5 * (s[0] + 1.0)));
384  }
385 
386 
387  //=================================================================
388  /// Southern edge of top left macro element \f$ s \in [-1,1] \f$
389  //=================================================================
391  const Vector<double>& s,
392  Vector<double>& f)
393  {
394  Vector<double> x(1);
395 
396  // Top left corner
397  Vector<double> r_top(2);
398  x[0] = Xi_hi;
399 
400  Wall_pt->position(t, x, r_top);
401 
402 
403  // Bottom right corner
404  Vector<double> r_bot(2);
405  x[0] = 0.0;
406 
407  Wall_pt->position(t, x, r_bot);
408 
409  f[0] = 0.5 * r_bot[0] * 0.5 * (s[0] + 1.0);
410  f[1] = 0.5 * r_top[1];
411  }
412 
413 
414  //=================================================================
415  /// Eastern edge of top left macro element \f$ s \in [-1,1] \f$
416  //=================================================================
418  const Vector<double>& s,
419  Vector<double>& f)
420  {
421  Vector<double> x(1);
422 
423  // Top left corner
424  Vector<double> r_top(2);
425  x[0] = Xi_hi;
426 
427  Wall_pt->position(t, x, r_top);
428 
429  // Bottom right corner
430  Vector<double> r_bot(2);
431  x[0] = Xi_lo;
432 
433  Wall_pt->position(t, x, r_bot);
434 
435  // Halfway along wall
436  Vector<double> r_half(2);
437  x[0] = Xi_lo + Fract_mid * (Xi_hi - Xi_lo);
438 
439  Wall_pt->position(t, x, r_half);
440 
441  f[0] = 0.5 * (r_bot[0] + s_squashed(0.5 * (s[0] + 1.0)) *
442  (2.0 * r_half[0] - r_bot[0]));
443  f[1] = 0.5 * (r_top[1] + s_squashed(0.5 * (s[0] + 1.0)) *
444  (2.0 * r_half[1] - r_top[1]));
445  }
446 
447 
448  //=================================================================
449  /// Northern edge of bottom right macro element
450  //=================================================================
452  const Vector<double>& s,
453  Vector<double>& f)
454  {
455  r_top_left_E(t, s, f);
456  }
457 
458  //=================================================================
459  /// Western edge of bottom right macro element
460  //=================================================================
462  const Vector<double>& s,
463  Vector<double>& f)
464  {
465  Vector<double> x(1);
466 
467  // Top left corner
468  Vector<double> r_top(2);
469  x[0] = Xi_hi;
470 
471  Wall_pt->position(t, x, r_top);
472 
473  // Bottom right corner
474  Vector<double> r_bot(2);
475  x[0] = Xi_lo;
476 
477  Wall_pt->position(t, x, r_bot);
478 
479  f[0] = 0.5 * r_bot[0];
480  f[1] = 0.5 * r_top[1] * 0.5 * (s[0] + 1.0);
481  }
482 
483  //=================================================================
484  /// Southern edge of bottom right macro element
485  //=================================================================
487  const Vector<double>& s,
488  Vector<double>& f)
489  {
490  Vector<double> x(1);
491 
492  // Bottom right corner
493  Vector<double> r_bot(2);
494  x[0] = Xi_lo;
495  Wall_pt->position(t, x, r_bot);
496 
497 
498  f[0] = 0.5 * r_bot[0] * (1.0 + s_squashed(0.5 * (s[0] + 1.0)));
499  f[1] = 0.0;
500  }
501 
502  //=================================================================
503  /// Eastern edge of bottom right macro element
504  //=================================================================
506  const Vector<double>& s,
507  Vector<double>& f)
508  {
509  Vector<double> x(1);
510 
511  // Coordinate along wall
512  x[0] = Xi_lo + (Fract_mid * (Xi_hi - Xi_lo) - Xi_lo) * (s[0] + 1.0) * 0.5;
513 
514  Wall_pt->position(t, x, f);
515  }
516 
517 
518  //=================================================================
519  /// Northern edge of central box
520  //=================================================================
521  void QuarterCircleSectorDomain::r_centr_N(const unsigned& t,
522  const Vector<double>& s,
523  Vector<double>& f)
524  {
525  r_top_left_S(t, s, f);
526  }
527 
528 
529  //=================================================================
530  /// Eastern edge of central box
531  //=================================================================
532  void QuarterCircleSectorDomain::r_centr_E(const unsigned& t,
533  const Vector<double>& s,
534  Vector<double>& f)
535  {
536  r_bot_right_W(t, s, f);
537  }
538 
539 
540  //=================================================================
541  /// Southern edge of central box
542  //=================================================================
543  void QuarterCircleSectorDomain::r_centr_S(const unsigned& t,
544  const Vector<double>& s,
545  Vector<double>& f)
546  {
547  Vector<double> x(1);
548 
549  // Bottom right corner
550  Vector<double> r_bot(2);
551  x[0] = Xi_lo;
552  Wall_pt->position(t, x, r_bot);
553 
554  f[0] = 0.5 * r_bot[0] * 0.5 * (s[0] + 1.0);
555  f[1] = 0.0;
556  }
557 
558 
559  //=================================================================
560  /// Western edge of central box
561  //=================================================================
562  void QuarterCircleSectorDomain::r_centr_W(const unsigned& t,
563  const Vector<double>& s,
564  Vector<double>& f)
565  {
566  Vector<double> x(1);
567 
568  // Top left corner
569  Vector<double> r_top(2);
570  x[0] = Xi_hi;
571  Wall_pt->position(t, x, r_top);
572 
573  f[0] = 0.0;
574  f[1] = 0.5 * r_top[1] * 0.5 * (s[0] + 1.0);
575  }
576 
577 
578 } // namespace oomph
579 
580 #endif
Circular sector as domain. Domain is bounded by curved boundary which is represented by a GeomObject....
void r_top_left_S(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left macro element zeta .
QuarterCircleSectorDomain(const QuarterCircleSectorDomain &)=delete
Broken copy constructor.
void r_centr_N(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
void r_centr_W(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
void r_bot_right_S(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of bottom right macro element zeta .
GeomObject * Wall_pt
Pointer to geometric object that represents the curved wall.
double s_squashed(const double &s)
Function that squashes the outer two macro elements towards the wall by mapping the input value of th...
double(* BLSquashFctPt)(const double &s)
Typedef for function pointer for function that squashes the outer two macro elements towards the wall...
void operator=(const QuarterCircleSectorDomain &)=delete
Broken assignment operator.
void r_top_left_W(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left macro element zeta .
void r_bot_right_W(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of bottom right macro element zeta .
void r_top_left_E(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left macro element zeta .
void macro_element_boundary(const unsigned &t, const unsigned &i_macro, const unsigned &i_direct, const Vector< double > &s, Vector< double > &f)
Vector representation of the i_macro-th macro element boundary i_direct (N/S/W/E) at time level t (t=...
void r_centr_S(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
double Xi_hi
Upper limit for the (1D) coordinates along the wall.
void r_bot_right_E(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of bottom right macro element zeta .
QuarterCircleSectorDomain(GeomObject *boundary_geom_object_pt, const double &xi_lo, const double &fract_mid, const double &xi_hi)
Constructor: Pass boundary object and start and end coordinates and fraction along boundary object wh...
BLSquashFctPt & bl_squash_fct_pt()
Function pointer for function that squashes the outer two macro elements towards the wall by mapping ...
void r_bot_right_N(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of bottom right macro element zeta .
static double default_BL_squash_fct(const double &s)
Default for function that squashes the outer two macro elements towards the wall by mapping the input...
double Fract_mid
Fraction along wall where outer ring is to be divided.
void r_centr_E(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of central box macro element zeta .
BLSquashFctPt BL_squash_fct_pt
Function pointer for function that squashes the outer two macro elements towards the wall by mapping ...
~QuarterCircleSectorDomain()
Destructor: empty; cleanup done in base class.
void r_top_left_N(const unsigned &t, const Vector< double > &zeta, Vector< double > &f)
Boundary of top left macro element zeta .
double Xi_lo
Lower limit for the (1D) coordinates along the wall.
////////////////////////////////////////////////////////////////////// //////////////////////////////...