collapsible_channel_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 //Include guards
31 #ifndef OOMPH_COLLAPSIBLE_CHANNEL_DOMAIN_HEADER
32 #define OOMPH_COLLAPSIBLE_CHANNEL_DOMAIN_HEADER
33 
34 // Generic oomph-lib includes
35 #include "../generic/quadtree.h"
36 #include "../generic/domain.h"
37 #include "../generic/geom_objects.h"
38 
39 namespace oomph
40 {
41 
42 //==================================================================
43 /// Collapsible channel domain
44 //==================================================================
45 class CollapsibleChannelDomain : public Domain
46 {
47 
48  public :
49 
50  /// \short Constructor: Pass the number of (macro-)elements,
51  /// the domain lengths in the x- and y-direction
52  /// and the pointer to the geometric object that specifies
53  /// the shape of the "collapsible" segment.
54  CollapsibleChannelDomain(const unsigned& nup,
55  const unsigned& ncollapsible,
56  const unsigned& ndown,
57  const unsigned& ny,
58  const double& lup,
59  const double& lcollapsible,
60  const double& ldown,
61  const double& ly,
62  GeomObject* wall_pt) :
65  {
66  Nup=nup;
68  Ndown=ndown;
69  Ny=ny;
70  Lup=lup;
71  Lcollapsible=lcollapsible;
72  Ldown=ldown;
73  Ly=ly;
75 
76  //Total number of macro elements
77  unsigned nmacro=(Nup+Ncollapsible+Ndown)*Ny;
78 
79  //Build the macro elements
80  Macro_element_pt.resize(nmacro);
81  for (unsigned i=0;i<nmacro;i++)
82  {
83  Macro_element_pt[i]= new QMacroElement<2>(this,i);
84  }
85  }
86 
87 
88  /// Destructor
90  {
91  // Kill macro elements
92  unsigned nelem=(Nup+Ncollapsible+Ndown)*Ny;
93  for (unsigned i=0;i<nelem;i++) delete Macro_element_pt[i];
94  }
95 
96 
97  /// Number of vertical columns of macro elements the upstream section
98  unsigned nup()
99  {
100  return Nup;
101  }
102 
103  /// Number of vertical clumns of macro elements in the "collapsible" segment
104  unsigned ncollapsible()
105  {
106  return Ncollapsible;
107  }
108 
109  /// Number of vertical columns of macro elements in the downstream section
110  unsigned ndown()
111  {
112  return Ndown;
113  }
114 
115  /// Number of macro-elements across the channel
116  unsigned ny()
117  {
118  return Ny;
119  }
120 
121  ///Length of upstream section
122  double l_up()
123  {
124  return Lup;
125  }
126 
127  ///\short Length of collapsible segment
128  double l_collapsible()
129  {
130  return Lcollapsible;
131  }
132 
133  ///Length of downstream section
134  double l_down()
135  {
136  return Ldown;
137  }
138 
139  /// Width of channel
140  double l_y()
141  {
142  return Ly;
143  }
144 
145  /// \short Access to pointer to the geometric object that parametrises
146  /// the collapsible wall
147  GeomObject*& wall_pt(){return Wall_pt;}
148 
149 
150  /// \short Access to pointer to the geometric object that parametrises
151  /// the collapsible wall (const version)
152  GeomObject* wall_pt() const {return Wall_pt;}
153 
154  /// \short Typedef for function pointer for function that squashes
155  /// the macro elements near the wall to help resolution of any
156  /// wall boundary layers.
157  typedef double (*BLSquashFctPt)(const double& s);
158 
159 
160  /// \short Default for function that squashes
161  /// the macro elements near the walls. Identity.
162  static double default_BL_squash_fct(const double& s)
163  {
164  return s;
165  }
166 
167  /// \short Function pointer for function that squashes
168  /// the macro elements near wall. Default mapping (identity)
169  /// leaves the y-coordinate of the nodal points unchanged.
171  {
172  return BL_squash_fct_pt;
173  }
174 
175  /// \short Function that squashes the macro elements near the wall.
176  /// Input argument should vary between 0 and 1; function should return
177  /// stretched/squashed coordinate in the same range. Default implementation
178  /// is the identity; can be overloaded by specifying a different
179  /// function pointer with bl_squash_fct_pt().
180  double s_squash(const double& s)
181  {
182  return BL_squash_fct_pt(s);
183  }
184 
185  /// \short Typedef for function pointer for function that implements
186  /// axial spacing of macro elements
187  typedef double (*AxialSpacingFctPt)(const double& xi);
188 
189  /// \short Function pointer for function that implements
190  /// axial spacing of macro elements
192  {
193  return Axial_spacing_fct_pt;
194  }
195 
196  /// \short Function that implements
197  /// axial spacing of macro elements
198  double axial_spacing_fct(const double& xi)
199  {
200  return Axial_spacing_fct_pt(xi);
201  }
202 
203 
204 /// \short Vector representation of the imacro-th macro element
205 /// boundary idirect (N/S/W/E) at time level t
206 /// (t=0: present; t>0: previous): \f$ {\bf r}({\bf zeta}) \f$
207 /// Note that the local coordinate \b zeta is a 1D
208 /// Vector rather than a scalar -- this is unavoidable because
209 /// this function implements the pure virtual function in the
210 /// Domain base class.
211 void macro_element_boundary(const unsigned& t,
212  const unsigned& imacro,
213  const unsigned& idirect,
214  const Vector<double>& zeta,
215  Vector<double>& r);
216 
217 
218  private :
219 
220  /// \short Northern boundary of the macro element imacro in the
221  /// upstream (part=0) or downstream (part=1) sections
222  void r_N_straight(const Vector<double>& zeta,
223  Vector<double>& r,
224  const unsigned& imacro,
225  const unsigned& part);
226 
227  /// \short Western boundary of the macro element imacro in the
228  /// upstream (part=0) or downstream (part=1) sections
229  void r_W_straight(const Vector<double>& zeta,
230  Vector<double>& r,
231  const unsigned& imacro,
232  const unsigned& part);
233 
234  /// \short Southern boundary of the macro element imacro in the
235  /// upstream (part=0) or downstream (part=1) sections
236  void r_S_straight(const Vector<double>& zeta,
237  Vector<double>& r,
238  const unsigned& imacro,
239  const unsigned& part);
240 
241  /// \short Eastern boundary of the macro element imacro in the
242  /// upstream (part=0) or downstream (part=1) sections
243  void r_E_straight(const Vector<double>& zeta,
244  Vector<double>& r,
245  const unsigned& imacro,
246  const unsigned& part );
247 
248  /// \short Northern boundary of the macro element imacro in the collapsible
249  /// section
250  void r_N_collapsible(const unsigned& t,
251  const Vector<double>& zeta,
252  Vector<double>& r,
253  const unsigned& imacro);
254 
255  /// \short Western boundary of the macro element imacro in the collapsible
256  /// section
257  void r_W_collapsible(const unsigned& t,
258  const Vector<double>& zeta,
259  Vector<double>& r,
260  const unsigned& imacro);
261 
262  /// \short Southern boundary of the macro element imacro in the collapsible
263  /// section
264  void r_S_collapsible(const unsigned& t,
265  const Vector<double>& zeta,
266  Vector<double>& r,
267  const unsigned& imacro );
268 
269  /// \short Eastern boundary of the macro element imacro in the collapsible
270  /// section
271  void r_E_collapsible(const unsigned& t,
272  const Vector<double>& zeta,
273  Vector<double>& r,
274  const unsigned& imacro );
275 
276 
277  /// \short Function pointer for function that squashes
278  /// the macro elements near the walls
280 
281  /// \short Function pointer for function that implements
282  /// axial spacing of macro elements
284 
285  /// \short Default for function that implements
286  /// axial spacing of macro elements
287  static double default_axial_spacing_fct(const double& xi)
288  {
289  return xi;
290  }
291 
292 
293  ///Number of vertical element columns in upstream section
294  unsigned Nup;
295 
296  ///Number of vertical element columns in "collapsible" section
297  unsigned Ncollapsible;
298 
299  ///Number of vertical element columns in downstream section
300  unsigned Ndown;
301 
302  ///Number of macro elements across channel
303  unsigned Ny;
304 
305  ///x-length in the upstream part of the channel
306  double Lup;
307 
308  ///x-length in the "collapsible" part of the channel
309  double Lcollapsible;
310 
311  ///x-length in the downstream part of the channel
312  double Ldown;
313 
314  ///Width
315  double Ly;
316 
317  ///Pointer to the geometric object that parametrises the collapsible wall
318  GeomObject* Wall_pt;
319 };
320 
321 
322 /////////////////////////////////////////////////////////////////////////
323 /////////////////////////////////////////////////////////////////////////
324 /////////////////////////////////////////////////////////////////////////
325 
326 //===================================================================
327 /// Vector representation of the imacro-th macro element
328 /// boundary idirect (N/S/W/E) at time level t
329 /// (t=0: present; t>0: previous): \f$ {\bf r}({\bf zeta}) \f$
330 /// Note that the local coordinate \b zeta is a 1D
331 /// Vector rather than a scalar -- this is unavoidable because
332 /// this function implements the pure virtual function in the
333 /// Domain base class.
334 //=================================================================
336  const unsigned& imacro,
337  const unsigned& idirect,
338  const Vector<double>& zeta,
339  Vector<double>& r)
340 {
341  using namespace QuadTreeNames;
342 
343 #ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
344  // Warn about time argument being moved to the front
345  OomphLibWarning(
346  "Order of function arguments has changed between versions 0.8 and 0.85",
347  "CollapsibleChannelDomain::macro_element_boundary(...)",
348  OOMPH_EXCEPTION_LOCATION);
349 #endif
350 
351  // Total number of vertical columns of (macro-)elements
352  unsigned n_x=Nup+Ncollapsible+Ndown;
353 
354  //Which direction?
355  if (idirect==N)
356  {
357  // Upstream part of the channel
358  if ((imacro%n_x)<Nup)
359  {
360  r_N_straight(zeta, r, imacro, 0);
361  }
362  // Downstream part of channel
363  else if ((imacro%n_x)>=Nup+Ncollapsible)
364  {
365  r_N_straight(zeta, r, imacro, 1);
366  }
367  // Collapsible segment
368  else if ( ((imacro%n_x)<Nup+Ncollapsible) && ((imacro%n_x)>=Nup) )
369  {
370  r_N_collapsible(t,zeta, r, imacro);
371  }
372  else
373  {
374  std::ostringstream error_stream;
375  error_stream << "Never get here! imacro, idirect: " << imacro
376  << " " << idirect << std::endl;
377 
378  throw OomphLibError(
379  error_stream.str(),
380  OOMPH_CURRENT_FUNCTION,
381  OOMPH_EXCEPTION_LOCATION);
382  }
383  }
384  else if (idirect==S)
385  {
386  // Upstream part
387  if ((imacro%n_x)<Nup)
388  {
389  r_S_straight(zeta, r, imacro, 0);
390  }
391  // Downstream part
392  else if ((imacro%n_x)>=Nup+Ncollapsible)
393  {
394  r_S_straight(zeta, r, imacro, 1);
395  }
396  // "Collapsible" bit
397  else if ( ((imacro%n_x)<Nup+Ncollapsible) && ((imacro%n_x)>=Nup) )
398  {
399  r_S_collapsible(t, zeta, r, imacro);
400  }
401  else
402  {
403  std::ostringstream error_stream;
404  error_stream << "Never get here! imacro, idirect: " << imacro
405  << " " << idirect << std::endl;
406 
407  throw OomphLibError(
408  error_stream.str(),
409  OOMPH_CURRENT_FUNCTION,
410  OOMPH_EXCEPTION_LOCATION);
411  }
412  }
413  else if (idirect==E)
414  {
415  // Upstream bit
416  if ((imacro%n_x)<Nup)
417  {
418  r_E_straight(zeta, r, imacro, 0);
419  }
420  // Downstream bit
421  else if ((imacro%n_x)>=Nup+Ncollapsible)
422  {
423  r_E_straight(zeta, r, imacro, 1);
424  }
425  // "Collapsible" bit
426  else if ( ((imacro%n_x)<Nup+Ncollapsible) && ((imacro%n_x)>=Nup) )
427  {
428  r_E_collapsible(t, zeta, r, imacro);
429  }
430  else
431  {
432  std::ostringstream error_stream;
433  error_stream << "Never get here! imacro, idirect: " << imacro
434  << " " << idirect << std::endl;
435 
436  throw OomphLibError(
437  error_stream.str(),
438  OOMPH_CURRENT_FUNCTION,
439  OOMPH_EXCEPTION_LOCATION);
440  }
441  }
442 
443  else if (idirect==W)
444  {
445  // Upstream bit
446  if ((imacro%n_x)<Nup)
447  {
448  r_W_straight(zeta, r, imacro, 0);
449  }
450  // Downstream bit
451  else if ((imacro%n_x)>=Nup+Ncollapsible)
452  {
453  r_W_straight(zeta, r, imacro, 1);
454  }
455  // "Collapsible" bit
456  else if ( ((imacro%n_x)<Nup+Ncollapsible) && ((imacro%n_x)>=Nup) )
457  {
458  r_W_collapsible(t,zeta, r, imacro);
459  }
460  else
461  {
462  std::ostringstream error_stream;
463  error_stream << "Never get here! imacro, idirect: " << imacro
464  << " " << idirect << std::endl;
465 
466  throw OomphLibError(
467  error_stream.str(),
468  OOMPH_CURRENT_FUNCTION,
469  OOMPH_EXCEPTION_LOCATION);
470  }
471  }
472 
473 };
474 
475 
476 
477 
478 
479 //===========================================================================
480 /// \short Western edge of the macro element in the upstream (part=0)
481 /// or downstream (part=1) parts of the channel; \f$ \zeta \in [-1,1] \f$
482 //===========================================================================
483 void CollapsibleChannelDomain::r_W_straight(const Vector<double>& zeta,
484  Vector<double>& r,
485  const unsigned& imacro,
486  const unsigned& part)
487 {
488 
489  //Determines the "coordinates" of the macro-element
490  unsigned x=unsigned(imacro%(Nup+Ncollapsible+Ndown));
491  unsigned y=unsigned(double(imacro)/double(Nup+Ncollapsible+Ndown));
492 
493  // Where are we?
494  switch(part)
495  {
496  case 0: //in the upstream part of the channel
497 
498  //Parametrize the boundary
499  r[0]=axial_spacing_fct(double(x)*(Lup/double(Nup)));
500  r[1]=(double(y)+(0.5*(1.0+zeta[0])))*(Ly/double(Ny));
501 
502  // Map it via squash fct
503  r[1]=Ly*s_squash(r[1]/Ly);
504 
505  break;
506 
507  case 1 : //in the downstream part of the channel
508 
509  //Parametrizes the boundary
510  r[0]=axial_spacing_fct(
511  double(x-Nup-Ncollapsible)*(Ldown/double(Ndown))+Lup+Lcollapsible);
512  r[1]=(double(y)+(0.5*(1.0+zeta[0])))*(Ly/double(Ny));
513 
514  // Map it via squash fct
515  r[1]=Ly*s_squash(r[1]/Ly);
516 
517  break;
518 
519  default:
520 
521  std::ostringstream error_stream;
522  error_stream << "Never get here! part=" << part << std::endl;
523 
524  throw OomphLibError(error_stream.str(),
525  OOMPH_CURRENT_FUNCTION,
526  OOMPH_EXCEPTION_LOCATION);
527  }
528 
529 }
530 
531 
532 
533 //===========================================================================
534 /// \short Eastern edge of the macro element in the straight parts
535 /// of the channel; \f$ \zeta \in [-1,1] \f$
536 /// part=0 in the upstream part, part=1 in the downstream part
537 //===========================================================================
538 void CollapsibleChannelDomain::r_E_straight(const Vector<double>& zeta,
539  Vector<double>& r,
540  const unsigned& imacro,
541  const unsigned& part)
542 {
543 
544  //Determines the "coordinates" of the macro-element
545  unsigned x=unsigned(imacro%(Nup+Ncollapsible+Ndown));
546  unsigned y=unsigned(double(imacro)/double(Nup+Ncollapsible+Ndown));
547 
548  // Where are we?
549  switch(part)
550  {
551 
552  case 0: //in the upstream part of the channel
553 
554  //Parametrizes the boundary
555  r[0]=axial_spacing_fct((double(x)+1.0)*(Lup/double(Nup)));
556  r[1]=(double(y)+(0.5*(1.0+zeta[0])))*(Ly/double(Ny));
557 
558  // Map it via squash fct
559  r[1]=Ly*s_squash(r[1]/Ly);
560 
561  break;
562 
563  case 1: //in the downstream part of the channel
564 
565  //Parametrizes the boundary
566  r[0]=axial_spacing_fct((double(x-Nup-Ncollapsible)+1.0)*
567  (Ldown/double(Ndown))+Lup+Lcollapsible);
568  r[1]=(double(y)+(0.5*(1.0+zeta[0])))*(Ly/double(Ny));
569 
570  // Map it via squash fct
571  r[1]=Ly*s_squash(r[1]/Ly);
572 
573  break;
574 
575  default:
576 
577 
578  std::ostringstream error_stream;
579  error_stream << "Never get here! part=" << part << std::endl;
580 
581  throw OomphLibError(error_stream.str(),
582  OOMPH_CURRENT_FUNCTION,
583  OOMPH_EXCEPTION_LOCATION);
584  }
585 }
586 
587 //==========================================================================
588 /// \short Northern edge of the macro element in the straight parts of
589 /// the channel; \f$ \zeta \in [-1,1] \f$
590 /// part=0 in the left part, part=1 in the right part
591 //==========================================================================
592 void CollapsibleChannelDomain::r_N_straight(const Vector<double>& zeta,
593  Vector<double>& r,
594  const unsigned& imacro,
595  const unsigned& part)
596 {
597 
598  //Determines the "coordinates" of the macro-element
599  unsigned x=unsigned(imacro%(Nup+Ncollapsible+Ndown));
600  unsigned y=unsigned(double(imacro)/double(Nup+Ncollapsible+Ndown));
601 
602 
603  // Where are we?
604  switch(part)
605  {
606 
607  case 0: //in the upstream part of the channel
608 
609  //Parametrizes the boundary
610  r[0]=axial_spacing_fct((double(x)+(0.5*(1.0+zeta[0])))*(Lup/double(Nup)));
611  r[1]=(double(y)+1.0)*(Ly/double(Ny));
612 
613  // Map it via squash fct
614  r[1]=Ly*s_squash(r[1]/Ly);
615 
616  break;
617 
618  case 1: //in the downstream part of the channel
619 
620  //Parametrizes the boundary
621  r[0]=axial_spacing_fct((double(x-Nup-Ncollapsible)+
622  (0.5*(1.0+zeta[0])))*
623  (Ldown/double(Ndown))+Lup+Lcollapsible);
624  r[1]=(double(y)+1.0)*(Ly/double(Ny));
625 
626  // Map it via squash fct
627  r[1]=Ly*s_squash(r[1]/Ly);
628 
629  break;
630 
631  default:
632 
633 
634  std::ostringstream error_stream;
635  error_stream << "Never get here! part=" << part << std::endl;
636 
637  throw OomphLibError(error_stream.str(),
638  OOMPH_CURRENT_FUNCTION,
639  OOMPH_EXCEPTION_LOCATION);
640  }
641 }
642 
643 
644 //=========================================================================
645 /// \short Southern edge of the macro element in the straight parts of
646 /// the channel; \f$ \zeta \in [-1,1] \f$
647 /// part=0 in the left part, part=1 in the right part
648 //=========================================================================
649 void CollapsibleChannelDomain::r_S_straight(const Vector<double>& zeta,
650  Vector<double>& r,
651  const unsigned& imacro,
652  const unsigned& part )
653 {
654 
655  //Determines the "coordinates" of the macro-element
656  unsigned x=unsigned(imacro%(Nup+Ncollapsible+Ndown));
657  unsigned y=unsigned(double(imacro)/double(Nup+Ncollapsible+Ndown));
658 
659  // Where are we?
660  switch(part)
661  {
662 
663  case 0: //in the upstream bit
664 
665  //Parametrizes the boundary
666  r[0]=axial_spacing_fct((double(x)+(0.5*(1+zeta[0])))*(Lup/double(Nup)));
667  r[1]=double(y)*(Ly/double(Ny));
668 
669  // Map it via squash fct
670  r[1]=Ly*s_squash(r[1]/Ly);
671 
672  break;
673 
674  case 1: //in the downstream bit
675 
676  //Parametrizes the boundary
677  r[0]=axial_spacing_fct((double(x-Nup-Ncollapsible)+
678  (0.5*(1+zeta[0])))*
679  (Ldown/double(Ndown))+Lup+Lcollapsible);
680  r[1]=double(y)*(Ly/double(Ny));
681 
682  // Map it via squash fct
683  r[1]=Ly*s_squash(r[1]/Ly);
684 
685  break;
686 
687  default:
688 
689 
690  std::ostringstream error_stream;
691  error_stream << "Never get here! part=" << part << std::endl;
692 
693  throw OomphLibError(error_stream.str(),
694  OOMPH_CURRENT_FUNCTION,
695  OOMPH_EXCEPTION_LOCATION);
696 
697  }
698 }
699 
700 
701 
702 
703 
704 
705 //========================================================================
706 /// Western edge of the macro element in the collapsible part of the
707 /// channel; \f$ \zeta \in [-1,1] \f$.
708 //========================================================================
710  const Vector<double>& zeta,
711  Vector<double>& r,
712  const unsigned& imacro)
713 {
714 
715 
716  //Determines the "coordinates" of the macro-element
717  unsigned x=unsigned(imacro%(Nup+Ncollapsible+Ndown));
718  unsigned y=unsigned(double(imacro)/double(Nup+Ncollapsible+Ndown));
719 
720  // Vector of Lagrangian coordinates
721  Vector<double> xi(1);
722  xi[0]=double(x-Nup)*(Lcollapsible/double(Ncollapsible));
723 
724  // Position vector on upper wall:
725  Vector<double> r_wall(2);
726  Wall_pt->position(t,xi,r_wall);
727 
728  // Point will be located on straight line from bottom to top wall
729  double fract=(double(y)+(0.5*(1.0+zeta[0])))/double(Ny);
730 
731  // Map it via squash fct
732  fract=s_squash(fract);
733 
734  // x-cooordinate -- straight line from fixed position on the bottom
735  // wall to moving position on the top wall
736  r[0]=Lup+xi[0]+(r_wall[0]-(xi[0]+Lup))*fract;
737 
738  // y-coordinate
739  r[1]=r_wall[1]*fract;
740 
741 }
742 
743 
744 
745 //=========================================================================
746 /// Eastern edge of the macro element in the collapsible part of the
747 /// channel; \f$ \zeta \in [-1,1] \f$
748 //=========================================================================
750  const Vector<double>& zeta,
751  Vector<double>& r,
752  const unsigned& imacro)
753 {
754 
755  //Determines the "coordinates" of the macro-element
756  unsigned x=unsigned(imacro%(Nup+Ncollapsible+Ndown));
757  unsigned y=unsigned(double(imacro)/double(Nup+Ncollapsible+Ndown));
758 
759  // Vector of Lagrangian coordinates
760  Vector<double> xi(1);
761  xi[0]=(double(x-Nup)+1.0)*(Lcollapsible/double(Ncollapsible));
762 
763  // Position vector on upper wall:
764  Vector<double> r_wall(2);
765  Wall_pt->position(t, xi,r_wall);
766 
767  // Point will be located on straight line from bottom to top wall
768  double fract=(double(y)+(0.5*(1.0+zeta[0])))/double(Ny);
769 
770  // Map it via squash fct
771  fract=s_squash(fract);
772 
773  // x-cooordinate -- straight line from fixed position on the bottom
774  // wall to moving position on the top wall
775  r[0]=Lup+xi[0]+(r_wall[0]-(xi[0]+Lup))*fract;
776 
777  // y-coordinate
778  r[1]=r_wall[1]*fract;
779 }
780 
781 
782 
783 //==========================================================================
784 /// Northern edge of the macro element in the collapsible part of the
785 /// channel; \f$ \zeta \in [-1,1] \f$
786 //==========================================================================
788  const Vector<double>& zeta,
789  Vector<double>& r,
790  const unsigned& imacro)
791 {
792 
793  //Determines the "coordinates" of the macro-element
794  unsigned x=unsigned(imacro%(Nup+Ncollapsible+Ndown));
795  unsigned y=unsigned(double(imacro)/double(Nup+Ncollapsible+Ndown));
796 
797  // Vector of Lagrangian coordinates
798  Vector<double> xi(1);
799  xi[0]=(double(x-Nup)+(0.5*(1.0+zeta[0])))*(Lcollapsible/double(Ncollapsible));
800 
801  // Position vector on upper wall:
802  Vector<double> r_wall(2);
803  Wall_pt->position(t,xi,r_wall);
804 
805  // Point will be located on straight line from bottom to top wall
806  double fract=(double(y)+1.0)/double(Ny);
807 
808  // Map it via squash fct
809  fract=s_squash(fract);
810 
811  // x-cooordinate -- straight line from fixed position on the bottom
812  // wall to moving position on the top wall
813  r[0]=Lup+xi[0]+(r_wall[0]-(xi[0]+Lup))*fract;
814 
815  // y-coordinate
816  r[1]=r_wall[1]*fract;
817 
818 }
819 
820 
821 //========================================================================
822 /// Southern edge of the macro element in the collapsible part of the
823 /// channel; \f$ \zeta \in [-1,1] \f$
824 //========================================================================
826  const Vector<double>& zeta,
827  Vector<double>& r,
828  const unsigned& imacro)
829 {
830 
831  //Determines the "coordinates" of the macro-element
832  unsigned x=unsigned(imacro%(Nup+Ncollapsible+Ndown));
833  unsigned y=unsigned(double(imacro)/double(Nup+Ncollapsible+Ndown));
834 
835  // Vector of Lagrangian coordinates
836  Vector<double> xi(1);
837  xi[0]=(double(x-Nup)+(0.5*(1.0+zeta[0])))*(Lcollapsible/double(Ncollapsible));
838 
839  // Position vector on upper wall:
840  Vector<double> r_wall(2);
841  Wall_pt->position(t,xi,r_wall);
842 
843  // Point will be located on straight line from bottom to top wall
844  double fract=double(y)/double(Ny);
845 
846  // Map it via squash fct
847  fract=s_squash(fract);
848 
849  // x-cooordinate -- straight line from fixed position on the bottom
850  // wall to moving position on the top wall
851  r[0]=Lup+xi[0]+(r_wall[0]-(xi[0]+Lup))*fract;
852 
853  // y-coordinate
854  r[1]=r_wall[1]*fract;
855 
856 
857 }
858 
859 
860 
861 
862 
863 
864 }
865 
866 #endif
void r_E_straight(const Vector< double > &zeta, Vector< double > &r, const unsigned &imacro, const unsigned &part)
Eastern boundary of the macro element imacro in the upstream (part=0) or downstream (part=1) sections...
unsigned Ny
Number of macro elements across channel.
double s_squash(const double &s)
Function that squashes the macro elements near the wall. Input argument should vary between 0 and 1; ...
double l_down()
Length of downstream section.
void r_S_collapsible(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &imacro)
Southern boundary of the macro element imacro in the collapsible section.
void r_W_straight(const Vector< double > &zeta, Vector< double > &r, const unsigned &imacro, const unsigned &part)
Western boundary of the macro element imacro in the upstream (part=0) or downstream (part=1) sections...
unsigned Ndown
Number of vertical element columns in downstream section.
unsigned ndown()
Number of vertical columns of macro elements in the downstream section.
CollapsibleChannelDomain(const unsigned &nup, const unsigned &ncollapsible, const unsigned &ndown, const unsigned &ny, const double &lup, const double &lcollapsible, const double &ldown, const double &ly, GeomObject *wall_pt)
Constructor: Pass the number of (macro-)elements, the domain lengths in the x- and y-direction and th...
double(* AxialSpacingFctPt)(const double &xi)
Typedef for function pointer for function that implements axial spacing of macro elements.
static double default_BL_squash_fct(const double &s)
Default for function that squashes the macro elements near the walls. Identity.
double l_collapsible()
Length of collapsible segment.
BLSquashFctPt BL_squash_fct_pt
Function pointer for function that squashes the macro elements near the walls.
void r_W_collapsible(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &imacro)
Western boundary of the macro element imacro in the collapsible section.
BLSquashFctPt & bl_squash_fct_pt()
Function pointer for function that squashes the macro elements near wall. Default mapping (identity) ...
void r_N_collapsible(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &imacro)
Northern boundary of the macro element imacro in the collapsible section.
double Lcollapsible
x-length in the "collapsible" part of the channel
unsigned Nup
Number of vertical element columns in upstream section.
unsigned ncollapsible()
Number of vertical clumns of macro elements in the "collapsible" segment.
unsigned Ncollapsible
Number of vertical element columns in "collapsible" section.
unsigned ny()
Number of macro-elements across the channel.
void macro_element_boundary(const unsigned &t, const unsigned &imacro, const unsigned &idirect, const Vector< double > &zeta, Vector< double > &r)
Vector representation of the imacro-th macro element boundary idirect (N/S/W/E) at time level t (t=0:...
double Lup
x-length in the upstream part of the channel
static double default_axial_spacing_fct(const double &xi)
Default for function that implements axial spacing of macro elements.
double(* BLSquashFctPt)(const double &s)
Typedef for function pointer for function that squashes the macro elements near the wall to help reso...
unsigned nup()
Number of vertical columns of macro elements the upstream section.
GeomObject * wall_pt() const
Access to pointer to the geometric object that parametrises the collapsible wall (const version) ...
AxialSpacingFctPt & axial_spacing_fct_pt()
Function pointer for function that implements axial spacing of macro elements.
void r_N_straight(const Vector< double > &zeta, Vector< double > &r, const unsigned &imacro, const unsigned &part)
Northern boundary of the macro element imacro in the upstream (part=0) or downstream (part=1) section...
double l_up()
Length of upstream section.
GeomObject *& wall_pt()
Access to pointer to the geometric object that parametrises the collapsible wall. ...
void r_E_collapsible(const unsigned &t, const Vector< double > &zeta, Vector< double > &r, const unsigned &imacro)
Eastern boundary of the macro element imacro in the collapsible section.
void r_S_straight(const Vector< double > &zeta, Vector< double > &r, const unsigned &imacro, const unsigned &part)
Southern boundary of the macro element imacro in the upstream (part=0) or downstream (part=1) section...
double Ldown
x-length in the downstream part of the channel
GeomObject * Wall_pt
Pointer to the geometric object that parametrises the collapsible wall.
AxialSpacingFctPt Axial_spacing_fct_pt
Function pointer for function that implements axial spacing of macro elements.
double axial_spacing_fct(const double &xi)
Function that implements axial spacing of macro elements.