polar_fluid_traction_elements.h
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// Header file for elements that are used to integrate fluid tractions
27// This includes the guts (i.e. equations) because we want to inline them
28// for faster operation, although it slows down the compilation!
29
30#ifndef OOMPH_POLAR_FLUID_TRACTION_ELEMENTS_HEADER
31#define OOMPH_POLAR_FLUID_TRACTION_ELEMENTS_HEADER
32
33// Config header generated by autoconfig
34#ifdef HAVE_CONFIG_H
35#include <config.h>
36#endif
37
38
39// OOMPH-LIB headers
40#include "../generic/Qelements.h"
41
42namespace oomph
43{
44 //======================================================================
45 /// A class for elements that allow the imposition of an applied traction
46 /// to the Navier--Stokes equations
47 /// The geometrical information can be read from the FaceGeometery<ELEMENT>
48 /// class and and thus, we can be generic enough without the need to have
49 /// a separate equations class
50 //======================================================================
51 template<class ELEMENT>
52 class PolarNavierStokesTractionElement : public virtual FaceGeometry<ELEMENT>,
53 public virtual FaceElement
54 {
55 private:
56 /// Pointer to an imposed traction function
57 void (*Traction_fct_pt)(const double& time,
58 const Vector<double>& x,
59 Vector<double>& result);
60
61 /// The highest dimension of the problem
62 unsigned Dim;
63
64 protected:
65 /// Access function that returns the local equation numbers
66 /// for velocity components.
67 /// u_local_eqn(n,i) = local equation number or < 0 if pinned.
68 /// The default is to asssume that n is the local node number
69 /// and the i-th velocity component is the i-th unknown stored at the node.
70 virtual inline int u_local_eqn(const unsigned& n, const unsigned& i)
71 {
72 return nodal_local_eqn(n, i);
73 }
74
75 /// Function to compute the shape and test functions and to return
76 /// the Jacobian of mapping
77 inline double shape_and_test_at_knot(const unsigned& ipt,
78 Shape& psi,
79 Shape& test) const
80 {
81 // Find number of nodes
82 unsigned n_node = nnode();
83 // Calculate the shape functions
84 shape_at_knot(ipt, psi);
85 // Set the test functions to be the same as the shape functions
86 for (unsigned i = 0; i < n_node; i++)
87 {
88 test[i] = psi[i];
89 }
90 // Return the value of the jacobian
91 return J_eulerian_at_knot(ipt);
92 }
93
94
95 /// Function to calculate the traction applied to the fluid
96 void get_traction(double time,
97 const Vector<double>& x,
98 Vector<double>& result)
99 {
100 // If the function pointer is zero return zero
101 if (Traction_fct_pt == 0)
102 {
103 // Loop over dimensions and set body forces to zero
104 for (unsigned i = 0; i < Dim; i++)
105 {
106 result[i] = 0.0;
107 }
108 }
109 // Otherwise call the function
110 else
111 {
112 (*Traction_fct_pt)(time, x, result);
113 }
114 }
115
116 /// This function returns the residuals for the
117 /// traction function.
118 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
120 DenseMatrix<double>& jacobian,
121 DenseMatrix<double>& mass_matrix,
122 unsigned flag);
123 /// Pointer to the angle alpha
124 double* Alpha_pt;
125
126 /// Pointer to the Data item that stores the external pressure
128
129 /// The Data that contains the traded pressure is stored
130 /// as external Data for the element. Which external Data item is it?
132
133 // Traction elements need to know whether they're at the inlet or outlet
134 // as the unit outward normal has a differing sign dependent on
135 // the boundary
136 // -1=inlet, 1=outlet
138
139 // Pointer to homotopy parameter
140 double Eta;
141
142 public:
143 /// Alpha
144 const double& alpha() const
145 {
146 return *Alpha_pt;
147 }
148
149 /// Pointer to Alpha
150 double*& alpha_pt()
151 {
152 return Alpha_pt;
153 }
154
155 /// Function for setting up external pressure
157 {
158 // Set external pressure pointer
159 Pext_data_pt = pext_data_pt;
160
161 // Add to the element's external data so it gets included
162 // in the black-box local equation numbering scheme
164 }
165
166 /// Boundary
167 const int boundary() const
168 {
169 return Boundary;
170 }
171
172 /// Function to set boundary
173 void set_boundary(int bound)
174 {
175 Boundary = bound;
176 }
177
178 /// Eta
179 const double get_eta() const
180 {
181 return Eta;
182 }
183
184 /// Function to set Eta
185 void set_eta(double eta)
186 {
187 Eta = eta;
188 }
189
190 /// Constructor, which takes a "bulk" element and the value of the index
191 /// and its limit
193 const int& face_index)
194 : FaceGeometry<ELEMENT>(), FaceElement()
195 {
196 // Attach the geometrical information to the element. N.B. This function
197 // also assigns nbulk_value from the required_nvalue of the bulk element
198 element_pt->build_face_element(face_index, this);
199
200#ifdef PARANOID
201 {
202 // Check that the element is not a refineable 3d element
203 ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(element_pt);
204 // If it's three-d
205 if (elem_pt->dim() == 3)
206 {
207 // Is it refineable
208 RefineableElement* ref_el_pt =
209 dynamic_cast<RefineableElement*>(elem_pt);
210 if (ref_el_pt != 0)
211 {
212 if (this->has_hanging_nodes())
213 {
214 throw OomphLibError("This flux element will not work correctly "
215 "if nodes are hanging\n",
216 OOMPH_CURRENT_FUNCTION,
217 OOMPH_EXCEPTION_LOCATION);
218 }
219 }
220 }
221 }
222#endif
223
224 // Set the body force function pointer to zero
225 Traction_fct_pt = 0;
226
227 // Set the external pressure pointer to be zero
228 Pext_data_pt = 0;
229
230 // Set the dimension from the dimension of the first node
231 Dim = this->node_pt(0)->ndim();
232
233 // Set Eta to one by default
234 Eta = 1.0;
235 }
236
237 /// Destructor should not delete anything
239
240 // Access function for the imposed traction pointer
241 void (*&traction_fct_pt())(const double& t,
242 const Vector<double>& x,
243 Vector<double>& result)
244 {
245 return Traction_fct_pt;
246 }
247
248 /// This function returns just the residuals
250 {
251 // Call the generic residuals function with flag set to 0
252 // using a dummy matrix argument
256 0);
257 }
258
259 /// This function returns the residuals and the jacobian
261 DenseMatrix<double>& jacobian)
262 {
263 // Call the generic routine with the flag set to 1
265 residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
266 }
267
268 /// Compute the element's residual Vector and the jacobian matrix
269 /// Plus the mass matrix especially for eigenvalue problems
271 Vector<double>& residuals,
272 DenseMatrix<double>& jacobian,
273 DenseMatrix<double>& mass_matrix)
274 {
275 // Call the generic routine with the flag set to 2
277 residuals, jacobian, GeneralisedElement::Dummy_matrix, 2);
278 }
279
280 /// Overload the output function
281 void output(std::ostream& outfile)
282 {
283 FiniteElement::output(outfile);
284 }
285
286 /// Output function: x,y,[z],u,v,[w],p in tecplot format
287 void output(std::ostream& outfile, const unsigned& nplot)
288 {
289 FiniteElement::output(outfile, nplot);
290 }
291
292 /// local velocities
293 double u(const unsigned& l, const unsigned& i)
294 {
295 return nodal_value(l, i);
296 }
297
298 /// local position
299 double x(const unsigned& l, const unsigned& i)
300 {
301 return nodal_position(l, i);
302 }
303 };
304
305
306 /// ////////////////////////////////////////////////////////////////////
307 /// ////////////////////////////////////////////////////////////////////
308 /// ////////////////////////////////////////////////////////////////////
309
310
311 //============================================================================
312 /// Function that returns the residuals for the imposed traction Navier_Stokes
313 /// equations
314 //============================================================================
315 template<class ELEMENT>
318 DenseMatrix<double>& jacobian,
319 DenseMatrix<double>& mass_matrix,
320 unsigned flag)
321 {
322 // Find out how many nodes there are
323 unsigned n_node = nnode();
324
325 // Get continuous time from timestepper of first node
326 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
327
328 // Set up memory for the shape and test functions
329 Shape psif(n_node), testf(n_node);
330
331 // Set the value of n_intpt
332 unsigned n_intpt = integral_pt()->nweight();
333
334 // Get Alpha
335 const double Alpha = alpha();
336
337 // Storage for external pressure
338 double pext = 0.0;
339
340 // Get boundary multiplier
341 // This is necessary because the sign of the traction is
342 // dependent on the boundary
343 const int multiplier = boundary();
344
345 // Get the homotopy parameter (if necessary)
346 const double eta = get_eta();
347
348 // Integers to store local equation numbers
349 int local_eqn = 0, local_unknown = 0, pext_local_eqn = 0,
350 pext_local_unknown = 0;
351
352 /// /////////////////////////////////////NEW//////////////////////////////////////////
353
354 // Get local equation number of external pressure
355 // Note that if we have not passed an external pressure pointer to this
356 // element (and at the same time added it to the element's external data)
357 // than this will be -1 to indicate that it is not a degree of freedom here
358 if (Pext_data_pt == 0)
359 {
360 pext_local_eqn = -1;
361 }
362 else
363 {
364 // If at a non-zero degree of freedom add in the entry
365 pext_local_eqn = external_local_eqn(External_data_number_of_Pext, 0);
366
367 // Get external pressure
368 pext = Pext_data_pt->value(0);
369 }
370
371 // The local unkown number of pext will be the same
372 pext_local_unknown = pext_local_eqn;
373
374 /// /////////////////////////////////////NEW//////////////////////////////////////////
375
376 // Loop over the integration points
377 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
378 {
379 // Get the integral weight
380 double w = integral_pt()->weight(ipt);
381
382 // Find the shape and test functions and return the Jacobian
383 // of the mapping
384 double J = shape_and_test_at_knot(ipt, psif, testf);
385
386 // Premultiply the weights and the Jacobian
387 double W = w * J;
388
389 // Need to find position to feed into Traction function
390 Vector<double> interpolated_x(Dim);
391 Vector<double> interpolated_u(Dim);
392
393 // Initialise to zero
394 for (unsigned i = 0; i < Dim; i++)
395 {
396 interpolated_x[i] = 0.0;
397 interpolated_u[i] = 0.0;
398 }
399
400 // Calculate velocities and derivatives:
401 // Loop over nodes
402 for (unsigned l = 0; l < n_node; l++)
403 {
404 // Loop over directions
405 for (unsigned i = 0; i < Dim; i++)
406 {
407 // Get the nodal value
408 interpolated_u[i] += this->nodal_value(l, i) * psif[l];
409 interpolated_x[i] += this->nodal_position(l, i) * psif[l];
410 }
411 }
412
413 // Get the user-defined traction terms
414 Vector<double> traction(Dim);
415 get_traction(time, interpolated_x, traction);
416
417 // Now add to the appropriate equations
418
419 // Loop over the test functions
420 for (unsigned l = 0; l < n_node; l++)
421 {
422 // Only alter u velocity component
423 {
424 unsigned i = 0;
425 local_eqn = u_local_eqn(l, i);
426 /*IF it's not a boundary condition*/
427 if (local_eqn >= 0)
428 {
429 // Add the user-defined traction terms
430 residuals[local_eqn] -= multiplier * eta * 3.0 *
431 (interpolated_u[i] / interpolated_x[0]) *
432 testf[l] * interpolated_x[0] * Alpha * W;
433
434 /// /////////////////////////////////////NEW//////////////////////////////////////////
435
436 // Plus additional external pressure contribution at inlet
437 // This is zero if we haven't passed a Pext_data_pt to the element
438 residuals[local_eqn] +=
439 pext * testf[l] * interpolated_x[0] * Alpha * W;
440
441 /// /////////////////////////////////////NEW//////////////////////////////////////////
442
443 // CALCULATE THE JACOBIAN
444 if (flag)
445 {
446 // Loop over the velocity shape functions again
447 for (unsigned l2 = 0; l2 < n_node; l2++)
448 {
449 // We only have an i2=0 contribution
450 unsigned i2 = 0;
451 {
452 // If at a non-zero degree of freedom add in the entry
453 local_unknown = u_local_eqn(l2, i2);
454 if (local_unknown >= 0)
455 {
456 // Add contribution to Elemental Matrix
457 jacobian(local_eqn, local_unknown) -=
458 multiplier * eta * 3.0 * (psif[l2] / interpolated_x[0]) *
459 testf[l] * interpolated_x[0] * Alpha * W;
460
461 } // End of (Jacobian's) if not boundary condition statement
462 } // End of i2=0 section
463 } // End of l2 loop
464
465 /// /////////////////////////////////////NEW//////////////////////////////////////////
466 // Add pext's contribution to these residuals
467 // This only needs to be done once hence why it is outside the l2
468 // loop
469 if (pext_local_unknown >= 0)
470 {
471 // Add contribution to Elemental Matrix
472 jacobian(local_eqn, pext_local_unknown) +=
473 testf[l] * interpolated_x[0] * Alpha * W;
474 }
475 /// /////////////////////////////////////NEW//////////////////////////////////////////
476
477 } /*End of Jacobian calculation*/
478
479 } // end of if not boundary condition statement
480 } // End of i=0 section
481
482 } // End of loop over shape functions
483
484 /// /////////////////////////////////////NEW//////////////////////////////////////////
485
486 /// The additional residual for the mass flux
487 /// (ie. the extra equation for pext)
488 /// This is an integral equation along the whole boundary
489 /// It lies outside the loop over shape functions above
490 {
491 /*IF it's not a boundary condition*/
492 if (pext_local_eqn >= 0)
493 {
494 // Add the user-defined traction terms
495 residuals[pext_local_eqn] +=
496 interpolated_u[0] * interpolated_x[0] * Alpha * W;
497
498 // No longer necessary due to my FluxCosntraint element
499 // Now take off a fraction of the desired mass flux
500 // Divided by number of elements and number of int points in each
501 // HACK
502 // residuals[pext_local_eqn] -= (1.0/(30.*3.));
503 // HACK
504
505 // CALCULATE THE JACOBIAN
506 if (flag)
507 {
508 // Loop over the velocity shape functions again
509 for (unsigned l2 = 0; l2 < n_node; l2++)
510 {
511 // We only have an i2=0 contribution
512 unsigned i2 = 0;
513 {
514 // If at a non-zero degree of freedom add in the entry
515 local_unknown = u_local_eqn(l2, i2);
516 if (local_unknown >= 0)
517 {
518 // Add contribution to Elemental Matrix
519 jacobian(pext_local_eqn, local_unknown) +=
520 psif[l2] * interpolated_x[0] * Alpha * W;
521
522 } // End of (Jacobian's) if not boundary condition statement
523 } // End of i2=0 section
524
525 } // End of l2 loop
526 } /*End of Jacobian calculation*/
527
528 } // end of if not boundary condition statement
529 } // End of additional residual for the mass flux
530
531 /// /////////////////////////////////////NEW//////////////////////////////////////////
532
533 } // End of loop over integration points
534
535 } // End of fill_in_generic_residual_contribution
536
537} // End of namespace oomph
538
539#endif
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:293
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5328
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
Definition: elements.h:3050
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement....
Definition: elements.cc:5132
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3220
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2470
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:227
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) external data object to the element and return its index (i....
Definition: elements.cc:307
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
An OomphLibError object which should be thrown when an run-time error is encountered....
A class for elements that allow the imposition of an applied traction to the Navier–Stokes equations ...
void output(std::ostream &outfile)
Overload the output function.
void get_traction(double time, const Vector< double > &x, Vector< double > &result)
Function to calculate the traction applied to the fluid.
double * Alpha_pt
Pointer to the angle alpha.
~PolarNavierStokesTractionElement()
Destructor should not delete anything.
PolarNavierStokesTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the jacobian.
void fill_in_generic_residual_contribution(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
This function returns the residuals for the traction function. flag=1(or 0): do (or don't) compute th...
unsigned External_data_number_of_Pext
The Data that contains the traded pressure is stored as external Data for the element....
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to an imposed traction function.
void(*&)(const double &t, const Vector< double > &x, Vector< double > &result) traction_fct_pt()
unsigned Dim
The highest dimension of the problem.
Data * Pext_data_pt
Pointer to the Data item that stores the external pressure.
double x(const unsigned &l, const unsigned &i)
local position
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping.
double u(const unsigned &l, const unsigned &i)
local velocities
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Compute the element's residual Vector and the jacobian matrix Plus the mass matrix especially for eig...
void set_external_pressure_data(Data *pext_data_pt)
Function for setting up external pressure.
void set_eta(double eta)
Function to set Eta.
void set_boundary(int bound)
Function to set boundary.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
virtual int u_local_eqn(const unsigned &n, const unsigned &i)
Access function that returns the local equation numbers for velocity components. u_local_eqn(n,...
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:60
//////////////////////////////////////////////////////////////////// ////////////////////////////////...