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