channel_with_leaflet_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: 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//====================================================================
30 #ifndef OOMPH_CHANNEL_WITH_LEAFLET_DOMAIN_HEADER
31 #define OOMPH_CHANNEL_WITH_LEAFLET_DOMAIN_HEADER
32 
33 
34 // Generic includes
35 #include "../generic/geom_objects.h"
36 #include "../generic/domain.h"
37 #include "../generic/macro_element.h"
38 
39 namespace oomph
40 {
41 
42 //===========================================================
43 /// Rectangular domain with a leaflet blocking the lower
44 /// half
45 //===========================================================
47 {
48 
49 public:
50 
51  ///\short Constructor: Pass pointer to GeomObject that represents the leaflet,
52  /// the length of the domain to left and right of the leaflet, the
53  /// height of the leaflet and the overall height of the channel,
54  /// the number of element columns to the left and right of the leaflet,
55  /// the number of rows of elements from the bottom of the channel to
56  /// the end of the leaflet, the number of rows of elements above the
57  /// end of the leaflet.
59  const double& lright, const double& hleaflet,
60  const double& htot,
61  const unsigned& nleft, const unsigned& nright,
62  const unsigned& ny1, const unsigned& ny2)
63  {
64  // Copy assignments into private storage
65  Lleft = lleft;
66  Lright=lright;
68  Htot=htot;
69  Nleft = nleft;
70  Nright = nright;
71  Ny1 = ny1;
72  Ny2 = ny2;
74 
75  // Store origin of leaflet for fast reference
76  Vector<double> zeta(1);
77  zeta[0]=0.0;
78  Vector<double> r(2);
79  Leaflet_pt->position(zeta,r);
80  X_0 = r[0];
81 
82  // Number of macro elements
83  unsigned nmacro = (Ny1+Ny2)*(Nleft+Nright);
84  Macro_element_pt.resize(nmacro);
85 
86  // Build the 2D macro elements
87  for (unsigned i=0;i<nmacro;i++)
88  {
90  }
91 
92  }
93 
94 
95  /// Destructor: Kill macro elements
97  {
98  unsigned nmacro = (Ny1+Ny2)*(Nleft+Nright);
99  for (unsigned i=0;i<nmacro;i++){delete Macro_element_pt[i];}
100  }
101 
102 
103  /// Total height of domain (width of channel)
104  double htot(){return Htot;}
105 
106  /// Height of leaflet
107  double hleaflet(){return Hleaflet;}
108 
109  /// Length of domain to the left of leaflet
110  double lleft(){return Lleft;}
111 
112  /// Length of domain to the right of leaflet
113  double lright(){return Lright;}
114 
115  /// Pointer to the wall
117 
118  /// Parametrisation of macro element boundaries
119  void macro_element_boundary(const unsigned& t,
120  const unsigned& imacro,
121  const unsigned& idirect,
122  const Vector<double>& zeta,
123  Vector<double>& r);
124 
125 protected :
126 
127  /// Helper function
128  void macro_bound_I_N(const unsigned& t, const Vector<double>& zeta,
129  Vector<double>& r, const unsigned& i,const unsigned& j);
130 
131  /// Helper function
132  void macro_bound_I_S(const unsigned& t,const Vector<double>& zeta,
133  Vector<double>& r,const unsigned& i,const unsigned& j);
134 
135  /// Helper function
136  void macro_bound_I_W(const unsigned& t,const Vector<double>& zeta,
137  Vector<double>& r,const unsigned& i,const unsigned& j);
138 
139  /// Helper function
140  void macro_bound_I_E(const unsigned& t,const Vector<double>& zeta,
141  Vector<double>& r,const unsigned& i,
142  const unsigned& j);
143 
144  /// Helper function
145  void macro_bound_II_N(const unsigned& t, const Vector<double>& zeta,
146  Vector<double>& r,
147  const unsigned& i,const unsigned& j);
148 
149  /// Helper function
150  void macro_bound_II_S(const unsigned& t,const Vector<double>& zeta,
151  Vector<double>& r,const unsigned& i,const unsigned& j);
152 
153  /// Helper function
154  void macro_bound_II_W(const unsigned& t,const Vector<double>& zeta,
155  Vector<double>& r,const unsigned& i,const unsigned& j);
156 
157  /// Helper function
158  void macro_bound_II_E(const unsigned& t,const Vector<double>& zeta,
159  Vector<double>& r,const unsigned& i,
160  const unsigned& j);
161 
162  /// Helper function
163  void macro_bound_III_N(const unsigned& t, const Vector<double>& zeta,
164  Vector<double>& r,
165  const unsigned& i,const unsigned& j);
166 
167  /// Helper function
168  void macro_bound_III_S(const unsigned& t,const Vector<double>& zeta,
169  Vector<double>& r,const unsigned& i,const unsigned& j);
170 
171  /// Helper function
172  void macro_bound_III_W(const unsigned& t,const Vector<double>& zeta,
173  Vector<double>& r,const unsigned& i,const unsigned& j);
174 
175  /// Helper function
176  void macro_bound_III_E(const unsigned& t,const Vector<double>& zeta,
177  Vector<double>& r,const unsigned& i,
178  const unsigned& j);
179 
180  /// Helper function
181  void macro_bound_IV_N(const unsigned& t, const Vector<double>& zeta,
182  Vector<double>& r,
183  const unsigned& i,const unsigned& j);
184 
185  /// Helper function
186  void macro_bound_IV_S(const unsigned& t,const Vector<double>& zeta,
187  Vector<double>& r,const unsigned& i,const unsigned& j);
188 
189  /// Helper function
190  void macro_bound_IV_W(const unsigned& t,const Vector<double>& zeta,
191  Vector<double>& r,const unsigned& i,const unsigned& j);
192 
193  /// Helper function
194  void macro_bound_IV_E(const unsigned& t,const Vector<double>& zeta,
195  Vector<double>& r,const unsigned& i,
196  const unsigned& j);
197  /// Helper function
198  void slanted_bound_up(const unsigned& t,
199  const Vector<double>& zeta,
200  Vector<double>& r);
201 
202 
203 
204  /// Length of the domain to the right of the leaflet
205  double Lright;
206 
207  /// Length of the domain to the left of the leaflet
208  double Lleft;
209 
210  /// Lagrangian coordinate at end of leaflet
211  double Hleaflet;
212 
213  /// Total width of the channel
214  double Htot;
215 
216  ///Number of macro element columnns to the right of the leaflet
217  unsigned Nright;
218 
219  ///Number of macro element columns to the left of the leaflet
220  unsigned Nleft;
221 
222  ///Number of macro element rows up to the end of the leaflet
223  unsigned Ny1;
224 
225  ///Number of macro element rows above the leaflet
226  unsigned Ny2;
227 
228  /// \short Center of the domain : origin of the leaflet, extracted
229  /// from GeomObject and stored for fast access.
230  double X_0;
231 
232  /// Pointer to leaflet
234 
235 };
236 
237 
238 
239 //===================================================================
240 /// Parametrisation of macro element boundaries
241 //===================================================================
243  const unsigned& t,
244  const unsigned& imacro,
245  const unsigned& idirect,
246  const Vector<double>& zeta,
247  Vector<double>& r)
248 {
249 
250 #ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
251  // Warn about time argument being moved to the front
253  "Order of function arguments has changed between versions 0.8 and 0.85",
254  "ChannelWithLeafletDomain::macro_element_boundary(...)",
255  OOMPH_EXCEPTION_LOCATION);
256 #endif
257 
258  using namespace QuadTreeNames;
259 
260  //Number of the line and of the colum in the whole domain
261  //Beware : iline starts on zero
262  unsigned iline = int(floor(double(imacro)/double(Nleft+Nright)));
263  unsigned icol = imacro%(Nright+Nleft);
264 
265  //Number of the line and of the colum in the part considered
266  unsigned i,j;
267 
268  //Left low part of the domain : Part I
269  //------------------------------------
270  if( (iline<Ny1) && (icol<Nleft) )
271  {
272  i = iline;
273  j=icol;
274 
275  switch(idirect)
276  {
277  case N:
278  macro_bound_I_N(t,zeta,r,i,j);
279  break;
280  case S:
281  macro_bound_I_S(t,zeta,r,i,j);
282  break;
283  case W:
284  macro_bound_I_W(t,zeta,r,i,j);
285  break;
286  case E :
287  macro_bound_I_E(t,zeta,r,i,j);
288  break;
289  }
290  }
291  //Right low part of the domain : Part II
292  //--------------------------------------
293  else if( (iline<Ny1) && (icol>=Nleft) )
294  {
295  i = iline;
296  j=icol-Nleft;
297 
298  switch(idirect)
299  {
300  case N:
301  macro_bound_II_N(t,zeta,r,i,j);
302  break;
303  case S:
304  macro_bound_II_S(t,zeta,r,i,j);
305  break;
306  case W:
307  macro_bound_II_W(t,zeta,r,i,j);
308  break;
309  case E :
310  macro_bound_II_E(t,zeta,r,i,j);
311  break;
312  }
313  }
314  //Left upper part of the domain : Part III
315  //----------------------------------------
316  else if( (iline>=Ny1) && (icol<Nleft) )
317  {
318  i = iline-Ny1;
319  j=icol;
320 
321  switch(idirect)
322  {
323  case N:
324  macro_bound_III_N(t,zeta,r,i,j);
325  break;
326  case S:
327  macro_bound_III_S(t,zeta,r,i,j);
328  break;
329  case W:
330  macro_bound_III_W(t,zeta,r,i,j);
331  break;
332  case E :
333  macro_bound_III_E(t,zeta,r,i,j);
334  break;
335  }
336  }
337  //Right upper part of the domain : Part IV
338  //-----------------------------------------
339  else if( (iline>=Ny1) && (icol>=Nleft) )
340  {
341  i = iline-Ny1;
342  j=icol-Nleft;
343 
344  switch(idirect)
345  {
346  case N:
347  macro_bound_IV_N(t,zeta,r,i,j);
348  break;
349  case S:
350  macro_bound_IV_S(t,zeta,r,i,j);
351  break;
352  case W:
353  macro_bound_IV_W(t,zeta,r,i,j);
354  break;
355  case E :
356  macro_bound_IV_E(t,zeta,r,i,j);
357  break;
358  }
359  }
360 }
361 
362 
363 
364 /////////////////////////////////////////////////////////////////////
365 /////////////////////////////////////////////////////////////////////
366 // Helper functions for region I (lower left region)
367 /////////////////////////////////////////////////////////////////////
368 /////////////////////////////////////////////////////////////////////
369 
370 
371 
372 //=====================================================================
373 /// Helper function for eastern boundary in lower left region
374 //=====================================================================
376  const Vector<double>& zeta,
377  Vector<double>& r,
378  const unsigned& i,
379  const unsigned& j)
380 {
381 
382  //Find x,y on the wall corresponding to the position of the macro element
383 
384  //xi_wall varies from xi0 to xi1 on the wall
385  double xi0, xi1;
386  xi0 = double(i)*Hleaflet/double(Ny1);
387  xi1 = double(i+1)* Hleaflet/double(Ny1);
388 
389  Vector<double> xi_wall(1);
390  xi_wall[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
391 
392  Vector<double> r_wall(2);
393  Leaflet_pt->position(t,xi_wall,r_wall);
394 
395  //Find x,y on a vertical line corresponding
396  //to the position of the macro element
397 
398  //the vertical line goes from y0 to y1
399  double y0,y1;
400  y0 = double(i)*Hleaflet/double(Ny1);
401  y1 = double(i+1)*Hleaflet/double(Ny1);
402 
403  Vector<double> r_vert(2);
404  r_vert[0] = -Lleft+X_0;
405  r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
406 
407 //Parameter with value 0 in -Lleft and value 1 on the wall.
408  double s = double(j+1)/double(Nleft);
409 
410 ///Final expression of r
411  r[0] = r_vert[0] + s*( r_wall[0]-r_vert[0]);
412  r[1] = r_vert[1] + s*( r_wall[1]-r_vert[1]);
413 }
414 
415 
416 //=====================================================================
417 /// Helper function for western boundary in lower left region
418 //=====================================================================
420  const Vector<double>& zeta,
421  Vector<double>& r,
422  const unsigned& i,
423  const unsigned& j)
424 {
425  //Find x,y on the wall corresponding to the position of the macro element
426 
427  //xi_wall varies from xi0 to xi1 on the wall
428  double xi0, xi1;
429  xi0 = double(i)*Hleaflet/double(Ny1);
430  xi1 = double(i+1)*Hleaflet/double(Ny1);
431 
432  Vector<double> xi_wall(1);
433  xi_wall[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
434 
435  Vector<double> r_wall(2);
436  Leaflet_pt->position(t,xi_wall,r_wall);
437 
438  // Find x,y on a vertical line corresponding
439  //to the position of the macro element
440 
441  //the vertical line goes from y0 to y1
442  double y0,y1;
443  y0 = double(i)*Hleaflet/double(Ny1);
444  y1 = double(i+1)*Hleaflet/double(Ny1);
445 
446  Vector<double> r_vert(2);
447  r_vert[0] = -Lleft+X_0;
448  r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
449 
450  //Parameter with value 0 in -Lleft and value 1 on the wall.
451  double s = double(j)/double(Nleft);
452 
453  ///Final expression of r
454  r[0] = r_vert[0] + s*( r_wall[0]-r_vert[0]);
455  r[1] = r_vert[1] + s*( r_wall[1]-r_vert[1]);
456 
457 }
458 
459 
460 
461 //=====================================================================
462 /// Helper function for northern boundary in lower left region
463 //=====================================================================
465  const Vector<double>& zeta,
466  Vector<double>& r,
467  const unsigned& i,
468  const unsigned& j)
469 {
470  //Find the coordinates of the two corners of the north boundary
471  Vector<double> xi(1);
472  Vector<double> r_left(2);
473  Vector<double> r_right(2);
474  xi[0]=1;
475  macro_bound_I_W(t,xi,r_left,i,j);
476  macro_bound_I_E(t,xi,r_right,i,j);
477 
478 //Connect those two points with a straight line
479  r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
480  r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
481 
482 }
483 
484 
485 //=====================================================================
486 /// Helper function for southern boundary in lower left region
487 //=====================================================================
489  const Vector<double>& zeta,
490  Vector<double>& r,
491  const unsigned& i,
492  const unsigned& j)
493 {
494 ///Find the coordinates of the two corners of the south boundary
495  Vector<double> xi(1);
496  Vector<double> r_left(2);
497  Vector<double> r_right(2);
498  xi[0]=-1.0;
499  macro_bound_I_W(t,xi,r_left,i,j);
500  macro_bound_I_E(t,xi,r_right,i,j);
501 
502 //Connect those two points with a straight line
503  r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
504  r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
505 }
506 
507 
508 
509 
510 
511 
512 /////////////////////////////////////////////////////////////////////
513 /////////////////////////////////////////////////////////////////////
514 // Helper functions for region II (lower right region)
515 /////////////////////////////////////////////////////////////////////
516 /////////////////////////////////////////////////////////////////////
517 
518 
519 
520 //=====================================================================
521 /// Helper function for eastern boundary in lower right region
522 //=====================================================================
524  const Vector<double>& zeta,
525  Vector<double>& r,
526  const unsigned& i,
527  const unsigned& j)
528 
529 {
530  //Find x,y on the wall corresponding to the position of the macro element
531 
532  //xi_wall varies from xi0 to xi1 on the wall
533  double xi0, xi1;
534  xi0 = double(i)*Hleaflet/double(Ny1);
535  xi1 = double(i+1)*Hleaflet/double(Ny1);
536 
537  Vector<double> xi_wall(1);
538  xi_wall[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
539 
540  Vector<double> r_wall(2);
541  Leaflet_pt->position(t,xi_wall,r_wall);
542 
543  // Find x,y on a vertical line corresponding
544  // to the position of the macro element
545 
546  //the vertical line goes from y0 to y1
547  double y0,y1;
548  y0 = double(i)*Hleaflet/double(Ny1);
549  y1 = double(i+1)*Hleaflet/double(Ny1);
550 
551  Vector<double> r_vert(2);
552  r_vert[0] = Lright+X_0;
553  r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
554 
555 //Parameter with value 0 in Lright and value 1 on the wall.
556  double s = double(Nright-j-1)/double(Nright); /***Change****/
557 
558 ///Final expression of r
559  r[0] = r_vert[0] + s*( r_wall[0]-r_vert[0]);
560  r[1] = r_vert[1] + s*( r_wall[1]-r_vert[1]);
561 }
562 
563 
564 
565 
566 //=====================================================================
567 /// Helper function for western boundary in lower right region
568 //=====================================================================
570  const Vector<double>& zeta,
571  Vector<double>& r,
572  const unsigned& i,
573  const unsigned& j)
574 {
575  //Abscissa of the origin of the boudary east
576 
577  //Find x,y on the wall corresponding to the position of the macro element
578 
579  //xi_wall varies from xi0 to xi1 on the wall
580  double xi0, xi1;
581  xi0 = double(i)*Hleaflet/double(Ny1);
582  xi1 = double(i+1)*Hleaflet/double(Ny1);
583 
584  Vector<double> xi_wall(1);
585  xi_wall[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
586 
587  Vector<double> r_wall(2);
588  Leaflet_pt->position(t,xi_wall,r_wall);
589 
590  // Find x,y on a vertical line corresponding
591  // to the position of the macro element
592 
593  //the vertical line goes from y0 to y1
594  double y0,y1;
595  y0 = double(i)*Hleaflet/double(Ny1);
596  y1 = double(i+1)*Hleaflet/double(Ny1);
597 
598  Vector<double> r_vert(2);
599  r_vert[0] = Lright+X_0;
600  r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
601 
602  //Parameter with value 0 in -Lleft and value 1 on the wall.
603  double s = double(Nright-j)/double(Nright); /***Change****/
604 
605  //Final expression of r
606  r[0] = r_vert[0] + s*( r_wall[0]-r_vert[0]);
607  r[1] = r_vert[1] + s*( r_wall[1]-r_vert[1]);
608 
609 }
610 
611 
612 
613 //=====================================================================
614 /// Helper function for northern boundary in lower right region
615 //=====================================================================
617  const Vector<double>& zeta,
618  Vector<double>& r,
619  const unsigned& i,
620  const unsigned& j)
621 {
622  //Find the coordinates of the two corners of the north boundary
623  Vector<double> xi(1);
624  Vector<double> r_left(2);
625  Vector<double> r_right(2);
626  xi[0]=1;
627  macro_bound_II_W(t,xi,r_left,i,j);
628  macro_bound_II_E(t,xi,r_right,i,j);
629 
630 //Connect those two points with a straight line
631  r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
632  r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
633 
634 }
635 
636 
637 
638 //=====================================================================
639 /// Helper function for southern boundary in lower right region
640 //=====================================================================
642  const Vector<double>& zeta,
643  Vector<double>& r,
644  const unsigned& i,
645  const unsigned& j)
646 {
647  //Find the coordinates of the two corners of the south boundary
648  Vector<double> xi(1);
649  Vector<double> r_left(2);
650  Vector<double> r_right(2);
651  xi[0]=-1.0;
652  macro_bound_II_W(t,xi,r_left,i,j);
653  macro_bound_II_E(t,xi,r_right,i,j);
654 
655 //Connect those two points with a straight line
656  r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
657  r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
658 }
659 
660 
661 
662 
663 
664 
665 
666 /////////////////////////////////////////////////////////////////////
667 /////////////////////////////////////////////////////////////////////
668 // Helper functions for region III (upper left region)
669 /////////////////////////////////////////////////////////////////////
670 /////////////////////////////////////////////////////////////////////
671 
672 
673 
674 //=====================================================================
675 /// Describe the line between the boundary north of the domain (at x=X_0)
676 /// and the top of the wall, when zeta goes from 0 to 1.
677 //=====================================================================
679  const Vector<double>& zeta,
680  Vector<double>& r)
681 {
682  // Coordinates of the point on the boundary beetween the upper
683  // and the lower part, in the same column, at the east.
684  Vector<double> xi(1);
685  xi[0]=Hleaflet;
686 
687  Vector<double> r_join(2);
688 
689  Leaflet_pt->position(t,xi,r_join);
690 
691  r[0] = r_join[0] + zeta[0]*(X_0-r_join[0]);
692  r[1] = r_join[1] + zeta[0]*(Htot-r_join[1]);
693 }
694 
695 
696 
697 //=====================================================================
698 /// Helper function for eastern boundary in upper left region
699 //=====================================================================
701  const Vector<double>& zeta,
702  Vector<double>& r,
703  const unsigned& i,
704  const unsigned& j)
705 {
706  // Find x,y on the slanted straight line (SSL) corresponding to
707  // the position of the macro element
708 
709  //xi_line varies from xi0 to xi1 on the SSL
710  double xi0, xi1;
711  xi0 = double(i)/double(Ny2);
712  xi1 = double(i+1)/double(Ny2);
713 
714  Vector<double> xi_line(1);
715  xi_line[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
716 
717  Vector<double> r_line(2);
718  slanted_bound_up(t,xi_line,r_line);
719 
720  // Find x,y on a vertical line corresponding
721  // to the position of the macro element
722 
723  //the vertical line goes from y0 to y1
724  double y0,y1;
725  y0 = double(i)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;
726  y1 = double(i+1)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;;
727 
728  Vector<double> r_vert(2);
729  r_vert[0] = -Lleft+X_0;
730  r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
731 
732 //Parameter with value 0 in Lright and value 1 on the wall.
733  double s = double(j+1)/double(Nleft); /***Change****/
734 
735 ///Final expression of r
736  r[0] = r_vert[0] + s*( r_line[0]-r_vert[0]);
737  r[1] = r_vert[1] + s*( r_line[1]-r_vert[1]);
738 }
739 
740 
741 
742 
743 //=====================================================================
744 /// Helper function for western boundary in upper left region
745 //=====================================================================
747  const Vector<double>& zeta,
748  Vector<double>& r,
749  const unsigned& i,
750  const unsigned& j)
751 {
752  // Find x,y on the slanted straight line (SSL) corresponding to
753  // the position of the macro element
754 
755  //xi_line varies from xi0 to xi1 on the SSL
756  double xi0, xi1;
757  xi0 = double(i)/double(Ny2);
758  xi1 = double(i+1)/double(Ny2);
759 
760  Vector<double> xi_line(1);
761  xi_line[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
762 
763  Vector<double> r_line(2);
764  slanted_bound_up(t,xi_line,r_line);
765 
766  // Find x,y on a vertical line corresponding
767  // to the position of the macro element
768 
769  // the vertical line goes from y0 to y1
770  double y0,y1;
771  y0 = double(i)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;
772  y1 = double(i+1)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;;
773 
774  Vector<double> r_vert(2);
775  r_vert[0] = -Lleft+X_0;
776  r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
777 
778 //Parameter with value 0 in Lright and value 1 on the wall.
779  double s = double(j)/double(Nleft); /***Change****/
780 
781  //Final expression of r
782  r[0] = r_vert[0] + s*( r_line[0]-r_vert[0]);
783  r[1] = r_vert[1] + s*( r_line[1]-r_vert[1]);
784 }
785 
786 
787 
788 
789 //=====================================================================
790 /// Helper function for northern boundary in upper left region
791 //=====================================================================
793  const Vector<double>& zeta,
794  Vector<double>& r,
795  const unsigned& i,
796  const unsigned& j)
797 {
798  //Find the coordinates of the two corners of the north boundary
799  Vector<double> xi(1);
800  Vector<double> r_left(2);
801  Vector<double> r_right(2);
802  xi[0]=1;
803  macro_bound_III_W(t,xi,r_left,i,j);
804  macro_bound_III_E(t,xi,r_right,i,j);
805 
806  //Connect those two points with a straight line
807  r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
808  r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
809 
810 }
811 
812 
813 //=====================================================================
814 /// Helper function for southern boundary in upper left region
815 //=====================================================================
817  const Vector<double>& zeta,
818  Vector<double>& r,
819  const unsigned& i,
820  const unsigned& j)
821 {
822  //Find the coordinates of the two corners of the south boundary
823  Vector<double> xi(1);
824  Vector<double> r_left(2);
825  Vector<double> r_right(2);
826  xi[0]=-1;
827  macro_bound_III_W(t,xi,r_left,i,j);
828  macro_bound_III_E(t,xi,r_right,i,j);
829 
830  //Connect those two points with a straight line
831  r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
832  r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
833 }
834 
835 
836 
837 
838 
839 
840 
841 /////////////////////////////////////////////////////////////////////
842 /////////////////////////////////////////////////////////////////////
843 // Helper functions for region IV (upper right region)
844 /////////////////////////////////////////////////////////////////////
845 /////////////////////////////////////////////////////////////////////
846 
847 
848 //=====================================================================
849 /// Helper function for eastern boundary in upper right region
850 //=====================================================================
852  const Vector<double>& zeta,
853  Vector<double>& r,
854  const unsigned& i,
855  const unsigned& j)
856 {
857  // Find x,y on the slanted straight line (SSL) corresponding to
858  // the position of the macro element
859 
860  //xi_line varies from xi0 to xi1 on the SSL
861  double xi0, xi1;
862  xi0 = double(i)/double(Ny2);
863  xi1 = double(i+1)/double(Ny2);
864 
865  Vector<double> xi_line(1);
866  xi_line[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
867 
868  Vector<double> r_line(2);
869  slanted_bound_up(t,xi_line,r_line);
870 
871  // Find x,y on a vertical line corresponding
872  // to the position of the macro element
873 
874  //the vertical line goes from y0 to y1
875  double y0,y1;
876  y0 = double(i)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;
877  y1 = double(i+1)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;;
878 
879  Vector<double> r_vert(2);
880  r_vert[0] = Lright+X_0;
881  r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
882 
883  //Parameter with value 0 in Lright and value 1 on the wall.
884  double s = double(Nright -j-1)/double(Nright); /***Change****/
885 
886  //Final expression of r
887  r[0] = r_vert[0] + s*( r_line[0]-r_vert[0]);
888  r[1] = r_vert[1] + s*( r_line[1]-r_vert[1]);
889 }
890 
891 
892 
893 //=====================================================================
894 /// Helper function for western boundary in upper right region
895 //=====================================================================
897  const Vector<double>& zeta,
898  Vector<double>& r,
899  const unsigned& i,
900  const unsigned& j)
901 {
902  // Find x,y on the slanted straight line (SSL) corresponding to
903  // the position of the macro element
904 
905  //xi_line varies from xi0 to xi1 on the SSL
906  double xi0, xi1;
907  xi0 = double(i)/double(Ny2);
908  xi1 = double(i+1)/double(Ny2);
909 
910  Vector<double> xi_line(1);
911  xi_line[0] = xi0 + (1.0+zeta[0])/2.0*(xi1-xi0);
912 
913  Vector<double> r_line(2);
914  slanted_bound_up(t,xi_line,r_line);
915 
916  // Find x,y on a vertical line corresponding
917  // to the position of the macro element
918 
919  // The vertical line goes from y0 to y1
920  double y0,y1;
921  y0 = double(i)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;
922  y1 = double(i+1)*(Htot-Hleaflet)/double(Ny2)+Hleaflet;;
923 
924  Vector<double> r_vert(2);
925  r_vert[0] = Lright+X_0;
926  r_vert[1] = y0 + (1.0+zeta[0])/2.0*(y1-y0);
927 
928  //Parameter with value 0 in Lright and value 1 on the wall.
929  double s = double(Nright-j)/double(Nright); /***Change****/
930 
931  //Final expression of r
932  r[0] = r_vert[0] + s*( r_line[0]-r_vert[0]);
933  r[1] = r_vert[1] + s*( r_line[1]-r_vert[1]);
934 }
935 
936 
937 
938 //=====================================================================
939 /// Helper function for northern boundary in upper right region
940 //=====================================================================
942  const Vector<double>& zeta,
943  Vector<double>& r,
944  const unsigned& i,
945  const unsigned& j)
946 {
947  // Find the coordinates of the two corners of the north boundary
948  Vector<double> xi(1);
949  Vector<double> r_left(2);
950  Vector<double> r_right(2);
951  xi[0]=1;
952  macro_bound_IV_W(t,xi,r_left,i,j);
953  macro_bound_IV_E(t,xi,r_right,i,j);
954 
955  //Connect those two points with a straight line
956  r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
957  r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
958 
959 }
960 
961 
962 //=====================================================================
963 /// Helper function for southern boundary in upper right region
964 //=====================================================================
966  const Vector<double>& zeta,
967  Vector<double>& r,
968  const unsigned& i,
969  const unsigned& j)
970 {
971  //Find the coordinates of the two corners of the south boundary
972  Vector<double> xi(1);
973  Vector<double> r_left(2);
974  Vector<double> r_right(2);
975  xi[0]=-1;
976  macro_bound_IV_W(t,xi,r_left,i,j);
977  macro_bound_IV_E(t,xi,r_right,i,j);
978 
979  //Connect those two points with a straight line
980  r[0] = r_left[0] + (1.0+zeta[0])/2.0*(r_right[0]-r_left[0]);
981  r[1] = r_left[1] + (1.0+zeta[0])/2.0*(r_right[1]-r_left[1]);
982 }
983 
984 
985 }
986 
987 
988 #endif
989 
void macro_bound_III_S(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
double Htot
Total width of the channel.
cstr elem_len * i
Definition: cfortran.h:607
double Lright
Length of the domain to the right of the leaflet.
void macro_bound_II_E(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
Vector< MacroElement * > Macro_element_pt
Vector of pointers to macro elements.
Definition: domain.h:237
void macro_bound_IV_W(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
void macro_bound_II_S(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
void macro_bound_IV_N(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
void macro_bound_I_E(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
void macro_element_boundary(const unsigned &t, const unsigned &imacro, const unsigned &idirect, const Vector< double > &zeta, Vector< double > &r)
Parametrisation of macro element boundaries.
char t
Definition: cfortran.h:572
double Hleaflet
Lagrangian coordinate at end of leaflet.
void macro_bound_III_W(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
void macro_bound_IV_E(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
unsigned Ny1
Number of macro element rows up to the end of the leaflet.
GeomObject * Leaflet_pt
Pointer to leaflet.
void slanted_bound_up(const unsigned &t, const Vector< double > &zeta, Vector< double > &r)
Helper function.
void macro_bound_IV_S(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
double Lleft
Length of the domain to the left of the leaflet.
ChannelWithLeafletDomain(GeomObject *leaflet_pt, const double &lleft, const double &lright, const double &hleaflet, const double &htot, const unsigned &nleft, const unsigned &nright, const unsigned &ny1, const unsigned &ny2)
Constructor: Pass pointer to GeomObject that represents the leaflet, the length of the domain to left...
static char t char * s
Definition: cfortran.h:572
double X_0
Center of the domain : origin of the leaflet, extracted from GeomObject and stored for fast access...
void macro_bound_II_N(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
unsigned Nright
Number of macro element columnns to the right of the leaflet.
GeomObject *& leaflet_pt()
Pointer to the wall.
unsigned Nleft
Number of macro element columns to the left of the leaflet.
~ChannelWithLeafletDomain()
Destructor: Kill macro elements.
void macro_bound_II_W(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
double htot()
Total height of domain (width of channel)
unsigned Ny2
Number of macro element rows above the leaflet.
void macro_bound_I_N(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
void macro_bound_I_S(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
void macro_bound_I_W(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
double lright()
Length of domain to the right of leaflet.
double lleft()
Length of domain to the left of leaflet.
Base class for Domains with curvilinear and/or time-dependent boundaries. Domain boundaries are typic...
Definition: domain.h:71
void macro_bound_III_N(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.
void macro_bound_III_E(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &i, const unsigned &j)
Helper function.