collapsible_channel_mesh.template.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifndef OOMPH_COLLAPSIBLE_CHANNEL_MESH_TEMPLATE_CC
27 #define OOMPH_COLLAPSIBLE_CHANNEL_MESH_TEMPLATE_CC
28 
29 // Include the headers file for collapsible channel
31 
32 
33 namespace oomph
34 {
35  //========================================================================
36  /// Constructor: Pass number of elements in upstream/collapsible/downstream
37  /// segment and across the channel; lengths of upstream/collapsible/downstream
38  /// segments and width of channel, pointer to GeomObject that defines
39  /// the collapsible segment and pointer to TimeStepper (defaults to the
40  /// default timestepper, Steady).
41  //========================================================================
42  template<class ELEMENT>
44  const unsigned& nup,
45  const unsigned& ncollapsible,
46  const unsigned& ndown,
47  const unsigned& ny,
48  const double& lup,
49  const double& lcollapsible,
50  const double& ldown,
51  const double& ly,
52  GeomObject* wall_pt,
53  TimeStepper* time_stepper_pt)
54  : SimpleRectangularQuadMesh<ELEMENT>(nup + ncollapsible + ndown,
55  ny,
56  lup + lcollapsible + ldown,
57  ly,
58  time_stepper_pt),
59  Nup(nup),
60  Ncollapsible(ncollapsible),
61  Ndown(ndown),
62  Ny(ny),
63  Wall_pt(wall_pt)
64  {
65  // Mesh can only be built with 2D Qelements.
66  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
67 
68  // Create the CollapsibleChannelDomain with the wall represented
69  // by the geometric object
71  nup, ncollapsible, ndown, ny, lup, lcollapsible, ldown, ly, wall_pt);
72 
73  // Total number of (macro/finite)elements
74  unsigned nmacro = (nup + ncollapsible + ndown) * ny;
75 
76  // Loop over all elements and set macro element pointer
77  for (unsigned e = 0; e < nmacro; e++)
78  {
79  // Get pointer to finite element
80  FiniteElement* el_pt = this->finite_element_pt(e);
81 
82  // Set pointer to macro element so the curvlinear boundaries
83  // of the mesh/domain can get picked up during adaptive
84  // mesh refinement
85  el_pt->set_macro_elem_pt(this->Domain_pt->macro_element_pt(e));
86  }
87 
88  // Update the nodal positions corresponding to their
89  // macro element representations.
90  this->node_update();
91 
92  // Update the boundary numbering scheme and set boundary coordinate
93  //-----------------------------------------------------------------
94 
95  // (Note: The original SimpleRectangularQuadMesh had four boundaries.
96  // We need to overwrite the boundary lookup scheme for the current
97  // mesh so that the collapsible segment becomes identifiable).
98  // While we're doing this, we're also setting up a boundary
99  // coordinate for the nodes located on the collapsible segment.
100  // The boundary coordinate can be used to setup FSI.
101 
102  // How many boundaries does the mesh have now?
103  unsigned nbound = this->nboundary();
104  for (unsigned b = 0; b < nbound; b++)
105  {
106  // Remove all nodes on this boundary from the mesh's lookup scheme
107  // and also delete the reverse lookup scheme held by the nodes
108  this->remove_boundary_nodes(b);
109  }
110 
111 #ifdef PARANOID
112  // Sanity check
113  unsigned nnod = this->nnode();
114  for (unsigned j = 0; j < nnod; j++)
115  {
116  if (this->node_pt(j)->is_on_boundary())
117  {
118  std::ostringstream error_message;
119  error_message << "Node " << j << "is still on boundary " << std::endl;
120 
121  throw OomphLibError(error_message.str(),
122  OOMPH_CURRENT_FUNCTION,
123  OOMPH_EXCEPTION_LOCATION);
124  }
125  }
126 #endif
127 
128  // Change the numbers of boundaries
129  this->set_nboundary(6);
130 
131  // Get the number of nodes along the element edge from first element
132  unsigned nnode_1d = this->finite_element_pt(0)->nnode_1d();
133 
134  // Vector of Lagrangian coordinates used as boundary coordinate
135  Vector<double> zeta(1);
136 
137  // Index of first element underneath the collapsible bit
138  unsigned first_collapsible = (ny - 1) * (nup + ncollapsible + ndown) + nup;
139 
140  // Zeta increment over elements (used for assignment of
141  // boundary coordinate)
142  double dzeta = lcollapsible / double(ncollapsible);
143 
144  // Manually loop over the elements near the boundaries and
145  // assign nodes to boundaries. Also set up boundary coordinate
146  unsigned nelem = this->nelement();
147  for (unsigned e = 0; e < nelem; e++)
148  {
149  // Bottom row of elements
150  if (e < nup + ncollapsible + ndown)
151  {
152  for (unsigned i = 0; i < nnode_1d; i++)
153  {
154  this->add_boundary_node(0, this->finite_element_pt(e)->node_pt(i));
155  }
156  }
157  // Upstream upper rigid bit
158  if ((e > (ny - 1) * (nup + ncollapsible + ndown) - 1) &&
159  (e < (ny - 1) * (nup + ncollapsible + ndown) + nup))
160  {
161  for (unsigned i = 0; i < nnode_1d; i++)
162  {
163  this->add_boundary_node(
164  4,
165  this->finite_element_pt(e)->node_pt((nnode_1d - 1) * nnode_1d + i));
166  }
167  }
168  // Collapsible bit
169  if ((e > (ny - 1) * (nup + ncollapsible + ndown) + nup - 1) &&
170  (e < (ny - 1) * (nup + ncollapsible + ndown) + nup + ncollapsible))
171  {
172  for (unsigned i = 0; i < nnode_1d; i++)
173  {
174  this->add_boundary_node(
175  3,
176  this->finite_element_pt(e)->node_pt((nnode_1d - 1) * nnode_1d + i));
177 
178  // What column of elements are we in?
179  unsigned ix = e - first_collapsible;
180 
181  // Zeta coordinate
182  zeta[0] =
183  double(ix) * dzeta + double(i) * dzeta / double(nnode_1d - 1);
184 
185  // Set boundary coordinate
186  this->finite_element_pt(e)
187  ->node_pt((nnode_1d - 1) * nnode_1d + i)
188  ->set_coordinates_on_boundary(3, zeta);
189  }
190  }
191  // Downstream upper rigid bit
192  if ((e >
193  (ny - 1) * (nup + ncollapsible + ndown) + nup + ncollapsible - 1) &&
194  (e < ny * (nup + ncollapsible + ndown)))
195  {
196  for (unsigned i = 0; i < nnode_1d; i++)
197  {
198  this->add_boundary_node(
199  2,
200  this->finite_element_pt(e)->node_pt((nnode_1d - 1) * nnode_1d + i));
201  }
202  }
203  // Left end
204  if (e % (nup + ncollapsible + ndown) == 0)
205  {
206  for (unsigned i = 0; i < nnode_1d; i++)
207  {
208  this->add_boundary_node(
209  5, this->finite_element_pt(e)->node_pt(i * nnode_1d));
210  }
211  }
212  // Right end
213  if (e % (nup + ncollapsible + ndown) == (nup + ncollapsible + ndown) - 1)
214  {
215  for (unsigned i = 0; i < nnode_1d; i++)
216  {
217  this->add_boundary_node(
218  1, this->finite_element_pt(e)->node_pt((i + 1) * nnode_1d - 1));
219  }
220  }
221  }
222 
223  // Re-setup lookup scheme that establishes which elements are located
224  // on the mesh boundaries (doesn't need to be wiped)
225  this->setup_boundary_element_info();
226 
227  // We have only bothered to parametrise boundary 3
228  this->Boundary_coordinate_exists[3] = true;
229  }
230 
231 
232  /// ////////////////////////////////////////////////////////////////////////
233  /// ////////////////////////////////////////////////////////////////////////
234  /// ////////////////////////////////////////////////////////////////////////
235 
236 
237  //=================================================================
238  /// Perform algebraic mesh update at time level t (t=0: present;
239  /// t>0: previous)
240  //=================================================================
241  template<class ELEMENT>
243  const unsigned& t, AlgebraicNode*& node_pt)
244  {
245 #ifdef PARANOID
246  // We're updating the nodal positions (!) at time level t
247  // and determine them by evaluating the wall GeomObject's
248  // position at that gime level. I believe this only makes sense
249  // if the t-th history value in the positional timestepper
250  // actually represents previous values (rather than some
251  // generalised quantity). Hence if this function is called with
252  // t>nprev_values(), we issue a warning and terminate the execution.
253  // It *might* be possible that the code still works correctly
254  // even if this condition is violated (e.g. if the GeomObject's
255  // position() function returns the appropriate "generalised"
256  // position value that is required by the timestepping scheme but it's
257  // probably worth flagging this up and forcing the user to manually switch
258  // off this warning if he/she is 100% sure that this is kosher.
259  if (t > node_pt->position_time_stepper_pt()->nprev_values())
260  {
261  std::string error_message =
262  "Trying to update the nodal position at a time level";
263  error_message += "beyond the number of previous values in the nodes'";
264  error_message += "position timestepper. This seems highly suspect!";
265  error_message += "If you're sure the code behaves correctly";
266  error_message += "in your application, remove this warning ";
267  error_message += "or recompile with PARNOID switched off.";
268 
269  std::string function_name = "AlgebraicCollapsibleChannelMesh::";
270  function_name += "algebraic_node_update()";
271 
272  throw OomphLibError(
273  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
274  }
275 #endif
276 
277  // Extract references for update by copy construction
278  Vector<double> ref_value(node_pt->vector_ref_value());
279 
280  // First reference value: Original x-position. Used as the start point
281  // for the lines connecting the nodes in the vertical direction
282  double x_bottom = ref_value[0];
283 
284  // Second reference value: Fractional position along
285  // straight line from the bottom (at the original x position)
286  // to the reference point on the upper wall
287  double fract = ref_value[1];
288 
289  // Third reference value: Reference local coordinate
290  // in GeomObject that represents the upper wall (local coordinate
291  // in finite element if the wall GeomObject is a finite element mesh)
292  Vector<double> s(1);
293  s[0] = ref_value[2];
294 
295  // Fourth reference value: zeta coordinate on the upper wall
296  // If the wall is a simple GeomObject, zeta[0]=s[0]
297  // but if it's a compound GeomObject (e.g. a finite element mesh)
298  // zeta scales during mesh refinement, whereas s[0] and the
299  // pointer to the geom object have to be re-computed.
300  // double zeta=ref_value[3]; // not needed here
301 
302  // Extract geometric objects for update by copy construction
303  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
304 
305  // Pointer to actual wall geom object (either the same as the wall object
306  // or the pointer to the actual finite element)
307  GeomObject* geom_obj_pt = geom_object_pt[0];
308 
309  // Get position vector to wall at previous timestep t!
310  Vector<double> r_wall(2);
311  geom_obj_pt->position(t, s, r_wall);
312 
313  // Assign new nodal coordinate
314  node_pt->x(t, 0) = x_bottom + fract * (r_wall[0] - x_bottom);
315  node_pt->x(t, 1) = fract * r_wall[1];
316  }
317 
318 
319  //=====start_setup=================================================
320  /// Setup algebraic mesh update -- assumes that mesh has
321  /// initially been set up with a flush upper wall
322  //=================================================================
323  template<class ELEMENT>
325  {
326  // Shorthand for some geometric data:
327  double l_up = this->domain_pt()->l_up();
328  double l_collapsible = this->domain_pt()->l_collapsible();
329 
330  // Loop over all nodes in mesh
331  unsigned nnod = this->nnode();
332  for (unsigned j = 0; j < nnod; j++)
333  {
334  // Get pointer to node -- recall that that Mesh::node_pt(...) has been
335  // overloaded in the AlgebraicMesh class to return a pointer to
336  // an AlgebraicNode.
337  AlgebraicNode* nod_pt = node_pt(j);
338 
339  // Get coordinates
340  double x = nod_pt->x(0);
341  double y = nod_pt->x(1);
342 
343  // Check if it's in the collapsible part:
344  if ((x >= l_up) && (x <= (l_up + l_collapsible)))
345  {
346  // Get zeta coordinate on the undeformed wall
347  Vector<double> zeta(1);
348  zeta[0] = x - l_up;
349 
350  // Get pointer to geometric (sub-)object and Lagrangian coordinate
351  // on that sub-object. For a wall that is represented by
352  // a single geom object, this simply returns the input.
353  // If the geom object consists of sub-objects (e.g.
354  // if it is a finite element mesh representing a wall,
355  // then we'll obtain the pointer to the finite element
356  // (in its incarnation as a GeomObject) and the
357  // local coordinate in that element.
358  GeomObject* geom_obj_pt;
359  Vector<double> s(1);
360  this->Wall_pt->locate_zeta(zeta, geom_obj_pt, s);
361 
362  // Get position vector to wall:
363  Vector<double> r_wall(2);
364  geom_obj_pt->position(s, r_wall);
365 
366  // Sanity check: Confirm that the wall is in its undeformed position
367 #ifdef PARANOID
368  if ((std::fabs(r_wall[0] - x) > 1.0e-15) &&
369  (std::fabs(r_wall[1] - y) > 1.0e-15))
370  {
371  std::ostringstream error_stream;
372  error_stream << "Wall must be in its undeformed position when\n"
373  << "algebraic node update information is set up!\n "
374  << "x-discrepancy: " << std::fabs(r_wall[0] - x)
375  << std::endl
376  << "y-discrepancy: " << std::fabs(r_wall[1] - y)
377  << std::endl;
378 
379  throw OomphLibError(error_stream.str(),
380  OOMPH_CURRENT_FUNCTION,
381  OOMPH_EXCEPTION_LOCATION);
382  }
383 #endif
384 
385 
386  // One geometric object is involved in update operation
387  Vector<GeomObject*> geom_object_pt(1);
388 
389  // The actual geometric object (If the wall is simple GeomObject
390  // this is the same as Wall_pt; if it's a compound GeomObject
391  // this points to the sub-object)
392  geom_object_pt[0] = geom_obj_pt;
393 
394  // The update function requires four parameters:
395  Vector<double> ref_value(4);
396 
397  // First reference value: Original x-position
398  ref_value[0] = r_wall[0];
399 
400  // Second reference value: fractional position along
401  // straight line from the bottom (at the original x position)
402  // to the point on the wall)
403  ref_value[1] = y / r_wall[1];
404 
405  // Third reference value: Reference local coordinate
406  // in wall element (local coordinate in FE if we're dealing
407  // with a wall mesh)
408  ref_value[2] = s[0];
409 
410  // Fourth reference value: zeta coordinate on wall
411  // If the wall is a simple GeomObject, zeta[0]=s[0]
412  // but if it's a compound GeomObject (e.g. a finite element mesh)
413  // zeta scales during mesh refinement, whereas s[0] and the
414  // pointer to the geom object have to be re-computed.
415  ref_value[3] = zeta[0];
416 
417  // Setup algebraic update for node: Pass update information
418  nod_pt->add_node_update_info(this, // mesh
419  geom_object_pt, // vector of geom objects
420  ref_value); // vector of ref. values
421  }
422  }
423 
424  } // end of setup_algebraic_node_update
425 
426 
427  /// /////////////////////////////////////////////////////////////////
428  /// /////////////////////////////////////////////////////////////////
429  /// /////////////////////////////////////////////////////////////////
430 
431 
432  //========start_update_node_update=================================
433  /// Update the geometric references that are used
434  /// to update node after mesh adaptation.
435  //=================================================================
436  template<class ELEMENT>
438  AlgebraicNode*& node_pt)
439  {
440  // Extract reference values for node update by copy construction
441  Vector<double> ref_value(node_pt->vector_ref_value());
442 
443  // First reference value: Original x-position
444  // double x_bottom=ref_value[0]; // not needed here
445 
446  // Second reference value: fractional position along
447  // straight line from the bottom (at the original x position)
448  // to the point on the wall)
449  // double fract=ref_value[1]; // not needed here
450 
451  // Third reference value: Reference local coordinate
452  // in GeomObject (local coordinate in finite element if the wall
453  // GeomObject is a finite element mesh)
454  // Vector<double> s(1);
455  // s[0]=ref_value[2]; // This needs to be re-computed!
456 
457  // Fourth reference value: intrinsic coordinate on the (possibly
458  // compound) wall.
459  double zeta = ref_value[3];
460 
461  // Extract geometric objects for update by copy construction
462  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
463 
464  // Pointer to actual wall geom object (either the same as wall object
465  // or the pointer to the actual finite element)
466  // GeomObject* geom_obj_pt=geom_object_pt[0]; // This needs to be
467  // re-computed!
468 
469  // Get zeta coordinate on wall (as vector)
470  Vector<double> zeta_wall(1);
471  zeta_wall[0] = zeta;
472 
473  // Get pointer to geometric (sub-)object and Lagrangian coordinate
474  // on that sub-object. For a wall that is represented by
475  // a single geom object, this simply returns the input.
476  // If the geom object consists of sub-objects (e.g.
477  // if it is a finite element mesh representing a wall,
478  // then we'll obtain the pointer to the finite element
479  // (in its incarnation as a GeomObject) and the
480  // local coordinate in that element.
481  Vector<double> s(1);
482  GeomObject* geom_obj_pt;
483  this->Wall_pt->locate_zeta(zeta_wall, geom_obj_pt, s);
484 
485  // Update the pointer to the (sub-)GeomObject within which the
486  // reference point is located. (If the wall is simple GeomObject
487  // this is the same as Wall_pt; if it's a compound GeomObject
488  // this points to the sub-object)
489  geom_object_pt[0] = geom_obj_pt;
490 
491  // First reference value: Original x-position
492  // ref_value[0]=r_wall[0]; // unchanged
493 
494  // Second reference value: fractional position along
495  // straight line from the bottom (at the original x position)
496  // to the point on the wall)
497  // ref_value[1]=y/r_wall[1]; // unchanged
498 
499  // Update third reference value: Reference local coordinate
500  // in wall element (local coordinate in FE if we're dealing
501  // with a wall mesh)
502  ref_value[2] = s[0];
503 
504  // Fourth reference value: zeta coordinate on wall
505  // If the wall is a simple GeomObject, zeta[0]=s[0]
506  // but if it's a compound GeomObject (e.g. a finite element mesh)
507  // zeta scales during mesh refinement, whereas s[0] and the
508  // pointer to the geom object have to be re-computed.
509  // ref_value[3]=zeta[0]; //unchanged
510 
511 
512  // Kill the existing node update info
513  node_pt->kill_node_update_info();
514 
515  // Setup algebraic update for node: Pass update information
516  node_pt->add_node_update_info(this, // mesh
517  geom_object_pt, // vector of geom objects
518  ref_value); // vector of ref. values
519  }
520 
521 
522 } // namespace oomph
523 #endif
void setup_algebraic_node_update()
Function to setup the algebraic node update.
void algebraic_node_update(const unsigned &t, AlgebraicNode *&node_pt)
Update nodal position at time level t (t=0: present; t>0: previous)
GeomObject *& wall_pt()
Access function to GeomObject representing wall.
CollapsibleChannelMesh(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, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in upstream/collapsible/ downstream segment and across the chann...
CollapsibleChannelDomain * Domain_pt
Pointer to domain.
void update_node_update(AlgebraicNode *&node_pt)
Update the node update data for specified node following any mesh adapation.
Simple rectangular 2D Quad mesh class. Nx : number of elements in the x direction.
const unsigned & ny() const
Access function for number of elements in y directions.
////////////////////////////////////////////////////////////////////// //////////////////////////////...