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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifndef OOMPH_FSI_DRIVEN_CAVITY_MESH_TEMPLATE_CC
31 #define OOMPH_FSI_DRIVEN_CAVITY_MESH_TEMPLATE_CC
32 
33 //Include the headers file for collapsible channel
35 
36 
37 namespace oomph
38 {
39 
40 //========================================================================
41 /// Constructor: Pass number of elements, lengths, pointer to GeomObject
42 /// that defines the collapsible segment and pointer to TimeStepper
43 /// (defaults to the default timestepper, Steady).
44 //========================================================================
45 template <class ELEMENT>
47  const unsigned& nx,
48  const unsigned& ny,
49  const double& lx,
50  const double& ly,
51  const double& gap_fraction,
52  GeomObject* wall_pt,
53  TimeStepper* time_stepper_pt)
54  : SimpleRectangularQuadMesh<ELEMENT>(nx,ny,lx,ly,time_stepper_pt),
55  Nx(nx), Ny(ny), Gap_fraction(gap_fraction), Wall_pt(wall_pt)
56 {
57  // Mesh can only be built with 2D Qelements.
58  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
59 
60  // Update the boundary numbering scheme and set boundary coordinate
61  //-----------------------------------------------------------------
62 
63  // (Note: The original SimpleRectangularQuadMesh had four boundaries.
64  // We need to overwrite the boundary lookup scheme for the current
65  // mesh so that the collapsible segment becomes identifiable).
66  // While we're doing this, we're also setting up a boundary
67  // coordinate for the nodes located on the collapsible segment.
68  // The boundary coordinate can be used to setup FSI.
69 
70  // How many boundaries does the mesh have now?
71  unsigned nbound=this->nboundary();
72  for (unsigned b=0;b<nbound;b++)
73  {
74  // Remove all nodes on this boundary from the mesh's lookup scheme
75  // and also delete the reverse lookup scheme held by the nodes
76  this->remove_boundary_nodes(b);
77  }
78 
79 #ifdef PARANOID
80  // Sanity check
81  unsigned nnod=this->nnode();
82  for (unsigned j=0;j<nnod;j++)
83  {
84  if (this->node_pt(j)->is_on_boundary())
85  {
86  std::ostringstream error_message;
87  error_message << "Node " << j << "is still on boundary " << std::endl;
88 
89  throw OomphLibError(error_message.str(),
90  OOMPH_CURRENT_FUNCTION,
91  OOMPH_EXCEPTION_LOCATION);
92  }
93  }
94 #endif
95 
96  //Change the numbers of boundaries
97  this->set_nboundary(6);
98 
99  // Get the number of nodes along the element edge from first element
100  unsigned nnode_1d=this->finite_element_pt(0)->nnode_1d();
101 
102  // Vector of Lagrangian coordinates used as boundary coordinate
103  Vector<double> zeta(1);
104 
105  // Zeta increment over elements (used for assignment of
106  // boundary coordinate)
107  double dzeta= lx/double(nx);
108 
109  // Manually loop over the elements near the boundaries and
110  // assign nodes to boundaries. Also set up boundary coordinate
111  unsigned nelem=this->nelement();
112  for (unsigned e=0;e<nelem;e++)
113  {
114  // Bottom row of elements
115  if (e<nx)
116  {
117  for (unsigned i=0;i<nnode_1d;i++)
118  {
119  this->add_boundary_node(0, this->finite_element_pt(e)->node_pt(i));
120  }
121  }
122  // Collapsible bit
123  if (e>((ny-1)*nx)-1)
124  {
125  for (unsigned i=0;i<nnode_1d;i++)
126  {
127  this->add_boundary_node(3,
128  this->finite_element_pt(e)->node_pt(2*nnode_1d+i));
129 
130  // What column of elements are we in?
131  unsigned ix=e-(ny-1)*nx;
132 
133  // Zeta coordinate
134  zeta[0]=double(ix)*dzeta+double(i)*dzeta/double(nnode_1d-1);
135 
136  // Set boundary coordinate
137  this->finite_element_pt(e)->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=this->finite_element_pt(e)->node_pt((i+1)*nnode_1d-1);
166 
167  // Rigid bit?
168  if (nod_pt->x(1)>=ly*Gap_fraction)
169  {
170  this->add_boundary_node(2,nod_pt);
171  }
172  // Free bit
173  else
174  {
175  this->add_boundary_node(1,nod_pt);
176  }
177  }
178  }
179  }
180 
181  // Re-setup lookup scheme that establishes which elements are located
182  // on the mesh boundaries (doesn't need to be wiped)
184 
185  //We have only bothered to parametrise boundary 3
186  this->Boundary_coordinate_exists[3] = true;
187 }
188 
189 
190 
191 
192 ///////////////////////////////////////////////////////////////////////////
193 ///////////////////////////////////////////////////////////////////////////
194 ///////////////////////////////////////////////////////////////////////////
195 
196 
197 
198 //=================================================================
199 /// Perform algebraic mesh update at time level t (t=0: present;
200 /// t>0: previous)
201 //=================================================================
202 template<class ELEMENT>
204  const unsigned& t, AlgebraicNode*& node_pt)
205 {
206 
207 
208 #ifdef PARANOID
209  // We're updating the nodal positions (!) at time level t
210  // and determine them by evaluating the wall GeomObject's
211  // position at that gime level. I believe this only makes sense
212  // if the t-th history value in the positional timestepper
213  // actually represents previous values (rather than some
214  // generalised quantity). Hence if this function is called with
215  // t>nprev_values(), we issue a warning and terminate the execution.
216  // It *might* be possible that the code still works correctly
217  // even if this condition is violated (e.g. if the GeomObject's
218  // position() function returns the appropriate "generalised"
219  // position value that is required by the timestepping scheme but it's
220  // probably worth flagging this up and forcing the user to manually switch
221  // off this warning if he/she is 100% sure that this is kosher.
222  if (t>node_pt->position_time_stepper_pt()->nprev_values())
223  {
224  std::string error_message =
225  "Trying to update the nodal position at a time level";
226  error_message += "beyond the number of previous values in the nodes'";
227  error_message += "position timestepper. This seems highly suspect!";
228  error_message += "If you're sure the code behaves correctly";
229  error_message += "in your application, remove this warning ";
230  error_message += "or recompile with PARNOID switched off.";
231 
232  std::string function_name =
233  "AlgebraicFSIDrivenCavityMesh::";
234  function_name += "algebraic_node_update()";
235 
236  throw OomphLibError(error_message,
237  OOMPH_CURRENT_FUNCTION,
238  OOMPH_EXCEPTION_LOCATION);
239  }
240 #endif
241 
242  // Extract references for update by copy construction
243  Vector<double> ref_value(node_pt->vector_ref_value());
244 
245  // First reference value: Original x-position. Used as the start point
246  // for the lines connecting the nodes in the vertical direction
247  double x_bottom=ref_value[0];
248 
249  // Second reference value: Fractional position along
250  // straight line from the bottom (at the original x position)
251  // to the reference point on the upper wall
252  double fract=ref_value[1];
253 
254  // Third reference value: Reference local coordinate
255  // in GeomObject that represents the upper wall (local coordinate
256  // in finite element if the wall GeomObject is a finite element mesh)
257  Vector<double> s(1);
258  s[0]=ref_value[2];
259 
260  // Fourth reference value: zeta coordinate on the upper wall
261  // If the wall is a simple GeomObject, zeta[0]=s[0]
262  // but if it's a compound GeomObject (e.g. a finite element mesh)
263  // zeta scales during mesh refinement, whereas s[0] and the
264  // pointer to the geom object have to be re-computed.
265  // double zeta=ref_value[3]; // not needed here
266 
267  // Extract geometric objects for update by copy construction
268  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
269 
270  // Pointer to actual wall geom object (either the same as the wall object
271  // or the pointer to the actual finite element)
272  GeomObject* geom_obj_pt=geom_object_pt[0];
273 
274  // Get position vector to wall at previous timestep t!
275  Vector<double> r_wall(2);
276  geom_obj_pt->position(t,s,r_wall);
277 
278  // Assign new nodal coordinate
279  node_pt->x(t,0)=x_bottom+fract*(r_wall[0]-x_bottom);
280  node_pt->x(t,1)=fract*r_wall[1];
281 
282 }
283 
284 
285 
286 
287 
288 //=====start_setup=================================================
289 /// Setup algebraic mesh update -- assumes that mesh has
290 /// initially been set up with a flush upper wall.
291 //=================================================================
292 template<class ELEMENT>
294 {
295 
296 
297  // Loop over all nodes in mesh
298  unsigned nnod=this->nnode();
299  for (unsigned j=0;j<nnod;j++)
300  {
301  // Get pointer to node -- recall that that Mesh::node_pt(...) has been
302  // overloaded in the AlgebraicMesh class to return a pointer to
303  // an AlgebraicNode.
304  AlgebraicNode* nod_pt=node_pt(j);
305 
306  // Get coordinates
307  double x=nod_pt->x(0);
308  double y=nod_pt->x(1);
309 
310  // Get zeta coordinate on the undeformed wall
311  Vector<double> zeta(1);
312  zeta[0]=x;
313 
314  // Get pointer to geometric (sub-)object and Lagrangian coordinate
315  // on that sub-object. For a wall that is represented by
316  // a single geom object, this simply returns the input.
317  // If the geom object consists of sub-objects (e.g.
318  // if it is a finite element mesh representing a wall,
319  // then we'll obtain the pointer to the finite element
320  // (in its incarnation as a GeomObject) and the
321  // local coordinate in that element.
322  GeomObject* geom_obj_pt;
323  Vector<double> s(1);
324  this->Wall_pt->locate_zeta(zeta,geom_obj_pt,s);
325 
326  // Get position vector to wall:
327  Vector<double> r_wall(2);
328  geom_obj_pt->position(s,r_wall);
329 
330  // Sanity check: Confirm that the wall is in its undeformed position
331 #ifdef PARANOID
332  if ((std::fabs(r_wall[0]-x)>1.0e-15)&&(std::fabs(r_wall[1]-y)>1.0e-15))
333  {
334  std::ostringstream error_stream;
335  error_stream
336  << "Wall must be in its undeformed position when\n"
337  << "algebraic node update information is set up!\n "
338  << "x-discrepancy: " << std::fabs(r_wall[0]-x) << std::endl
339  << "y-discrepancy: " << std::fabs(r_wall[1]-y) << std::endl;
340 
341  throw OomphLibError(
342  error_stream.str(),
343  OOMPH_CURRENT_FUNCTION,
344  OOMPH_EXCEPTION_LOCATION);
345  }
346 #endif
347 
348 
349  // One geometric object is involved in update operation
350  Vector<GeomObject*> geom_object_pt(1);
351 
352  // The actual geometric object (If the wall is simple GeomObject
353  // this is the same as Wall_pt; if it's a compound GeomObject
354  // this points to the sub-object)
355  geom_object_pt[0]=geom_obj_pt;
356 
357  // The update function requires four parameters:
358  Vector<double> ref_value(4);
359 
360  // First reference value: Original x-position
361  ref_value[0]=r_wall[0];
362 
363  // Second reference value: fractional position along
364  // straight line from the bottom (at the original x position)
365  // to the point on the wall)
366  ref_value[1]=y/r_wall[1];
367 
368  // Third reference value: Reference local coordinate
369  // in wall element (local coordinate in FE if we're dealing
370  // with a wall mesh)
371  ref_value[2]=s[0];
372 
373  // Fourth reference value: zeta coordinate on wall
374  // If the wall is a simple GeomObject, zeta[0]=s[0]
375  // but if it's a compound GeomObject (e.g. a finite element mesh)
376  // zeta scales during mesh refinement, whereas s[0] and the
377  // pointer to the geom object have to be re-computed.
378  ref_value[3]=zeta[0];
379 
380  // Setup algebraic update for node: Pass update information
381  nod_pt->add_node_update_info(
382  this, // mesh
383  geom_object_pt, // vector of geom objects
384  ref_value); // vector of ref. values
385  }
386 
387 
388 } //end of setup_algebraic_node_update
389 
390 
391 
392 
393 ////////////////////////////////////////////////////////////////////
394 ////////////////////////////////////////////////////////////////////
395 ////////////////////////////////////////////////////////////////////
396 
397 
398 
399 
400 
401 
402 //========start_update_node_update=================================
403 /// Update the geometric references that are used
404 /// to update node after mesh adaptation.
405 //=================================================================
406 template<class ELEMENT>
408  AlgebraicNode*& node_pt)
409 {
410  // Extract reference values for node update by copy construction
411  Vector<double> ref_value(node_pt->vector_ref_value());
412 
413  // First reference value: Original x-position
414  // double x_bottom=ref_value[0]; // not needed here
415 
416  // Second reference value: fractional position along
417  // straight line from the bottom (at the original x position)
418  // to the point on the wall)
419  // double fract=ref_value[1]; // not needed here
420 
421  // Third reference value: Reference local coordinate
422  // in GeomObject (local coordinate in finite element if the wall
423  // GeomObject is a finite element mesh)
424  // Vector<double> s(1);
425  // s[0]=ref_value[2]; // This needs to be re-computed!
426 
427  // Fourth reference value: intrinsic coordinate on the (possibly
428  // compound) wall.
429  double zeta=ref_value[3];
430 
431  // Extract geometric objects for update by copy construction
432  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
433 
434  // Pointer to actual wall geom object (either the same as wall object
435  // or the pointer to the actual finite element)
436  //GeomObject* geom_obj_pt=geom_object_pt[0]; // This needs to be re-computed!
437 
438  // Get zeta coordinate on wall (as vector)
439  Vector<double> zeta_wall(1);
440  zeta_wall[0]=zeta;
441 
442  // Get pointer to geometric (sub-)object and Lagrangian coordinate
443  // on that sub-object. For a wall that is represented by
444  // a single geom object, this simply returns the input.
445  // If the geom object consists of sub-objects (e.g.
446  // if it is a finite element mesh representing a wall,
447  // then we'll obtain the pointer to the finite element
448  // (in its incarnation as a GeomObject) and the
449  // local coordinate in that element.
450  Vector<double> s(1);
451  GeomObject* geom_obj_pt;
452  this->Wall_pt->locate_zeta(zeta_wall,geom_obj_pt,s);
453 
454  // Update the pointer to the (sub-)GeomObject within which the
455  // reference point is located. (If the wall is simple GeomObject
456  // this is the same as Wall_pt; if it's a compound GeomObject
457  // this points to the sub-object)
458  geom_object_pt[0]=geom_obj_pt;
459 
460  // First reference value: Original x-position
461  // ref_value[0]=r_wall[0]; // unchanged
462 
463  // Second reference value: fractional position along
464  // straight line from the bottom (at the original x position)
465  // to the point on the wall)
466  // ref_value[1]=y/r_wall[1]; // unchanged
467 
468  // Update third reference value: Reference local coordinate
469  // in wall element (local coordinate in FE if we're dealing
470  // with a wall mesh)
471  ref_value[2]=s[0];
472 
473  // Fourth reference value: zeta coordinate on wall
474  // If the wall is a simple GeomObject, zeta[0]=s[0]
475  // but if it's a compound GeomObject (e.g. a finite element mesh)
476  // zeta scales during mesh refinement, whereas s[0] and the
477  // pointer to the geom object have to be re-computed.
478  // ref_value[3]=zeta[0]; //unchanged
479 
480 
481  // Kill the existing node update info
482  node_pt->kill_node_update_info();
483 
484  // Setup algebraic update for node: Pass update information
485  node_pt->add_node_update_info(
486  this, // mesh
487  geom_object_pt, // vector of geom objects
488  ref_value); // vector of ref. values
489 
490 }
491 
492 
493 
494 }
495 #endif
void update_node_update(AlgebraicNode *&node_pt)
Update the node update data for specified node following any mesh adapation.
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:246
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
Definition: mesh.h:201
void kill_node_update_info(const int &id=0)
Erase algebraic node update information for id-th node update function. Id defaults to 0...
cstr elem_len * i
Definition: cfortran.h:607
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
char t
Definition: cfortran.h:572
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
double Gap_fraction
Fraction of the gap next to moving lid, relative to the height of the domain.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void algebraic_node_update(const unsigned &t, AlgebraicNode *&node_pt)
Update nodal position at time level t (t=0: present; t>0: previous)
void add_node_update_info(const int &id, AlgebraicMesh *mesh_pt, const Vector< GeomObject * > &geom_object_pt, const Vector< double > &ref_value, const bool &called_from_constructor=false)
Add algebraic update information for node: What's the ID of the mesh update function (typically used ...
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
e
Definition: cfortran.h:575
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2139
Vector< double > & vector_ref_value()
Return vector of reference values involved in default (usually first) update function.
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
static char t char * s
Definition: cfortran.h:572
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
void setup_boundary_element_info()
Definition: quad_mesh.h:86
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void setup_algebraic_node_update()
Function to setup the algebraic node update.
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:456
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...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806
virtual void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
A geometric object may be composed of may sub-objects (e.g. a finite-element representation of a boun...
Definition: geom_objects.h:352
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
Vector< GeomObject * > & vector_geom_object_pt(const int &id)
Return vector of geometric objects involved in id-th update function.
void remove_boundary_nodes()
Clear all pointers to boundary nodes.
Definition: mesh.cc:208