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-2022 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
33namespace 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)
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)
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
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
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.
////////////////////////////////////////////////////////////////////
Vector< GeomObject * > & vector_geom_object_pt(const int &id)
Return vector of geometric objects involved in id-th update function.
void kill_node_update_info(const int &id=0)
Erase algebraic node update information for id-th node update function. Id defaults to 0.
Vector< double > & vector_ref_value()
Return vector of reference values involved in default (usually first) update function.
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 ...
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...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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:2218
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
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:381
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:243
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:190
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
void remove_boundary_nodes()
Clear all pointers to boundary nodes.
Definition: mesh.cc:204
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2394
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
void setup_boundary_element_info()
Setup lookup schemes which establish whic elements are located next to mesh's boundaries (wrapper to ...
Definition: quad_mesh.h:73
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.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...