steady_axisym_advection_diffusion_elements.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// 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 //Non-inline functions for Axisymmetric Advection Diffusion elements
32 
33 namespace oomph
34 {
35 
36 /// 2D Steady Axisymmetric Advection Diffusion elements
37 
38 /// Default value for Peclet number
41 
42 
43 //======================================================================
44 /// \short Compute element residual Vector and/or element Jacobian matrix
45 ///
46 /// flag=1: compute both
47 /// flag=0: compute only residual Vector
48 ///
49 /// Pure version without hanging nodes
50 //======================================================================
53  DenseMatrix<double> &jacobian,
54  DenseMatrix<double> &mass_matrix,
55  unsigned flag)
56 {
57  //Find out how many nodes there are
58  unsigned n_node = nnode();
59 
60  //Get the nodal index at which the unknown is stored
61  unsigned u_nodal_index = u_index_axisym_adv_diff();
62 
63  //Set up memory for the shape and test functions
64  Shape psi(n_node), test(n_node);
65  DShape dpsidx(n_node,2), dtestdx(n_node,2);
66 
67  //Set the value of n_intpt
68  unsigned n_intpt = integral_pt()->nweight();
69 
70  //Set the Vector to hold local coordinates
71  Vector<double> s(2);
72 
73  //Get Physical Variables from Element
74  double scaled_peclet = pe();
75 
76  //Integers used to store the local equation number and local unknown
77  //indices for the residuals and jacobians
78  int local_eqn=0, local_unknown=0;
79 
80  //Loop over the integration points
81  for(unsigned ipt=0;ipt<n_intpt;ipt++)
82  {
83 
84  //Assign values of s
85  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
86 
87  //Get the integral weight
88  double w = integral_pt()->weight(ipt);
89 
90  //Call the derivatives of the shape and test functions
91  double J =
92  dshape_and_dtest_eulerian_at_knot_adv_diff(ipt,psi,dpsidx,test,dtestdx);
93 
94  //Premultiply the weights and the Jacobian
95  double W = w*J;
96 
97  //Calculate local values of the solution and its derivatives
98  //Allocate
99  double interpolated_u = 0.0;
101  Vector<double> interpolated_dudx(2,0.0);
102 
103  //Calculate function value and derivatives:
104  //-----------------------------------------
105  //Loop over nodes
106  for(unsigned l=0;l<n_node;l++)
107  {
108  //Get the value at the node
109  double u_value = raw_nodal_value(l,u_nodal_index);
110  interpolated_u += u_value*psi(l);
111  //Loop over the two coordinates directions
112  for(unsigned j=0;j<2;j++)
113  {
114  interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
115  interpolated_dudx[j] += u_value*dpsidx(l,j);
116  }
117  }
118 
119  //Get source function
120  //-------------------
121  double source;
122  get_source_axisym_adv_diff(ipt,interpolated_x,source);
123 
124 
125  //Get wind
126  //--------
127  Vector<double> wind(3);
128  get_wind_axisym_adv_diff(ipt,s,interpolated_x,wind);
129 
130  //r is the first position component
131  double r = interpolated_x[0];
132 
133  //TEMPERATURE EQUATION (Neglected viscous dissipation)
134  //Assemble residuals and Jacobian
135  //--------------------------------
136 
137  // Loop over the test functions
138  for(unsigned l=0;l<n_node;l++)
139  {
140  //Set the local equation number
141  local_eqn = nodal_local_eqn(l,u_nodal_index);
142 
143  /*IF it's not a boundary condition*/
144  if(local_eqn >= 0)
145  {
146  //Add body force/source term
147  residuals[local_eqn] += r*source*test(l)*W;
148 
149  //The Advection Diffusion bit itself
150  for(unsigned k=0;k<2;k++)
151  {
152  //Construct the contribution to the residuals
153  residuals[local_eqn] -=
154  r*interpolated_dudx[k]*(scaled_peclet*wind[k]*test(l) + dtestdx(l,k))*W;
155  }
156 
157  // Calculate the jacobian
158  //-----------------------
159  if(flag)
160  {
161  //Loop over the velocity shape functions again
162  for(unsigned l2=0;l2<n_node;l2++)
163  {
164  //Set the number of the unknown
165  local_unknown = nodal_local_eqn(l2,u_nodal_index);
166 
167  //If at a non-zero degree of freedom add in the entry
168  if(local_unknown >= 0)
169  {
170  //Add contribution to Elemental Matrix
171  for(unsigned i=0;i<2;i++)
172  {
173  //Assemble Jacobian term
174  jacobian(local_eqn,local_unknown) -=
175  r*dpsidx(l2,i)*(scaled_peclet*wind[i]*test(l) + dtestdx(l,i))*W;
176  }
177  }
178  }
179  }
180  }
181  }
182 
183 
184  } // End of loop over integration points
185 }
186 
187 
188 
189 
190 //======================================================================
191 /// Self-test: Return 0 for OK
192 //======================================================================
193 //template <unsigned DIM>
195 {
196 
197  bool passed = true;
198 
199  // Check lower-level stuff
200  if (FiniteElement::self_test()!=0)
201  {
202  passed = false;
203  }
204 
205  // Return verdict
206  if (passed)
207  {
208  return 0;
209  }
210  else
211  {
212  return 1;
213  }
214 
215 }
216 
217 
218 
219 //======================================================================
220 /// \short Output function:
221 ///
222 /// r,z,u,w_r,w_z
223 ///
224 /// nplot points in each coordinate direction
225 //======================================================================
227  const unsigned &nplot)
228 {
229  //Vector of local coordinates
230  Vector<double> s(2);
231 
232  //Tecplot header info
233  outfile << tecplot_zone_string(nplot);
234 
235  //Loop over plot points
236  unsigned num_plot_points = nplot_points(nplot);
237  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
238  {
239  //Get local coordinates of plot point
240  get_s_plot(iplot,nplot,s);
241 
242  //Get Eulerian coordinate of plot point
243  Vector<double> x(2);
244  interpolated_x(s,x);
245 
246  for(unsigned i=0;i<2;i++)
247  {
248  outfile << x[i] << " ";
249  }
250  outfile << interpolated_u_adv_diff(s) << " ";
251 
252  //Get the wind
253  Vector<double> wind(2);
254  //Dummy ipt argument needed... ?
255  unsigned ipt = 0;
256  get_wind_axisym_adv_diff(ipt,s,x,wind);
257  for(unsigned i=0;i<2;i++)
258  {
259  outfile << wind[i] << " ";
260  }
261  outfile << std::endl;
262 
263  }
264 
265  //Write tecplot footer (e.g. FE connectivity lists)
266  write_tecplot_zone_footer(outfile,nplot);
267 
268 }
269 
270 
271 //======================================================================
272 /// C-style output function:
273 ///
274 /// r,z,u
275 ///
276 /// nplot points in each coordinate direction
277 //======================================================================
278 //template <unsigned DIM>
280  const unsigned &nplot)
281 {
282  //Vector of local coordinates
283  Vector<double> s(2);
284 
285  //Tecplot header info
286  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
287 
288  //Loop over plot points
289  unsigned num_plot_points = nplot_points(nplot);
290  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
291  {
292 
293  //Get local coordinates of plot point
294  get_s_plot(iplot,nplot,s);
295 
296  for(unsigned i=0;i<2;i++)
297  {
298  fprintf(file_pt,"%g ",interpolated_x(s,i));
299 
300  }
301  fprintf(file_pt,"%g \n",interpolated_u_adv_diff(s));
302  }
303 
304  //Write tecplot footer (e.g. FE connectivity lists)
305  write_tecplot_zone_footer(file_pt,nplot);
306 
307 }
308 
309 //======================================================================
310  /// \short Output exact solution
311  ///
312  /// Solution is provided via function pointer.
313  /// Plot at a given number of plot points.
314  ///
315  /// r,z,u_exact
316 //======================================================================
317 //template <unsigned DIM>
319  const unsigned &nplot,
321  {
322 
323  //Vector of local coordinates
324  Vector<double> s(2);
325 
326  //Vector for coordintes
327  Vector<double> x(2);
328 
329  //Tecplot header info
330  outfile << tecplot_zone_string(nplot);
331 
332  //Exact solution Vector (here a scalar)
333  Vector<double> exact_soln(1);
334 
335  //Loop over plot points
336  unsigned num_plot_points=nplot_points(nplot);
337  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
338  {
339 
340  //Get local coordinates of plot point
341  get_s_plot(iplot,nplot,s);
342 
343  //Get x position as Vector
344  interpolated_x(s,x);
345 
346  //Get exact solution at this point
347  (*exact_soln_pt)(x,exact_soln);
348 
349  //Output x,y,...,u_exact
350  for(unsigned i=0;i<2;i++)
351  {
352  outfile << x[i] << " ";
353  }
354  outfile << exact_soln[0] << std::endl;
355  }
356 
357  //Write tecplot footer (e.g. FE connectivity lists)
358  write_tecplot_zone_footer(outfile,nplot);
359 
360  }
361 
362 
363 //======================================================================
364  /// \short Validate against exact solution
365  ///
366  /// Solution is provided via function pointer.
367  /// Plot error at a given number of plot points.
368  ///
369 //======================================================================
370 //template <unsigned DIM>
373  double& error,
374  double& norm)
375 {
376 
377  //Initialise
378  error=0.0;
379  norm=0.0;
380 
381  //Vector of local coordinates
382  Vector<double> s(2);
383 
384  //Vector for coordintes
385  Vector<double> x(2);
386 
387  //Find out how many nodes there are in the element
388  unsigned n_node = nnode();
389 
390  Shape psi(n_node);
391 
392  //Set the value of n_intpt
393  unsigned n_intpt = integral_pt()->nweight();
394 
395  //Tecplot header info
396  outfile << "ZONE" << std::endl;
397 
398  //Exact solution Vector (here a scalar)
399  Vector<double> exact_soln(1);
400 
401  //Loop over the integration points
402  for(unsigned ipt=0;ipt<n_intpt;ipt++)
403  {
404 
405  //Assign values of s
406  for(unsigned i=0;i<2;i++)
407  {
408  s[i] = integral_pt()->knot(ipt,i);
409  }
410 
411  //Get the integral weight
412  double w = integral_pt()->weight(ipt);
413 
414  //Get jacobian of mapping
415  double J=J_eulerian(s);
416 
417  //Premultiply the weights and the Jacobian
418  double W = w*J;
419 
420  //Get x position as Vector
421  interpolated_x(s,x);
422 
423  //Get FE function value
424  double u_fe=interpolated_u_adv_diff(s);
425 
426  //Get exact solution at this point
427  (*exact_soln_pt)(x,exact_soln);
428 
429  //Output x,y,...,error
430  for(unsigned i=0;i<2;i++)
431  {
432  outfile << x[i] << " ";
433  }
434  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
435 
436  //Add to error and norm
437  norm+=exact_soln[0]*exact_soln[0]*W;
438  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
439 
440  }
441 
442 }
443 
444 
445 //======================================================================
446 // Set the data for the number of Variables at each node
447 //======================================================================
448 template<unsigned NNODE_1D>
450 Initial_Nvalue = 1;
451 
452 //====================================================================
453 // Force build of templates
454 //====================================================================
455 
459 
460 }
virtual void get_wind_axisym_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...
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
virtual double dshape_and_dtest_eulerian_at_knot_adv_diff(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
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
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,z,u_exact at nplot^2 plot points.
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation. NOTE: Moved to cc file because of a possible compiler bug in gcc (yes, really!). The move to the cc file avoids inlining which appears to cause problems (only) when compiled with gcc and -O3; offensive "illegal read" is in optimised-out section of code and data that is allegedly illegal is readily readable (by other means) just before this function is called so I can't really see how we could possibly be responsible for this...
Definition: elements.cc:1657
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
QSteadyAxisymAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Axisymmetric Ad...
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
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
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:1383
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2446
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
static double Default_peclet_number
Static default value for the Peclet number.
virtual unsigned u_index_axisym_adv_diff() const
Broken assignment operator.
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 void get_source_axisym_adv_diff(const unsigned &ipt, const Vector< double > &x, double &source) const
Get source term at (Eulerian) position x. This function is virtual to allow overloading in multi-phys...
double interpolated_u_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
void output(std::ostream &outfile)
Output with default number of plot points.
virtual void fill_in_generic_residual_contribution_adv_diff(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add the element's contribution to its residual vector only (if flag=and/or element Jacobian matrix...