topologically_rectangular_domain.cc
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//====================================================================
27 
28 
29 namespace oomph
30 {
31  //=============================================================================
32  /// Constructor - domain boundaries are described with four boundary
33  /// function pointers describing the topology of the north, east, south, and
34  /// west boundaries
35  //=============================================================================
37  BoundaryFctPt north_pt,
38  BoundaryFctPt east_pt,
39  BoundaryFctPt south_pt,
40  BoundaryFctPt west_pt)
41  {
42  // domain comprises one macro element
43  Macro_element_pt.resize(1);
44 
45  // Create the macro element
46  Macro_element_pt[0] = new QMacroElement<2>(this, 0);
47 
48  // set boundary function pointers
49  North_boundary_fn_pt = north_pt;
50  East_boundary_fn_pt = east_pt;
51  South_boundary_fn_pt = south_pt;
52  West_boundary_fn_pt = west_pt;
53 
54  // by default derivative boundary function pointers are null
59 
60  // also by default second derivate boundary function pointers are null
65 
66  // paranoid self check to ensure that ends of boundaries meet
67 #ifdef PARANOID
68  // error message stream
69  std::ostringstream error_message;
70  bool error_flag = false;
71  // check NE corner
72  {
73  Vector<double> x_N(2);
74  (*North_boundary_fn_pt)(1.0, x_N);
75  Vector<double> x_E(2);
76  (*East_boundary_fn_pt)(1.0, x_E);
77  if (x_N[0] != x_E[0] || x_N[1] != x_E[1])
78  {
79  error_message << "North and East Boundaries do not meet at the "
80  << "North East Corner.\n"
81  << "North Boundary : x[0] = " << x_N[0] << "\n"
82  << " x[1] = " << x_N[1] << "\n"
83  << "East Boundary : x[0] = " << x_E[0] << "\n"
84  << " x[1] = " << x_E[1] << "\n\n";
85  error_flag = true;
86  }
87  }
88  // check SE corner
89  {
90  Vector<double> x_S(2);
91  (*South_boundary_fn_pt)(1.0, x_S);
92  Vector<double> x_E(2);
93  (*East_boundary_fn_pt)(-1.0, x_E);
94  if (x_S[0] != x_E[0] || x_S[1] != x_E[1])
95  {
96  error_message << "South and East Boundaries do not meet at the "
97  << "South East Corner.\n"
98  << "South Boundary : x[0] = " << x_S[0] << "\n"
99  << " x[1] = " << x_S[1] << "\n"
100  << "East Boundary : x[0] = " << x_E[0] << "\n"
101  << " x[1] = " << x_E[1] << "\n\n";
102  error_flag = true;
103  }
104  }
105  // check SW corner
106  {
107  Vector<double> x_S(2);
108  (*South_boundary_fn_pt)(-1.0, x_S);
109  Vector<double> x_W(2);
110  (*West_boundary_fn_pt)(-1.0, x_W);
111  if (x_S[0] != x_W[0] || x_S[1] != x_W[1])
112  {
113  error_message << "South and West Boundaries do not meet at the "
114  << "South West Corner.\n"
115  << "South Boundary : x[0] = " << x_S[0] << "\n"
116  << " x[1] = " << x_S[1] << "\n"
117  << "West Boundary : x[0] = " << x_W[0] << "\n"
118  << " x[1] = " << x_W[1] << "\n\n";
119  error_flag = true;
120  }
121  }
122  // check NW corner
123  {
124  Vector<double> x_N(2);
125  (*North_boundary_fn_pt)(-1.0, x_N);
126  Vector<double> x_W(2);
127  (*West_boundary_fn_pt)(1.0, x_W);
128  if (x_N[0] != x_W[0] || x_N[1] != x_W[1])
129  {
130  error_message << "North and West Boundaries do not meet at the "
131  << "North West Corner.\n"
132  << "North Boundary : x[0] = " << x_N[0] << "\n"
133  << " x[1] = " << x_N[1] << "\n"
134  << "West Boundary : x[0] = " << x_W[0] << "\n"
135  << " x[1] = " << x_W[1] << "\n\n";
136  error_flag = true;
137  }
138  }
139  if (error_flag)
140  {
141  throw OomphLibError(
142  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
143  }
144 #endif
145  }
146 
147  //=============================================================================
148  /// Constructor - takes length of domain in x and y direction as
149  /// arguements. Assumes domain is rectangular, and the south west (lower
150  /// left) corner is at 0,0.
151  //=============================================================================
153  const double& l_x, const double& l_y)
154  {
155  // domain comprises one macro element
156  Macro_element_pt.resize(1);
157 
158  // Create the macro element
159  Macro_element_pt[0] = new QMacroElement<2>(this, 0);
160 
161  // set the boundary function pointers to zero
166 
167  // resize x vectors
168  x_south_west.resize(2);
169  x_north_east.resize(2);
170 
171  // set south west corner to 0,0
172  x_south_west[0] = 0.0;
173  x_south_west[1] = 0.0;
174 
175  // set north east corner
176  x_north_east[0] = l_x;
177  x_north_east[1] = l_y;
178 
179  // by default derivative boundary function pointers are null
184 
185  // also by default second derivate boundary function pointers are null
190  }
191 
192  //=============================================================================
193  /// Constructor - takes the minimum and maximum coordinates of the
194  /// of an assumed rectanguler domain in the x and y direction
195  //=============================================================================
197  const double& x_min,
198  const double& x_max,
199  const double& y_min,
200  const double& y_max)
201  {
202  // domain comprises one macro element
203  Macro_element_pt.resize(1);
204 
205  // Create the macro element
206  Macro_element_pt[0] = new QMacroElement<2>(this, 0);
207 
208  // set the boundary function pointers to zero
213 
214  // resize x vectors
215  x_south_west.resize(2);
216  x_north_east.resize(2);
217 
218  // set vector values
219  x_south_west[0] = x_min;
220  x_south_west[1] = y_min;
221  x_north_east[0] = x_max;
222  x_north_east[1] = y_max;
223 
224  // by default derivative boundary function pointers are null
229 
230  // also by default second derivate boundary function pointers are null
235  }
236 
237  //=============================================================================
238  /// allows the boundary derivate function pointers to be set. To
239  /// compute the derivatives of the problem domain global coordinates (x_i) wrt
240  /// the macro element coordinates (m_i), dx_i/dm_t is required along the
241  /// domain boundaries (where dm_t is the macro element coordinate tangential
242  /// to the domain boundary). The derivatives dx_i/dm_t can either be
243  /// prescribed with function pointers, or if the function pointers are not
244  /// provided then dx_i/dm_t is computed with finite differencing.
245  /// Note - these functions are only required for domains contructed with
246  /// boundary function pointers
247  //=============================================================================
249  BoundaryFctPt d_north_pt,
250  BoundaryFctPt d_east_pt,
251  BoundaryFctPt d_south_pt,
252  BoundaryFctPt d_west_pt)
253  {
254  // set the boundary derivate function pointers
255  dNorth_boundary_fn_pt = d_north_pt;
256  dEast_boundary_fn_pt = d_east_pt;
257  dSouth_boundary_fn_pt = d_south_pt;
258  dWest_boundary_fn_pt = d_west_pt;
259  }
260 
261 
262  //=============================================================================
263  /// allows the boundary second derivate function pointers to be set.
264  /// To compute the second derivatives of the problem domain global
265  /// coordinates (x_i) wrt the macro element coordinates (m_i), d2x_i/dm_t^2
266  /// is required along the domain boundaries (where dm_t is the macro element
267  /// coordinate tangential to the domain boundary). The derivatives
268  /// d2x_i/dm_t^2 can either be prescribed with function pointers, or if the
269  /// function pointers are not provided then dx_i/dm_t is computed with finite
270  /// differencing.
271  /// Note - these functions are only required for domains contructed with
272  /// boundary function pointers
273  //=============================================================================
275  BoundaryFctPt d2_north_pt,
276  BoundaryFctPt d2_east_pt,
277  BoundaryFctPt d2_south_pt,
278  BoundaryFctPt d2_west_pt)
279  {
280  // set the boundary derivate function pointers
281  d2North_boundary_fn_pt = d2_north_pt;
282  d2East_boundary_fn_pt = d2_east_pt;
283  d2South_boundary_fn_pt = d2_south_pt;
284  d2West_boundary_fn_pt = d2_west_pt;
285  }
286 
287 
288  //=============================================================================
289  /// returns the global coordinate position (f) of macro element position s
290  /// on boundary i_direct (e.g. N/S/W/E in 2D) at time t (no time dependence)
291  //=============================================================================
293  const unsigned& t,
294  const unsigned& i_macro,
295  const unsigned& i_direct,
296  const Vector<double>& s,
297  Vector<double>& f)
298  {
299  // use quad tree edge names to label edge of domain
300  using namespace QuadTreeNames;
301 
302 #ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
303  // Warn about time argument being moved to the front
304  OomphLibWarning(
305  "Order of function arguments has changed between versions 0.8 and 0.85",
306  "TopologicallyRectangularDomain::macro_element_boundary(...)",
307  OOMPH_EXCEPTION_LOCATION);
308 #endif
309 
310  // north boundary
311  if (i_direct == N) r_N(s, f);
312  // east boundary
313  else if (i_direct == E)
314  r_E(s, f);
315  // south boundary
316  else if (i_direct == S)
317  r_S(s, f);
318  // west boundary
319  else if (i_direct == W)
320  r_W(s, f);
321  }
322 
323  //=============================================================================
324  /// returns the derivates of the global coordinate position (f) wrt to the
325  /// macro element coordinate at macro macro element position s on boundary
326  /// i_direct (e.g. N/S/W/E in 2D) at time t (no time dependence)
327  //=============================================================================
329  const unsigned& t,
330  const unsigned& i_macro,
331  const unsigned& i_direct,
332  const Vector<double>& s,
333  Vector<double>& f)
334  {
335  // use quad tree edge names to label edge of domain
336  using namespace QuadTreeNames;
337 
338 #ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
339  // Warn about time argument being moved to the front
340  OomphLibWarning(
341  "Order of function arguments has changed between versions 0.8 and 0.85",
342  "TopologicallyRectangularDomain::dmacro_element_boundary(...)",
343  OOMPH_EXCEPTION_LOCATION);
344 #endif
345 
346  // north boundary
347  if (i_direct == N) dr_N(s, f);
348  // east boundary
349  else if (i_direct == E)
350  dr_E(s, f);
351  // south boundary
352  else if (i_direct == S)
353  dr_S(s, f);
354  // west boundary
355  else if (i_direct == W)
356  dr_W(s, f);
357  }
358 
359  //=============================================================================
360  /// returns the second derivates of the global coordinate position (f) wrt to
361  /// the macro element coordinate at macro macro element position s on boundary
362  /// i_direct (e.g. N/S/W/E in 2D) at time t (no time dependence)
363  //=============================================================================
365  const unsigned& t,
366  const unsigned& i_macro,
367  const unsigned& i_direct,
368  const Vector<double>& s,
369  Vector<double>& f)
370  {
371  // use quad tree edge names to label edge of domain
372  using namespace QuadTreeNames;
373 
374 
375 #ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
376  // Warn about time argument being moved to the front
377  OomphLibWarning(
378  "Order of function arguments has changed between versions 0.8 and 0.85",
379  "TopologicallyRectangularDomain::d2macro_element_boundary(...)",
380  OOMPH_EXCEPTION_LOCATION);
381 #endif
382 
383  // north boundary
384  if (i_direct == N) d2r_N(s, f);
385  // east boundary
386  else if (i_direct == E)
387  d2r_E(s, f);
388  // south boundary
389  else if (i_direct == S)
390  d2r_S(s, f);
391  // west boundary
392  else if (i_direct == W)
393  d2r_W(s, f);
394  }
395 
396  //=============================================================================
397  /// takes the macro element coordinate position along the north
398  /// boundary and returns the global coordinate position along that boundary
399  //=============================================================================
400  void TopologicallyRectangularDomain::r_N(const Vector<double>& s,
401  Vector<double>& f)
402  {
403  if (North_boundary_fn_pt != 0)
404  {
405  (*North_boundary_fn_pt)(s[0], f);
406  }
407  else
408  {
409  f[0] =
410  x_south_west[0] + (s[0] + 1) / 2 * (x_north_east[0] - x_south_west[0]);
411  f[1] = x_north_east[1];
412  }
413  }
414 
415  //=============================================================================
416  /// takes the macro element coordinate position along the east
417  /// boundary and returns the global coordinate position along that boundary
418  //=============================================================================
419  void TopologicallyRectangularDomain::r_E(const Vector<double>& s,
420  Vector<double>& f)
421  {
422  if (East_boundary_fn_pt != 0)
423  {
424  (*East_boundary_fn_pt)(s[0], f);
425  }
426  else
427  {
428  f[0] = x_north_east[0];
429  f[1] =
430  x_south_west[1] + (s[0] + 1) / 2 * (x_north_east[1] - x_south_west[1]);
431  }
432  }
433 
434  //=============================================================================
435  /// takes the macro element coordinate position along the south
436  /// boundary and returns the global coordinate position along that boundary
437  //=============================================================================
438  void TopologicallyRectangularDomain::r_S(const Vector<double>& s,
439  Vector<double>& f)
440  {
441  if (South_boundary_fn_pt != 0)
442  {
443  (*South_boundary_fn_pt)(s[0], f);
444  }
445  else
446  {
447  f[0] =
448  x_south_west[0] + (s[0] + 1) / 2 * (x_north_east[0] - x_south_west[0]);
449  f[1] = x_south_west[1];
450  }
451  }
452 
453  //=============================================================================
454  /// takes the macro element coordinate position along the west
455  /// boundary and returns the global coordinate position along that boundary
456  /// access down boundary function pointer
457  //=============================================================================
458  void TopologicallyRectangularDomain::r_W(const Vector<double>& s,
459  Vector<double>& f)
460  {
461  if (West_boundary_fn_pt != 0)
462  {
463  (*West_boundary_fn_pt)(s[0], f);
464  }
465  else
466  {
467  f[0] = x_south_west[0];
468  f[1] =
469  x_south_west[1] + (s[0] + 1) / 2 * (x_north_east[1] - x_south_west[1]);
470  }
471  }
472 
473  //=============================================================================
474  /// takes the macro element coordinate position along the north
475  /// boundary and returns the derivates of the global coordinates with respect
476  /// to the boundary
477  //=============================================================================
478  void TopologicallyRectangularDomain::dr_N(const Vector<double>& s,
479  Vector<double>& dr)
480  {
481  // if N boundary fn provided
482  if (North_boundary_fn_pt != 0)
483  {
484  // if dN boundary fn provided
485  if (dNorth_boundary_fn_pt != 0)
486  {
487  (*dNorth_boundary_fn_pt)(s[0], dr);
488  }
489  // else compute by finite differencing
490  else
491  {
492  const double h = 10e-8;
493  Vector<double> x_N_left(2);
494  (*North_boundary_fn_pt)(s[0] - 0.5 * h, x_N_left);
495  Vector<double> x_N_right(2);
496  (*North_boundary_fn_pt)(s[0] + 0.5 * h, x_N_right);
497  dr[0] = (x_N_right[0] - x_N_left[0]) / h;
498  dr[1] = (x_N_right[1] - x_N_left[1]) / h;
499  }
500  }
501  // if N boundary fn not provided then mesh is rectangular
502  else
503  {
504  dr[0] = (x_north_east[0] - x_south_west[0]) / 2;
505  dr[1] = 0;
506  }
507  }
508 
509  //=============================================================================
510  /// takes the macro element coordinate position along the E
511  /// boundary and returns the derivates of the global coordinates with respect
512  /// to the boundary
513  //=============================================================================
514  void TopologicallyRectangularDomain::dr_E(const Vector<double>& s,
515  Vector<double>& dr)
516  {
517  // if E boundary fn provided
518  if (East_boundary_fn_pt != 0)
519  {
520  // if dE boundary fn provided
521  if (dEast_boundary_fn_pt != 0)
522  {
523  (*dEast_boundary_fn_pt)(s[0], dr);
524  }
525  // else compute by finite differencing
526  else
527  {
528  const double h = 10e-8;
529  Vector<double> x_E_down(2);
530  (*East_boundary_fn_pt)(s[0] - 0.5 * h, x_E_down);
531  Vector<double> x_E_up(2);
532  (*East_boundary_fn_pt)(s[0] + 0.5 * h, x_E_up);
533  dr[0] = (x_E_up[0] - x_E_down[0]) / h;
534  dr[1] = (x_E_up[1] - x_E_down[1]) / h;
535  }
536  }
537  // if E boundary fn not provided then mesh is rectangular
538  else
539  {
540  dr[0] = 0;
541  dr[1] = (x_north_east[1] - x_south_west[1]) / 2;
542  }
543  }
544 
545  //=============================================================================
546  /// takes the macro element coordinate position along the south
547  /// boundary and returns the derivates of the global coordinates with respect
548  /// to the boundary
549  //=============================================================================
550  void TopologicallyRectangularDomain::dr_S(const Vector<double>& s,
551  Vector<double>& dr)
552  {
553  // if S boundary fn provided
554  if (South_boundary_fn_pt != 0)
555  {
556  // if dS boundary fn provided
557  if (dSouth_boundary_fn_pt != 0)
558  {
559  (*dSouth_boundary_fn_pt)(s[0], dr);
560  }
561  // else compute by finite differencing
562  else
563  {
564  const double h = 10e-8;
565  Vector<double> x_N_left(2);
566  (*South_boundary_fn_pt)(s[0] - 0.5 * h, x_N_left);
567  Vector<double> x_N_right(2);
568  (*South_boundary_fn_pt)(s[0] + 0.5 * h, x_N_right);
569  dr[0] = (x_N_right[0] - x_N_left[0]) / h;
570  dr[1] = (x_N_right[1] - x_N_left[1]) / h;
571  }
572  }
573  // if S boundary fn not provided then mesh is rectangular
574  else
575  {
576  dr[0] = (x_north_east[0] - x_south_west[0]) / 2;
577  dr[1] = 0;
578  }
579  }
580 
581  //=============================================================================
582  /// takes the macro element coordinate position along the W
583  /// boundary and returns the derivates of the global coordinates with respect
584  /// to the boundary
585  //=============================================================================
586  void TopologicallyRectangularDomain::dr_W(const Vector<double>& s,
587  Vector<double>& dr)
588  {
589  // if W boundary fn provided
590  if (West_boundary_fn_pt != 0)
591  {
592  // if dW boundary fn provided
593  if (dWest_boundary_fn_pt != 0)
594  {
595  (*dWest_boundary_fn_pt)(s[0], dr);
596  }
597  // else compute by finite differencing
598  else
599  {
600  const double h = 10e-8;
601  Vector<double> x_W_down(2);
602  (*West_boundary_fn_pt)(s[0] - 0.5 * h, x_W_down);
603  Vector<double> x_W_up(2);
604  (*West_boundary_fn_pt)(s[0] + 0.5 * h, x_W_up);
605  dr[0] = (x_W_up[0] - x_W_down[0]) / h;
606  dr[1] = (x_W_up[1] - x_W_down[1]) / h;
607  }
608  }
609  // if E boundary fn not provided then mesh is rectangular
610  else
611  {
612  dr[0] = 0;
613  dr[1] = (x_north_east[1] - x_south_west[1]) / 2;
614  }
615  }
616 
617  //=============================================================================
618  /// takes the macro element coordinate position along the north
619  /// boundary and returns the second derivates of the global coordinates with
620  /// respect to the boundary
621  //=============================================================================
622  void TopologicallyRectangularDomain::d2r_N(const Vector<double>& s,
623  Vector<double>& d2r)
624  {
625  // if N boundary fn provided
626  if (North_boundary_fn_pt != 0)
627  {
628  // if d2N boundary fn provided
629  if (d2North_boundary_fn_pt != 0)
630  {
631  (*d2North_boundary_fn_pt)(s[0], d2r);
632  }
633  // else compute by finite differencing
634  else
635  {
636  // if dN boundary fn provided finite difference d2N from it
637  if (dNorth_boundary_fn_pt != 0)
638  {
639  const double h = 10e-8;
640  Vector<double> dx_N_left(2);
641  (*dNorth_boundary_fn_pt)(s[0] - 0.5 * h, dx_N_left);
642  Vector<double> dx_N_right(2);
643  (*dNorth_boundary_fn_pt)(s[0] + 0.5 * h, dx_N_right);
644  d2r[0] = (dx_N_right[0] - dx_N_left[0]) / h;
645  d2r[1] = (dx_N_right[1] - dx_N_left[1]) / h;
646  }
647  // else finite difference from N boundary fn
648  else
649  {
650  const double h = 10e-8;
651  Vector<double> N_left(2);
652  (*North_boundary_fn_pt)(s[0] - h, N_left);
653  Vector<double> N_right(2);
654  (*North_boundary_fn_pt)(s[0] + h, N_right);
655  Vector<double> N_centre(2);
656  (*North_boundary_fn_pt)(s[0], N_centre);
657  d2r[0] = (N_right[0] + N_left[0] - 2 * N_centre[0]) / (h * h);
658  d2r[1] = (N_right[1] + N_left[1] - 2 * N_centre[1]) / (h * h);
659  }
660  }
661  }
662  // if N boundary fn not provided then mesh is rectangular
663  else
664  {
665  d2r[0] = 0;
666  d2r[1] = 0;
667  }
668  }
669 
670  //=============================================================================
671  /// takes the macro element coordinate position along the east
672  /// boundary and returns the second derivates of the global coordinates with
673  /// respect to the boundary
674  //=============================================================================
675  void TopologicallyRectangularDomain::d2r_E(const Vector<double>& s,
676  Vector<double>& d2r)
677  {
678  // if E boundary fn provided
679  if (East_boundary_fn_pt != 0)
680  {
681  // if d2E boundary fn provided
682  if (d2East_boundary_fn_pt != 0)
683  {
684  (*d2East_boundary_fn_pt)(s[0], d2r);
685  }
686  // else compute by finite differencing
687  else
688  {
689  // if dE boundary fn provided finite difference d2E from it
690  if (dEast_boundary_fn_pt != 0)
691  {
692  const double h = 10e-8;
693  Vector<double> dx_E_lower(2);
694  (*dEast_boundary_fn_pt)(s[0] - 0.5 * h, dx_E_lower);
695  Vector<double> dx_E_upper(2);
696  (*dEast_boundary_fn_pt)(s[0] + 0.5 * h, dx_E_upper);
697  d2r[0] = (dx_E_upper[0] - dx_E_lower[0]) / h;
698  d2r[1] = (dx_E_upper[1] - dx_E_lower[1]) / h;
699  }
700  // else finite difference from E boundary fn
701  else
702  {
703  const double h = 10e-8;
704  Vector<double> E_left(2);
705  (*East_boundary_fn_pt)(s[0] - h, E_left);
706  Vector<double> E_right(2);
707  (*East_boundary_fn_pt)(s[0] + h, E_right);
708  Vector<double> E_centre(2);
709  (*East_boundary_fn_pt)(s[0], E_centre);
710  d2r[0] = (E_right[0] + E_left[0] - 2 * E_centre[0]) / (h * h);
711  d2r[1] = (E_right[1] + E_left[1] - 2 * E_centre[1]) / (h * h);
712  }
713  }
714  }
715  // if E boundary fn not provided then mesh is rectangular
716  else
717  {
718  d2r[0] = 0;
719  d2r[1] = 0;
720  }
721  }
722 
723  //=============================================================================
724  /// takes the macro element coordinate position along the south
725  /// boundary and returns the second derivates of the global coordinates with
726  /// respect to the boundary
727  //=============================================================================
728  void TopologicallyRectangularDomain::d2r_S(const Vector<double>& s,
729  Vector<double>& d2r)
730  {
731  // if S boundary fn provided
732  if (South_boundary_fn_pt != 0)
733  {
734  // if d2S boundary fn provided
735  if (d2South_boundary_fn_pt != 0)
736  {
737  (*d2South_boundary_fn_pt)(s[0], d2r);
738  }
739  // else compute by finite differencing
740  else
741  {
742  // if dS boundary fn provided finite difference d2S from it
743  if (dSouth_boundary_fn_pt != 0)
744  {
745  const double h = 10e-8;
746  Vector<double> dx_S_left(2);
747  (*dSouth_boundary_fn_pt)(s[0] - 0.5 * h, dx_S_left);
748  Vector<double> dx_S_right(2);
749  (*dSouth_boundary_fn_pt)(s[0] + 0.5 * h, dx_S_right);
750  d2r[0] = (dx_S_right[0] - dx_S_left[0]) / h;
751  d2r[1] = (dx_S_right[1] - dx_S_left[1]) / h;
752  }
753  // else finite difference from S boundary fn
754  else
755  {
756  const double h = 10e-8;
757  Vector<double> S_left(2);
758  (*South_boundary_fn_pt)(s[0] - h, S_left);
759  Vector<double> S_right(2);
760  (*South_boundary_fn_pt)(s[0] + h, S_right);
761  Vector<double> S_centre(2);
762  (*South_boundary_fn_pt)(s[0], S_centre);
763  d2r[0] = (S_right[0] + S_left[0] - 2 * S_centre[0]) / (h * h);
764  d2r[1] = (S_right[1] + S_left[1] - 2 * S_centre[1]) / (h * h);
765  }
766  }
767  }
768  // if S boundary fn not provided then mesh is rectangular
769  else
770  {
771  d2r[0] = 0;
772  d2r[1] = 0;
773  }
774  }
775 
776  //=============================================================================
777  /// takes the macro element coordinate position along the west
778  /// boundary and returns the second derivates of the global coordinates with
779  /// respect to the boundary
780  //=============================================================================
781  void TopologicallyRectangularDomain::d2r_W(const Vector<double>& s,
782  Vector<double>& d2r)
783  {
784  // if W boundary fn provided
785  if (West_boundary_fn_pt != 0)
786  {
787  // if d2W boundary fn provided
788  if (d2West_boundary_fn_pt != 0)
789  {
790  (*d2West_boundary_fn_pt)(s[0], d2r);
791  }
792  // else compute by finite differencing
793  else
794  {
795  // if dW boundary fn provided finite difference d2W from it
796  if (dWest_boundary_fn_pt != 0)
797  {
798  const double h = 10e-8;
799  Vector<double> dx_W_lower(2);
800  (*dWest_boundary_fn_pt)(s[0] - 0.5 * h, dx_W_lower);
801  Vector<double> dx_W_upper(2);
802  (*dWest_boundary_fn_pt)(s[0] + 0.5 * h, dx_W_upper);
803  d2r[0] = (dx_W_upper[0] - dx_W_lower[0]) / h;
804  d2r[1] = (dx_W_upper[1] - dx_W_lower[1]) / h;
805  }
806  // else finite difference from W boundary fn
807  else
808  {
809  const double h = 10e-8;
810  Vector<double> W_left(2);
811  (*West_boundary_fn_pt)(s[0] - h, W_left);
812  Vector<double> W_right(2);
813  (*West_boundary_fn_pt)(s[0] + h, W_right);
814  Vector<double> W_centre(2);
815  (*West_boundary_fn_pt)(s[0], W_centre);
816  d2r[0] = (W_right[0] + W_left[0] - 2 * W_centre[0]) / (h * h);
817  d2r[1] = (W_right[1] + W_left[1] - 2 * W_centre[1]) / (h * h);
818  }
819  }
820  }
821  // if W boundary fn not provided then mesh is rectangular
822  else
823  {
824  d2r[0] = 0;
825  d2r[1] = 0;
826  }
827  }
828 } // namespace oomph
void r_E(const Vector< double > &s, Vector< double > &f)
takes the macro element coordinate position along the east boundary and returns the global coordinate...
void d2macro_element_boundary(const unsigned &t, const unsigned &i_macro, const unsigned &i_direct, const Vector< double > &s, Vector< double > &f)
returns the second derivates of the global coordinate position (f) wrt to the macro element coordinat...
BoundaryFctPt d2North_boundary_fn_pt
Function pointer to prescribe the second derivates of global coordinates wrt to the macro element coo...
void dmacro_element_boundary(const unsigned &t, const unsigned &i_macro, const unsigned &i_direct, const Vector< double > &s, Vector< double > &f)
returns the derivates of the global coordinate position (f) wrt to the macro element coordinate at ma...
void dr_E(const Vector< double > &s, Vector< double > &dr)
takes the macro element coordinate position along the E boundary and returns the derivates of the glo...
void dr_S(const Vector< double > &s, Vector< double > &dr)
takes the macro element coordinate position along the south boundary and returns the derivates of the...
Vector< double > x_south_west
coordinate position of south west corner of domain (only used if boundary functions are not used)
BoundaryFctPt d2East_boundary_fn_pt
Function pointer to prescribe the second derivates of global coordinates wrt to the macro element coo...
void d2r_W(const Vector< double > &s, Vector< double > &d2r)
takes the macro element coordinate position along the west boundary and returns the second derivates ...
BoundaryFctPt South_boundary_fn_pt
Function pointer to prescribe the north boundary of this topologically rectangular domain.
void macro_element_boundary(const unsigned &t, const unsigned &i_macro, const unsigned &i_direct, const Vector< double > &s, Vector< double > &f)
returns the global coordinate position (f) of macro element position s on boundary i_direct (e....
BoundaryFctPt dEast_boundary_fn_pt
Function pointer to prescribe the derivates of global coordinates wrt to the macro element coordinate...
BoundaryFctPt d2West_boundary_fn_pt
Function pointer to prescribe the second derivates of global coordinates wrt to the macro element coo...
void set_boundary_second_derivative_functions(BoundaryFctPt d2_north_pt, BoundaryFctPt d2_east_pt, BoundaryFctPt d2_south_pt, BoundaryFctPt d2_west_pt)
allows the boundary second derivate function pointers to be set. To compute the second derivatives of...
void r_N(const Vector< double > &s, Vector< double > &f)
takes the macro element coordinate position along the north boundary and returns the global coordinat...
void d2r_N(const Vector< double > &s, Vector< double > &d2r)
takes the macro element coordinate position along the north boundary and returns the second derivates...
BoundaryFctPt d2South_boundary_fn_pt
Function pointer to prescribe the second derivates of global coordinates wrt to the macro element coo...
void r_S(const Vector< double > &s, Vector< double > &f)
takes the macro element coordinate position along the south boundary and returns the global coordinat...
BoundaryFctPt North_boundary_fn_pt
Function pointer to prescribe the north boundary of this topologically rectangular domain.
void set_boundary_derivative_functions(BoundaryFctPt d_north_pt, BoundaryFctPt d_east_pt, BoundaryFctPt d_south_pt, BoundaryFctPt d_west_pt)
allows the boundary derivate function pointers to be set. To compute the derivatives of the problem d...
BoundaryFctPt dWest_boundary_fn_pt
Function pointer to prescribe the derivates of global coordinates wrt to the macro element coordinate...
void dr_W(const Vector< double > &s, Vector< double > &dr)
takes the macro element coordinate position along the W boundary and returns the derivates of the glo...
void r_W(const Vector< double > &s, Vector< double > &f)
takes the macro element coordinate position along the west boundary and returns the global coordinate...
void d2r_E(const Vector< double > &s, Vector< double > &d2r)
takes the macro element coordinate position along the east boundary and returns the second derivates ...
BoundaryFctPt West_boundary_fn_pt
Function pointer to prescribe the west boundary of this topologically rectangular domain.
BoundaryFctPt East_boundary_fn_pt
Function pointer to prescribe the east boundary of this topologically rectangular domain.
TopologicallyRectangularDomain(BoundaryFctPt north_pt, BoundaryFctPt east_pt, BoundaryFctPt south_pt, BoundaryFctPt west_pt)
Constructor - domain boundaries are described with four boundary function pointers describing the top...
void dr_N(const Vector< double > &s, Vector< double > &dr)
takes the macro element coordinate position along the north boundary and returns the derivates of the...
void d2r_S(const Vector< double > &s, Vector< double > &d2r)
takes the macro element coordinate position along the south boundary and returns the second derivates...
Vector< double > x_north_east
coordinate position of north east corner of domain (only used if boundary functions are not used)
BoundaryFctPt dSouth_boundary_fn_pt
Function pointer to prescribe the derivates of global coordinates wrt to the macro element coordinate...
BoundaryFctPt dNorth_boundary_fn_pt
Function pointer to prescribe the derivates of global coordinates wrt to the macro element coordinate...
////////////////////////////////////////////////////////////////////// //////////////////////////////...