fsi.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#include "fsi.h"
27#include "integral.h"
28#include "algebraic_elements.h"
29
30namespace oomph
31{
32 //======================================================================
33 /// Namespace for "global" FSI functions
34 //======================================================================
35 namespace FSI_functions
36 {
37 /// Strouhal number St = a/(UT) for application of no slip condition.
38 /// Initialised to 1.0.
40
41 /// Apply no-slip condition for N.St. on a moving wall node
42 /// u = St dR/dt, where the Strouhal number St = a/(UT) is defined by
43 /// FSI_functions::Strouhal_for_no_slip and is initialised to 1.0.
44 /// Note: This requires the x,y,[z] velocity components to be stored
45 /// in nodal values 0,1,[2]. This is the default for all currently
46 /// existing Navier-Stokes elements. If you use any others,
47 /// use this function at your own risk.
49 {
50 // Find out spatial dimension of node
51 unsigned ndim = node_pt->ndim();
52 // Assign velocity
53 for (unsigned i = 0; i < ndim; i++)
54 {
55 node_pt->set_value(i,
56 Strouhal_for_no_slip * (node_pt->dposition_dt(i)));
57 }
58 }
59
60 } // namespace FSI_functions
61
62
63 //====================================================================
64 /// Flag that allows the suppression of warning messages
65 //====================================================================
67
68 //====================================================================
69 /// Function to describe the local dofs of the element. The ostream
70 /// specifies the output stream to which the description
71 /// is written; the string stores the currently
72 /// assembled output that is ultimately written to the
73 /// output stream by Data::describe_dofs(...); it is typically
74 /// built up incrementally as we descend through the
75 /// call hierarchy of this function when called from
76 /// Problem::describe_dofs(...)
77 //====================================================================
79 std::ostream& out, const std::string& current_string) const
80 {
81 // Call the standard finite element classification.
82 FiniteElement::describe_local_dofs(out, current_string);
83 describe_solid_local_dofs(out, current_string);
84 // Find the number of external field data
85 const unsigned n_external_field_data = nexternal_interaction_field_data();
86 // Now loop over the field data again to assign local equation numbers
87 for (unsigned i = 0; i < n_external_field_data; i++)
88 {
89 std::stringstream conversion;
90 conversion << " of External Interaction Field Data " << i
91 << current_string;
92 std::string in(conversion.str());
94 }
95
96 // Find the number of external geometric data
97 unsigned n_external_geom_data = nexternal_interaction_geometric_data();
98
99 // Now loop over the field data again assign local equation numbers
100 for (unsigned i = 0; i < n_external_geom_data; i++)
101 {
102 std::stringstream conversion;
103 conversion << " of External Interaction Geometric Data " << i
104 << current_string;
105 std::string in(conversion.str());
107 }
108 }
109
110 //================================================================
111 /// Compete build of element: Assign storage -- pass the Eulerian
112 /// dimension of the "adjacent" fluid elements and the
113 /// number of local coordinates required to parametrise
114 /// the wall element. E.g. for a FSIKirchhoffLoveBeam,
115 /// bounding a 2D fluid domain ndim_fluid=2 and nlagr_solid=1.
116 /// Note that, by default, we're only providing storage for fluid
117 /// loading on the "front" of the element. Call
118 /// FSIWallElement::enable_fluid_loading_on_both_sides()
119 /// to enable fluid loading on the back, too.
120 //====================================================================
121 void FSIWallElement::setup_fsi_wall_element(const unsigned& nlagr_solid,
122 const unsigned& ndim_fluid)
123 {
124 // Set storage for underlying GeomObject
125 set_nlagrangian_and_ndim(nlagr_solid, ndim_fluid);
126
127 // Set source element storage - one interaction
128 this->set_ninteraction(1);
129 }
130
131
132 //==================================================================
133 /// Allow element to be loaded by fluid on both
134 /// sides. (Resizes containers for lookup schemes and initialises
135 /// data associated with elements at the "back" of the FSIWallElement
136 /// to NULL.
137 //==================================================================
139 {
140 // Both faces are loaded
142
143 // Set source element storage - two interactions
144 this->set_ninteraction(2);
145 }
146
147
148 //==================================================================
149 /// Get the contribution to the load vector provided by
150 /// the adjacent fluid element: Pass number of integration point
151 /// in solid element,
152 /// and the unit normal vector (pointing into the fluid on the "front"
153 /// of the FSIWallElement) and return the load vector.
154 /// Note that the load is non-dimensionalised on the wall-stress scale,
155 /// i.e. it is obtained by computing the traction (on the fluid stress-scale)
156 /// from the adjacent fluid element and then multiplying it by
157 /// the stress-scale-ratio \f$ Q. \f$.
158 /// If element is loaded on both sides, the fluid load computed
159 /// from the reverse element is \b subtracted, because it's
160 /// computed with the same normal vector.
161 //==================================================================
162 void FSIWallElement::fluid_load_vector(const unsigned& intpt,
163 const Vector<double>& N,
164 Vector<double>& load)
165 {
166 // Size of load vector
167 unsigned n_load = load.size();
168
169 // Initialise
170 for (unsigned i = 0; i < n_load; i++) load[i] = 0.0;
171
172 // Create vector for fluid load
173 Vector<double> fluid_load(n_load);
174
175 // Loop over front and back if required: Get number of fluid-loaded faces
176 unsigned n_loaded_face = 2;
177 if (Only_front_is_loaded_by_fluid) n_loaded_face = 1;
178
179 for (unsigned face = 0; face < n_loaded_face; face++)
180 {
181 // Get the local coordinate in the fluid element (copy
182 // operation for Vector)
183 Vector<double> s_adjacent(external_element_local_coord(face, intpt));
184
185 // Call the load function for adjacent element
186 FSIFluidElement* el_f_pt =
187 dynamic_cast<FSIFluidElement*>(external_element_pt(face, intpt));
188 if (el_f_pt != 0)
189 {
190 el_f_pt->get_load(s_adjacent, N, fluid_load);
191 }
192 else
193 {
194#ifdef PARANOID
196 {
197 std::ostringstream warning_stream;
198 warning_stream
199 << "Info: No adjacent element set in FSIWallElement.\n\n"
200 << "Note: you can disable this message by setting \n "
201 << "FSIWallElement::Dont_warn_about_missing_adjacent_fluid_elements"
202 << "\n to true or recompiling without PARANOID.\n";
203 OomphLibWarning(warning_stream.str(),
204 "FSIWallElement::fluid_load_vector()",
205 OOMPH_EXCEPTION_LOCATION);
206 }
207#endif
208 }
209
210 // Adjust sign of fluid traction. If normal on the front
211 // points into the fluid, the normal at the back points away
212 // from it.
213 double sign = 1.0;
214 if (face == 1) sign = -1.0;
215
216 // The load is scaled by the stress-scale ratio Q
217 for (unsigned i = 0; i < n_load; i++)
218 {
219 load[i] += fluid_load[i] * sign * q();
220 }
221
222 } // end of loop over faces
223 }
224
225
226 //==================================================================
227 /// Update the nodal positions in all fluid elements that affect
228 /// the traction on this FSIWallElement
229 //==================================================================
231 {
232 // Don't update elements repeatedly
233 std::map<FSIFluidElement*, bool> done;
234
235 // Number of integration points
236 unsigned n_intpt = integral_pt()->nweight();
237
238 // Loop over front and back if required: Get number of fluid-loaded faces
239 unsigned n_loaded_face = 2;
240 if (Only_front_is_loaded_by_fluid) n_loaded_face = 1;
241 for (unsigned face = 0; face < n_loaded_face; face++)
242 {
243 // Loop over all integration points in wall element
244 for (unsigned iint = 0; iint < n_intpt; iint++)
245 {
246 // Get fluid element that affects this integration point
247 FSIFluidElement* el_f_pt =
248 dynamic_cast<FSIFluidElement*>(external_element_pt(face, iint));
249
250 // Is there an adjacent fluid element?
251 if (el_f_pt != 0)
252 {
253 // Have we updated its positions yet?
254 if (!done[el_f_pt])
255 {
256 // Update nodal positions
257 el_f_pt->node_update();
258 done[el_f_pt] = true;
259 }
260 }
261 }
262 }
263 }
264
265
266 //=================================================================
267 /// Static default value for the ratio of stress scales
268 /// used in the fluid and solid equations (default is 1.0)
269 //=================================================================
271
272
273 //=======================================================================
274 /// Overload the function that must return all field data involved
275 /// in the interactions from the external (fluid) element. It allows
276 /// the velocity degrees of freedom to be ignored if we want to
277 /// ignore the shear stresses when computing the Jacobian.
278 //=======================================================================
280 Vector<std::set<FiniteElement*>> const& external_elements_pt,
281 std::set<std::pair<Data*, unsigned>>& paired_interaction_data)
282 {
283 // Loop over each inteaction
284 const unsigned n_interaction = this->ninteraction();
285 for (unsigned i = 0; i < n_interaction; i++)
286 {
287 // Loop over each element in the set
288 for (std::set<FiniteElement*>::const_iterator it =
289 external_elements_pt[i].begin();
290 it != external_elements_pt[i].end();
291 it++)
292 {
293 // Cast the element the specific fluid element
294 FSIFluidElement* external_fluid_el_pt =
295 dynamic_cast<FSIFluidElement*>(*it);
296
297 // Only add pressure load if so desired
299 {
300 // Add the "pressure" load data
301 external_fluid_el_pt->identify_pressure_data(paired_interaction_data);
302 }
303 else
304 {
305 // Add the "direct" load data (usually velocities and pressures)
306 // to the set
307 external_fluid_el_pt->identify_load_data(paired_interaction_data);
308 }
309 }
310 } // End of loop over interactions
311 }
312
313 //===================================================================
314 /// Function that must return all geometric data involved
315 /// in the desired interactions from the external element
316 //===================================================================
318 Vector<std::set<FiniteElement*>> const& external_elements_pt,
319 std::set<Data*>& external_geometric_data_pt)
320 {
321 // If we are ignoring the shear stress, then we can ignore all
322 // external geometric data
324 {
325 return;
326 }
327 else
328 {
329 // Loop over each inteaction
330 const unsigned n_interaction = this->ninteraction();
331 for (unsigned i = 0; i < n_interaction; i++)
332 {
333 // Loop over each element in the set
334 for (std::set<FiniteElement*>::const_iterator it =
335 external_elements_pt[i].begin();
336 it != external_elements_pt[i].end();
337 it++)
338 {
339 (*it)->identify_geometric_data(external_geometric_data_pt);
340 }
341 }
342 }
343 }
344
345
346} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
virtual void describe_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the dofs of the Node. The ostream specifies the output stream to which the descr...
Definition: nodes.cc:939
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
unsigned nexternal_interaction_geometric_data() const
Return the number of geometric Data items that affect the external interactions in this element: i....
unsigned nexternal_interaction_field_data() const
Return the number of Data items that affect the external interactions in this element....
Data ** External_interaction_field_data_pt
/ Storage for pointers to external field Data that affect the interactions in the elemenet
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
Data ** External_interaction_geometric_data_pt
/ Storage for pointers to external geometric Data that affect the interactions in the elemenet
unsigned ninteraction() const
Return the number of interactions in the element.
/////////////////////////////////////////////////////////////////////////
Definition: fsi.h:63
virtual void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)=0
Add to the set paired_load_data pairs containing.
virtual void get_load(const Vector< double > &s, const Vector< double > &N, Vector< double > &load)=0
Compute the load vector that is applied by current element (at its local coordinate s) onto the adjac...
virtual void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)=0
Add to the set paired_pressure_data pairs containing.
void enable_fluid_loading_on_both_sides()
Allow element to be loaded by fluid on both sides. (Resizes containers for lookup schemes and initial...
Definition: fsi.cc:138
void identify_all_field_data_for_external_interaction(Vector< std::set< FiniteElement * > > const &external_elements_pt, std::set< std::pair< Data *, unsigned > > &paired_iteraction_data)
Overload the function that must return all field data involved in the interactions from the external ...
Definition: fsi.cc:279
void identify_all_geometric_data_for_external_interaction(Vector< std::set< FiniteElement * > > const &external_elements_pt, std::set< Data * > &external_geometric_data_pt)
Function that must return all geometric data involved in the desired interactions from the external e...
Definition: fsi.cc:317
void setup_fsi_wall_element(const unsigned &nlagr_solid, const unsigned &ndim_fluid)
Setup: Assign storage – pass the Eulerian dimension of the "adjacent" fluid elements and the number o...
Definition: fsi.cc:121
void fluid_load_vector(const unsigned &intpt, const Vector< double > &N, Vector< double > &load)
Get FE Jacobian by systematic finite differencing w.r.t. nodal positition Data, internal and external...
Definition: fsi.cc:162
bool Ignore_shear_stress_in_jacobian
Set this flag to true to ignore shear stress component of load when calculating the Jacobian,...
Definition: fsi.h:547
const double & q() const
Return the ratio of the stress scales used to non-dimensionalise the fluid and solid equations....
Definition: fsi.h:241
static double Default_Q_Value
Static default value for the ratio of stress scales used in the fluid and solid equations (default is...
Definition: fsi.h:528
static bool Dont_warn_about_missing_adjacent_fluid_elements
Static flag that allows the suppression of warning messages.
Definition: fsi.h:208
void node_update_adjacent_fluid_elements()
Update the nodal positions in all fluid elements that affect the traction on this FSIWallElement.
Definition: fsi.cc:230
bool Only_front_is_loaded_by_fluid
Is the element exposed to (and hence loaded by) fluid only on its "front"? True by default....
Definition: fsi.h:535
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
Definition: fsi.cc:78
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element[s]. The ostream specifies the output stream to whi...
Definition: elements.cc:1709
virtual void node_update()
Update the positions of all nodes in the element using each node update function. The default impleme...
Definition: elements.cc:5072
void set_nlagrangian_and_ndim(const unsigned &n_lagrangian, const unsigned &n_dim)
Set # of Lagrangian and Eulerian coordinates.
Definition: geom_objects.h:183
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
double dposition_dt(const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt, either directly or via hanging node representatio...
Definition: nodes.cc:2659
An OomphLibWarning object which should be created as a temporary object to issue a warning....
void describe_solid_local_dofs(std::ostream &out, const std::string &current_string) const
Classifies dofs locally for solid specific aspects.
Definition: elements.cc:6874
double Strouhal_for_no_slip
Strouhal number St = a/(UT) for application of no slip condition. Initialised to 1....
Definition: fsi.cc:39
void apply_no_slip_on_moving_wall(Node *node_pt)
Apply no-slip condition for N.St. on a moving wall node u = St dR/dt, where the Strouhal number St = ...
Definition: fsi.cc:48
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...