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