supg_advection_diffusion_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#ifndef OOMPH_SUPG_ADV_DIFF_ELEMENTS_HEADER
27#define OOMPH_SUPG_ADV_DIFF_ELEMENTS_HEADER
28
29#include "../advection_diffusion/refineable_advection_diffusion_elements.h"
30
31namespace oomph
32{
33 //======================================================================
34 /// QSUPGAdvectionDiffusionElement<DIM,NNODE_1D> elements are
35 /// SUPG-stabilised Advection Diffusion elements with
36 /// NNODE_1D nodal points in each coordinate direction. Inherits
37 /// from QAdvectionDiffusionElement and overwrites their
38 /// test functions
39 ///
40 //======================================================================
41 template<unsigned DIM, unsigned NNODE_1D>
43 : public virtual QAdvectionDiffusionElement<DIM, NNODE_1D>
44 {
45 public:
46 /// Constructor: Call constructors for underlying
47 /// QAdvectionDiffusion element. Initialise stabilisation parameter
48 /// to zero
50 : QAdvectionDiffusionElement<DIM, NNODE_1D>()
51 {
52 Tau_SUPG = 0.0;
53 }
54
55 /// Get stabilisation parameter for the element
56 double get_Tau_SUPG()
57 {
58 return Tau_SUPG;
59 }
60
61
62 /// Set stabilisation parameter for the element to zero
64 {
65 Tau_SUPG = 0.0;
66 }
67
68
69 /// Compute stabilisation parameter for the element
71 {
72 // Find out how many nodes there are
73 unsigned n_node = this->nnode();
74
75 // Set up memory for the shape functions and their derivatives
76 Shape psi(n_node), test(n_node);
77 DShape dpsidx(n_node, DIM);
78
79 // Evaluate everything at the element centroid
80 Vector<double> s(DIM, 0.0);
81
82 // Call the geometrical shape functions and derivatives
84
85 // Calculate Eulerian coordinates
87
88 // Loop over nodes
89 for (unsigned l = 0; l < n_node; l++)
90 {
91 // Loop over directions (we're in 2D)
92 for (unsigned j = 0; j < DIM; j++)
93 {
94 interpolated_x[j] += this->nodal_position(l, j) * psi[l];
95 }
96 }
97
98 // Element size: Choose the max. diagonal
99 double h = 0;
100 if (DIM == 1)
101 {
102 h = std::fabs(this->vertex_node_pt(1)->x(0) -
103 this->vertex_node_pt(0)->x(0));
104 }
105 else if (DIM == 2)
106 {
107 h =
108 pow(this->vertex_node_pt(3)->x(0) - this->vertex_node_pt(0)->x(0),
109 2) +
110 pow(this->vertex_node_pt(3)->x(1) - this->vertex_node_pt(0)->x(1), 2);
111 double h1 =
112 pow(this->vertex_node_pt(2)->x(0) - this->vertex_node_pt(1)->x(0),
113 2) +
114 pow(this->vertex_node_pt(2)->x(1) - this->vertex_node_pt(1)->x(1), 2);
115 if (h1 > h) h = h1;
116 h = sqrt(h);
117 }
118 else if (DIM == 3)
119 {
120 // diagonals are from nodes 0-7, 1-6, 2-5, 3-4
121 unsigned n1 = 0;
122 unsigned n2 = 7;
123 for (unsigned i = 0; i < 4; i++)
124 {
125 double h1 =
126 pow(this->vertex_node_pt(n1)->x(0) - this->vertex_node_pt(n2)->x(0),
127 2) +
128 pow(this->vertex_node_pt(n1)->x(1) - this->vertex_node_pt(n2)->x(1),
129 2) +
130 pow(this->vertex_node_pt(n1)->x(2) - this->vertex_node_pt(n2)->x(2),
131 2);
132 if (h1 > h) h = h1;
133 n1++;
134 n2--;
135 }
136 h = sqrt(h);
137 }
138
139 // Get wind
140 Vector<double> wind(DIM);
141 // Dummy ipt argument?
142 unsigned ipt = 0;
143 this->get_wind_adv_diff(ipt, s, interpolated_x, wind);
144 double abs_wind = 0;
145 for (unsigned j = 0; j < DIM; j++)
146 {
147 abs_wind += wind[j] * wind[j];
148 }
149 abs_wind = sqrt(abs_wind);
150
151 // Mesh Peclet number
152 double Pe_mesh = 0.5 * this->pe() * h * abs_wind;
153
154 // Stabilisation parameter
155 if (Pe_mesh > 1.0)
156 {
157 Tau_SUPG = h / (2.0 * abs_wind) * (1.0 - 1.0 / Pe_mesh);
158 }
159 else
160 {
161 Tau_SUPG = 0.0;
162 }
163 }
164
165
166 /// Output function:
167 /// x,y,u,w_x,w_y,tau_supg or x,y,z,u,w_x,w_y,w_z,tau_supg
168 /// nplot points in each coordinate direction
169 void output(std::ostream& outfile, const unsigned& nplot)
170 {
171 // Vector of local coordinates
172 Vector<double> s(DIM);
173
174 // Tecplot header info
175 outfile << this->tecplot_zone_string(nplot);
176
177 // Loop over plot points
178 unsigned num_plot_points = this->nplot_points(nplot);
179 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
180 {
181 // Get local coordinates of plot point
182 this->get_s_plot(iplot, nplot, s);
183
184 // Get Eulerian coordinate of plot point
185 Vector<double> x(DIM);
186 this->interpolated_x(s, x);
187
188 for (unsigned i = 0; i < DIM; i++)
189 {
190 outfile << x[i] << " ";
191 }
192 outfile << this->interpolated_u_adv_diff(s) << " ";
193
194 // Get the wind
195 Vector<double> wind(DIM);
196 // Dummy ipt argument
197 unsigned ipt = 0;
198 this->get_wind_adv_diff(ipt, s, x, wind);
199 for (unsigned i = 0; i < DIM; i++)
200 {
201 outfile << wind[i] << " ";
202 }
203
204 // Output stabilisation parameter
205 outfile << Tau_SUPG << std::endl;
206 }
207
208 // Write tecplot footer (e.g. FE connectivity lists)
209 this->write_tecplot_zone_footer(outfile, nplot);
210 }
211
212 /// Output at default number of plot points
213 void output(std::ostream& outfile)
214 {
215 FiniteElement::output(outfile);
216 }
217
218 /// C-style output
219 void output(FILE* file_pt)
220 {
221 FiniteElement::output(file_pt);
222 }
223
224 /// C_style output at n_plot points
225 void output(FILE* file_pt, const unsigned& n_plot)
226 {
227 FiniteElement::output(file_pt, n_plot);
228 }
229
230
231 protected:
232 /// Shape, test functions & derivs. w.r.t. to global coords. Return
233 /// Jacobian.
235 Shape& psi,
236 DShape& dpsidx,
237 Shape& test,
238 DShape& dtestdx) const;
239
240
241 /// Shape, test functions & derivs. w.r.t. to global coords. Return
242 /// Jacobian.
243 double dshape_and_dtest_eulerian_at_knot_adv_diff(const unsigned& ipt,
244 Shape& psi,
245 DShape& dpsidx,
246 Shape& test,
247 DShape& dtestdx) const;
248
249 /// SUPG stabilisation parameter
250 double Tau_SUPG;
251 };
252
253
254 /// //////////////////////////////////////////////////////////////////////
255 /// //////////////////////////////////////////////////////////////////////
256 /// //////////////////////////////////////////////////////////////////////
257
258
259 //======================================================================
260 /// Refineable version of QSUPGAdvectionDiffusionElement.
261 /// Inherit from the standard QSUPGAdvectionDiffusionElement and the
262 /// appropriate refineable geometric element and the refineable equations.
263 //======================================================================
264 template<unsigned DIM, unsigned NNODE_1D>
266 : public QSUPGAdvectionDiffusionElement<DIM, NNODE_1D>,
267 public virtual RefineableAdvectionDiffusionEquations<DIM>,
268 public virtual RefineableQElement<DIM>
269 {
270 public:
271 /// Constructor: Pass refinement level to refineable quad element
272 /// (default 0 = root)
276 RefineableQElement<DIM>(),
277 QSUPGAdvectionDiffusionElement<DIM, NNODE_1D>()
278 {
279 }
280
281
282 /// Broken copy constructor
285 delete;
286
287 /// Broken assignment operator
290
291 /// Number of continuously interpolated values: 1
293 {
294 return 1;
295 }
296
297 /// Number of vertex nodes in the element
298 unsigned nvertex_node() const
299 {
301 }
302
303 /// Pointer to the j-th vertex node in the element
304 Node* vertex_node_pt(const unsigned& j) const
305 {
307 }
308
309 /// Rebuild from sons: empty
310 void rebuild_from_sons(Mesh*& mesh_pt) {}
311
312 /// Order of recovery shape functions for Z2 error estimation:
313 /// Same order as shape functions.
315 {
316 return (NNODE_1D - 1);
317 }
318
319 /// Perform additional hanging node procedures for variables
320 /// that are not interpolated by all nodes. Empty.
322 };
323
324
325} // namespace oomph
326
327#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
double interpolated_u_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
virtual void get_wind_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Get wind at (Eulerian) position x and/or local coordinate s. This function is virtual to allow overlo...
const double & pe() const
Peclet number.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
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
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
virtual unsigned nvertex_node() const
Return the number of vertex nodes in this element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2491
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
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
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
virtual Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2500
A general mesh class.
Definition: mesh.h:67
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
QSUPGAdvectionDiffusionElement<DIM,NNODE_1D> elements are SUPG-stabilised Advection Diffusion element...
double dshape_and_dtest_eulerian_adv_diff(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
void compute_stabilisation_parameter()
Compute stabilisation parameter for the element.
void output(std::ostream &outfile)
Output at default number of plot points.
void switch_off_stabilisation()
Set stabilisation parameter for the element to zero.
double get_Tau_SUPG()
Get stabilisation parameter for the element.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,u,w_x,w_y,tau_supg or x,y,z,u,w_x,w_y,w_z,tau_supg nplot points in each coordina...
void output(FILE *file_pt, const unsigned &n_plot)
C_style output at n_plot points.
QSUPGAdvectionDiffusionElement()
Constructor: Call constructors for underlying QAdvectionDiffusion element. Initialise stabilisation p...
double dshape_and_dtest_eulerian_at_knot_adv_diff(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
A version of the Advection Diffusion equations that can be used with non-uniform mesh refinement....
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
////////////////////////////////////////////////////////////////////// //////////////////////////////...
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned nvertex_node() const
Number of vertex nodes in the element.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void operator=(const RefineableQSUPGAdvectionDiffusionElement< DIM, NNODE_1D > &)=delete
Broken assignment operator.
RefineableQSUPGAdvectionDiffusionElement()
Constructor: Pass refinement level to refineable quad element (default 0 = root)
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 1.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
RefineableQSUPGAdvectionDiffusionElement(const RefineableQSUPGAdvectionDiffusionElement< DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...