fsi_driven_cavity_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_FSI_DRIVEN_CAVITY_MESH_TEMPLATE_CC
27 #define OOMPH_FSI_DRIVEN_CAVITY_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, lengths, pointer to GeomObject
37  /// that defines the collapsible segment and pointer to TimeStepper
38  /// (defaults to the default timestepper, Steady).
39  //========================================================================
40  template<class ELEMENT>
42  const unsigned& nx,
43  const unsigned& ny,
44  const double& lx,
45  const double& ly,
46  const double& gap_fraction,
47  GeomObject* wall_pt,
48  TimeStepper* time_stepper_pt)
49  : SimpleRectangularQuadMesh<ELEMENT>(nx, ny, lx, ly, time_stepper_pt),
50  Nx(nx),
51  Ny(ny),
52  Gap_fraction(gap_fraction),
53  Wall_pt(wall_pt)
54  {
55  // Mesh can only be built with 2D Qelements.
56  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
57 
58  // Update the boundary numbering scheme and set boundary coordinate
59  //-----------------------------------------------------------------
60 
61  // (Note: The original SimpleRectangularQuadMesh had four boundaries.
62  // We need to overwrite the boundary lookup scheme for the current
63  // mesh so that the collapsible segment becomes identifiable).
64  // While we're doing this, we're also setting up a boundary
65  // coordinate for the nodes located on the collapsible segment.
66  // The boundary coordinate can be used to setup FSI.
67 
68  // How many boundaries does the mesh have now?
69  unsigned nbound = this->nboundary();
70  for (unsigned b = 0; b < nbound; b++)
71  {
72  // Remove all nodes on this boundary from the mesh's lookup scheme
73  // and also delete the reverse lookup scheme held by the nodes
74  this->remove_boundary_nodes(b);
75  }
76 
77 #ifdef PARANOID
78  // Sanity check
79  unsigned nnod = this->nnode();
80  for (unsigned j = 0; j < nnod; j++)
81  {
82  if (this->node_pt(j)->is_on_boundary())
83  {
84  std::ostringstream error_message;
85  error_message << "Node " << j << "is still on boundary " << std::endl;
86 
87  throw OomphLibError(error_message.str(),
88  OOMPH_CURRENT_FUNCTION,
89  OOMPH_EXCEPTION_LOCATION);
90  }
91  }
92 #endif
93 
94  // Change the numbers of boundaries
95  this->set_nboundary(6);
96 
97  // Get the number of nodes along the element edge from first element
98  unsigned nnode_1d = this->finite_element_pt(0)->nnode_1d();
99 
100  // Vector of Lagrangian coordinates used as boundary coordinate
101  Vector<double> zeta(1);
102 
103  // Zeta increment over elements (used for assignment of
104  // boundary coordinate)
105  double dzeta = lx / double(nx);
106 
107  // Manually loop over the elements near the boundaries and
108  // assign nodes to boundaries. Also set up boundary coordinate
109  unsigned nelem = this->nelement();
110  for (unsigned e = 0; e < nelem; e++)
111  {
112  // Bottom row of elements
113  if (e < nx)
114  {
115  for (unsigned i = 0; i < nnode_1d; i++)
116  {
117  this->add_boundary_node(0, this->finite_element_pt(e)->node_pt(i));
118  }
119  }
120  // Collapsible bit
121  if (e > ((ny - 1) * nx) - 1)
122  {
123  for (unsigned i = 0; i < nnode_1d; i++)
124  {
125  this->add_boundary_node(
126  3, this->finite_element_pt(e)->node_pt(2 * nnode_1d + i));
127 
128  // What column of elements are we in?
129  unsigned ix = e - (ny - 1) * nx;
130 
131  // Zeta coordinate
132  zeta[0] =
133  double(ix) * dzeta + double(i) * dzeta / double(nnode_1d - 1);
134 
135  // Set boundary coordinate
136  this->finite_element_pt(e)
137  ->node_pt(2 * nnode_1d + i)
138  ->set_coordinates_on_boundary(3, zeta);
139  }
140  }
141  // Left end
142  if (e % (nx) == 0)
143  {
144  for (unsigned i = 0; i < nnode_1d; i++)
145  {
146  Node* nod_pt = this->finite_element_pt(e)->node_pt(i * nnode_1d);
147 
148  // Rigid bit?
149  if (nod_pt->x(1) >= ly * Gap_fraction)
150  {
151  this->add_boundary_node(4, nod_pt);
152  }
153  // Free bit
154  else
155  {
156  this->add_boundary_node(5, nod_pt);
157  }
158  }
159  }
160  // Right end
161  if (e % nx == nx - 1)
162  {
163  for (unsigned i = 0; i < nnode_1d; i++)
164  {
165  Node* nod_pt =
166  this->finite_element_pt(e)->node_pt((i + 1) * nnode_1d - 1);
167 
168  // Rigid bit?
169  if (nod_pt->x(1) >= ly * Gap_fraction)
170  {
171  this->add_boundary_node(2, nod_pt);
172  }
173  // Free bit
174  else
175  {
176  this->add_boundary_node(1, nod_pt);
177  }
178  }
179  }
180  }
181 
182  // Re-setup lookup scheme that establishes which elements are located
183  // on the mesh boundaries (doesn't need to be wiped)
184  this->setup_boundary_element_info();
185 
186  // We have only bothered to parametrise boundary 3
187  this->Boundary_coordinate_exists[3] = true;
188  }
189 
190 
191  /// ////////////////////////////////////////////////////////////////////////
192  /// ////////////////////////////////////////////////////////////////////////
193  /// ////////////////////////////////////////////////////////////////////////
194 
195 
196  //=================================================================
197  /// Perform algebraic mesh update at time level t (t=0: present;
198  /// t>0: previous)
199  //=================================================================
200  template<class ELEMENT>
202  const unsigned& t, AlgebraicNode*& node_pt)
203  {
204 #ifdef PARANOID
205  // We're updating the nodal positions (!) at time level t
206  // and determine them by evaluating the wall GeomObject's
207  // position at that gime level. I believe this only makes sense
208  // if the t-th history value in the positional timestepper
209  // actually represents previous values (rather than some
210  // generalised quantity). Hence if this function is called with
211  // t>nprev_values(), we issue a warning and terminate the execution.
212  // It *might* be possible that the code still works correctly
213  // even if this condition is violated (e.g. if the GeomObject's
214  // position() function returns the appropriate "generalised"
215  // position value that is required by the timestepping scheme but it's
216  // probably worth flagging this up and forcing the user to manually switch
217  // off this warning if he/she is 100% sure that this is kosher.
218  if (t > node_pt->position_time_stepper_pt()->nprev_values())
219  {
220  std::string error_message =
221  "Trying to update the nodal position at a time level";
222  error_message += "beyond the number of previous values in the nodes'";
223  error_message += "position timestepper. This seems highly suspect!";
224  error_message += "If you're sure the code behaves correctly";
225  error_message += "in your application, remove this warning ";
226  error_message += "or recompile with PARNOID switched off.";
227 
228  std::string function_name = "AlgebraicFSIDrivenCavityMesh::";
229  function_name += "algebraic_node_update()";
230 
231  throw OomphLibError(
232  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
233  }
234 #endif
235 
236  // Extract references for update by copy construction
237  Vector<double> ref_value(node_pt->vector_ref_value());
238 
239  // First reference value: Original x-position. Used as the start point
240  // for the lines connecting the nodes in the vertical direction
241  double x_bottom = ref_value[0];
242 
243  // Second reference value: Fractional position along
244  // straight line from the bottom (at the original x position)
245  // to the reference point on the upper wall
246  double fract = ref_value[1];
247 
248  // Third reference value: Reference local coordinate
249  // in GeomObject that represents the upper wall (local coordinate
250  // in finite element if the wall GeomObject is a finite element mesh)
251  Vector<double> s(1);
252  s[0] = ref_value[2];
253 
254  // Fourth reference value: zeta coordinate on the upper wall
255  // If the wall is a simple GeomObject, zeta[0]=s[0]
256  // but if it's a compound GeomObject (e.g. a finite element mesh)
257  // zeta scales during mesh refinement, whereas s[0] and the
258  // pointer to the geom object have to be re-computed.
259  // double zeta=ref_value[3]; // not needed here
260 
261  // Extract geometric objects for update by copy construction
262  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
263 
264  // Pointer to actual wall geom object (either the same as the wall object
265  // or the pointer to the actual finite element)
266  GeomObject* geom_obj_pt = geom_object_pt[0];
267 
268  // Get position vector to wall at previous timestep t!
269  Vector<double> r_wall(2);
270  geom_obj_pt->position(t, s, r_wall);
271 
272  // Assign new nodal coordinate
273  node_pt->x(t, 0) = x_bottom + fract * (r_wall[0] - x_bottom);
274  node_pt->x(t, 1) = fract * r_wall[1];
275  }
276 
277 
278  //=====start_setup=================================================
279  /// Setup algebraic mesh update -- assumes that mesh has
280  /// initially been set up with a flush upper wall.
281  //=================================================================
282  template<class ELEMENT>
284  {
285  // Loop over all nodes in mesh
286  unsigned nnod = this->nnode();
287  for (unsigned j = 0; j < nnod; j++)
288  {
289  // Get pointer to node -- recall that that Mesh::node_pt(...) has been
290  // overloaded in the AlgebraicMesh class to return a pointer to
291  // an AlgebraicNode.
292  AlgebraicNode* nod_pt = node_pt(j);
293 
294  // Get coordinates
295  double x = nod_pt->x(0);
296  double y = nod_pt->x(1);
297 
298  // Get zeta coordinate on the undeformed wall
299  Vector<double> zeta(1);
300  zeta[0] = x;
301 
302  // Get pointer to geometric (sub-)object and Lagrangian coordinate
303  // on that sub-object. For a wall that is represented by
304  // a single geom object, this simply returns the input.
305  // If the geom object consists of sub-objects (e.g.
306  // if it is a finite element mesh representing a wall,
307  // then we'll obtain the pointer to the finite element
308  // (in its incarnation as a GeomObject) and the
309  // local coordinate in that element.
310  GeomObject* geom_obj_pt;
311  Vector<double> s(1);
312  this->Wall_pt->locate_zeta(zeta, geom_obj_pt, s);
313 
314  // Get position vector to wall:
315  Vector<double> r_wall(2);
316  geom_obj_pt->position(s, r_wall);
317 
318  // Sanity check: Confirm that the wall is in its undeformed position
319 #ifdef PARANOID
320  if ((std::fabs(r_wall[0] - x) > 1.0e-15) &&
321  (std::fabs(r_wall[1] - y) > 1.0e-15))
322  {
323  std::ostringstream error_stream;
324  error_stream << "Wall must be in its undeformed position when\n"
325  << "algebraic node update information is set up!\n "
326  << "x-discrepancy: " << std::fabs(r_wall[0] - x)
327  << std::endl
328  << "y-discrepancy: " << std::fabs(r_wall[1] - y)
329  << std::endl;
330 
331  throw OomphLibError(
332  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
333  }
334 #endif
335 
336 
337  // One geometric object is involved in update operation
338  Vector<GeomObject*> geom_object_pt(1);
339 
340  // The actual geometric object (If the wall is simple GeomObject
341  // this is the same as Wall_pt; if it's a compound GeomObject
342  // this points to the sub-object)
343  geom_object_pt[0] = geom_obj_pt;
344 
345  // The update function requires four parameters:
346  Vector<double> ref_value(4);
347 
348  // First reference value: Original x-position
349  ref_value[0] = r_wall[0];
350 
351  // Second reference value: fractional position along
352  // straight line from the bottom (at the original x position)
353  // to the point on the wall)
354  ref_value[1] = y / r_wall[1];
355 
356  // Third reference value: Reference local coordinate
357  // in wall element (local coordinate in FE if we're dealing
358  // with a wall mesh)
359  ref_value[2] = s[0];
360 
361  // Fourth reference value: zeta coordinate on wall
362  // If the wall is a simple GeomObject, zeta[0]=s[0]
363  // but if it's a compound GeomObject (e.g. a finite element mesh)
364  // zeta scales during mesh refinement, whereas s[0] and the
365  // pointer to the geom object have to be re-computed.
366  ref_value[3] = zeta[0];
367 
368  // Setup algebraic update for node: Pass update information
369  nod_pt->add_node_update_info(this, // mesh
370  geom_object_pt, // vector of geom objects
371  ref_value); // vector of ref. values
372  }
373 
374 
375  } // end of setup_algebraic_node_update
376 
377 
378  /// /////////////////////////////////////////////////////////////////
379  /// /////////////////////////////////////////////////////////////////
380  /// /////////////////////////////////////////////////////////////////
381 
382 
383  //========start_update_node_update=================================
384  /// Update the geometric references that are used
385  /// to update node after mesh adaptation.
386  //=================================================================
387  template<class ELEMENT>
389  AlgebraicNode*& node_pt)
390  {
391  // Extract reference values for node update by copy construction
392  Vector<double> ref_value(node_pt->vector_ref_value());
393 
394  // First reference value: Original x-position
395  // double x_bottom=ref_value[0]; // not needed here
396 
397  // Second reference value: fractional position along
398  // straight line from the bottom (at the original x position)
399  // to the point on the wall)
400  // double fract=ref_value[1]; // not needed here
401 
402  // Third reference value: Reference local coordinate
403  // in GeomObject (local coordinate in finite element if the wall
404  // GeomObject is a finite element mesh)
405  // Vector<double> s(1);
406  // s[0]=ref_value[2]; // This needs to be re-computed!
407 
408  // Fourth reference value: intrinsic coordinate on the (possibly
409  // compound) wall.
410  double zeta = ref_value[3];
411 
412  // Extract geometric objects for update by copy construction
413  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
414 
415  // Pointer to actual wall geom object (either the same as wall object
416  // or the pointer to the actual finite element)
417  // GeomObject* geom_obj_pt=geom_object_pt[0]; // This needs to be
418  // re-computed!
419 
420  // Get zeta coordinate on wall (as vector)
421  Vector<double> zeta_wall(1);
422  zeta_wall[0] = zeta;
423 
424  // Get pointer to geometric (sub-)object and Lagrangian coordinate
425  // on that sub-object. For a wall that is represented by
426  // a single geom object, this simply returns the input.
427  // If the geom object consists of sub-objects (e.g.
428  // if it is a finite element mesh representing a wall,
429  // then we'll obtain the pointer to the finite element
430  // (in its incarnation as a GeomObject) and the
431  // local coordinate in that element.
432  Vector<double> s(1);
433  GeomObject* geom_obj_pt;
434  this->Wall_pt->locate_zeta(zeta_wall, geom_obj_pt, s);
435 
436  // Update the pointer to the (sub-)GeomObject within which the
437  // reference point is located. (If the wall is simple GeomObject
438  // this is the same as Wall_pt; if it's a compound GeomObject
439  // this points to the sub-object)
440  geom_object_pt[0] = geom_obj_pt;
441 
442  // First reference value: Original x-position
443  // ref_value[0]=r_wall[0]; // unchanged
444 
445  // Second reference value: fractional position along
446  // straight line from the bottom (at the original x position)
447  // to the point on the wall)
448  // ref_value[1]=y/r_wall[1]; // unchanged
449 
450  // Update third reference value: Reference local coordinate
451  // in wall element (local coordinate in FE if we're dealing
452  // with a wall mesh)
453  ref_value[2] = s[0];
454 
455  // Fourth reference value: zeta coordinate on wall
456  // If the wall is a simple GeomObject, zeta[0]=s[0]
457  // but if it's a compound GeomObject (e.g. a finite element mesh)
458  // zeta scales during mesh refinement, whereas s[0] and the
459  // pointer to the geom object have to be re-computed.
460  // ref_value[3]=zeta[0]; //unchanged
461 
462 
463  // Kill the existing node update info
464  node_pt->kill_node_update_info();
465 
466  // Setup algebraic update for node: Pass update information
467  node_pt->add_node_update_info(this, // mesh
468  geom_object_pt, // vector of geom objects
469  ref_value); // vector of ref. values
470  }
471 
472 
473 } // namespace oomph
474 #endif
void algebraic_node_update(const unsigned &t, AlgebraicNode *&node_pt)
Update nodal position at time level t (t=0: present; t>0: previous)
void setup_algebraic_node_update()
Function to setup the algebraic node update.
double Gap_fraction
Fraction of the gap next to moving lid, relative to the height of the domain.
FSIDrivenCavityMesh(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly, const double &gap_fraction, GeomObject *wall_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements, number of elements, fractional height of the gap above the movi...
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.
const unsigned & nx() const
Access function for number of elements in x directions.
////////////////////////////////////////////////////////////////////// //////////////////////////////...