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