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