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 Advection Diffusion elements
32 
33 namespace oomph
34 {
35 
36 ///2D Advection Diffusion elements
37 
38 
39 /// Default value for Peclet number
40 template<unsigned DIM>
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 //======================================================================
51 template <unsigned DIM>
54  DenseMatrix<double> &jacobian,
56  &mass_matrix,
57  unsigned flag)
58 {
59  //Find out how many nodes there are
60  unsigned n_node = nnode();
61 
62  //Get the nodal index at which the unknown is stored
63  unsigned u_nodal_index = u_index_adv_diff();
64 
65  //Set up memory for the shape and test functions
66  Shape psi(n_node), test(n_node);
67  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
68 
69  //Set the value of n_intpt
70  unsigned n_intpt = integral_pt()->nweight();
71 
72  //Set the Vector to hold local coordinates
73  Vector<double> s(DIM);
74 
75  //Get Peclet number
76  double peclet = pe();
77 
78  //Get the Peclet*Strouhal number
79  double peclet_st = pe_st();
80 
81  //Integers used to store the local equation number and local unknown
82  //indices for the residuals and jacobians
83  int local_eqn=0, local_unknown=0;
84 
85  //Loop over the integration points
86  for(unsigned ipt=0;ipt<n_intpt;ipt++)
87  {
88 
89  //Assign values of s
90  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
91 
92  //Get the integral weight
93  double w = integral_pt()->weight(ipt);
94 
95  //Call the derivatives of the shape and test functions
96  double J =
97  dshape_and_dtest_eulerian_at_knot_adv_diff(ipt,psi,dpsidx,test,dtestdx);
98 
99  //Premultiply the weights and the Jacobian
100  double W = w*J;
101 
102  //Calculate local values of the solution and its derivatives
103  //Allocate
104  double interpolated_u=0.0;
105  double dudt=0.0;
106  Vector<double> interpolated_x(DIM,0.0);
107  Vector<double> interpolated_dudx(DIM,0.0);
108  Vector<double> mesh_velocity(DIM,0.0);
109 
110 
111  //Calculate function value and derivatives:
112  //-----------------------------------------
113  // Loop over nodes
114  for(unsigned l=0;l<n_node;l++)
115  {
116  //Get the value at the node
117  double u_value = raw_nodal_value(l,u_nodal_index);
118  interpolated_u += u_value*psi(l);
119  dudt += du_dt_adv_diff(l)*psi(l);
120  // Loop over directions
121  for(unsigned j=0;j<DIM;j++)
122  {
123  interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
124  interpolated_dudx[j] += u_value*dpsidx(l,j);
125  }
126  }
127 
128  // Mesh velocity?
129  if (!ALE_is_disabled)
130  {
131  for(unsigned l=0;l<n_node;l++)
132  {
133  for(unsigned j=0;j<DIM;j++)
134  {
135  mesh_velocity[j] += raw_dnodal_position_dt(l,j)*psi(l);
136  }
137  }
138  }
139 
140 
141  //Get source function
142  //-------------------
143  double source;
144  get_source_adv_diff(ipt,interpolated_x,source);
145 
146 
147  //Get wind
148  //--------
149  Vector<double> wind(DIM);
150  get_wind_adv_diff(ipt,s,interpolated_x,wind);
151 
152  // Assemble residuals and Jacobian
153  //--------------------------------
154 
155  // Loop over the test functions
156  for(unsigned l=0;l<n_node;l++)
157  {
158  //Set the local equation number
159  local_eqn = nodal_local_eqn(l,u_nodal_index);
160 
161  /*IF it's not a boundary condition*/
162  if(local_eqn >= 0)
163  {
164  // Add body force/source term and time derivative
165  residuals[local_eqn] -= (peclet_st*dudt + source)*test(l)*W;
166 
167  // The Advection Diffusion bit itself
168  for(unsigned k=0;k<DIM;k++)
169  {
170  //Terms that multiply the test function
171  double tmp = peclet*wind[k];
172  //If the mesh is moving need to subtract the mesh velocity
173  if(!ALE_is_disabled) {tmp -= peclet_st*mesh_velocity[k];}
174  //Now construct the contribution to the residuals
175  residuals[local_eqn] -=
176  interpolated_dudx[k]*(tmp*test(l) + dtestdx(l,k))*W;
177  }
178 
179  // Calculate the jacobian
180  //-----------------------
181  if(flag)
182  {
183  //Loop over the velocity shape functions again
184  for(unsigned l2=0;l2<n_node;l2++)
185  {
186  //Set the number of the unknown
187  local_unknown = nodal_local_eqn(l2,u_nodal_index);
188 
189  //If at a non-zero degree of freedom add in the entry
190  if(local_unknown >= 0)
191  {
192  //Mass matrix term
193  jacobian(local_eqn,local_unknown)
194  -= peclet_st*test(l)*psi(l2)*
195  node_pt(l2)->time_stepper_pt()->weight(1,0)*W;
196 
197  //Add the mass matrix term
198  if(flag==2)
199  {
200  mass_matrix(local_eqn,local_unknown)
201  += peclet_st*test(l)*psi(l2)*W;
202  }
203 
204  //Add contribution to Elemental Matrix
205  for(unsigned i=0;i<DIM;i++)
206  {
207  //Temporary term used in assembly
208  double tmp = peclet*wind[i];
209  if(!ALE_is_disabled) tmp -= peclet_st*mesh_velocity[i];
210  //Now assemble Jacobian term
211  jacobian(local_eqn,local_unknown)
212  -= dpsidx(l2,i)*(tmp*test(l)+dtestdx(l,i))*W;
213  }
214  }
215  }
216  }
217  }
218  }
219 
220 
221  } // End of loop over integration points
222 }
223 
224 
225 
226 
227 //======================================================================
228 /// Self-test: Return 0 for OK
229 //======================================================================
230 template <unsigned DIM>
232 {
233 
234  bool passed=true;
235 
236  // Check lower-level stuff
237  if (FiniteElement::self_test()!=0)
238  {
239  passed=false;
240  }
241 
242  // Return verdict
243  if (passed)
244  {
245  return 0;
246  }
247  else
248  {
249  return 1;
250  }
251 
252 }
253 
254 
255 
256 //======================================================================
257 /// \short Output function:
258 ///
259 /// x,y,u,w_x,w_y or x,y,z,u,w_x,w_y,w_z
260 ///
261 /// nplot points in each coordinate direction
262 //======================================================================
263 template <unsigned DIM>
264 void AdvectionDiffusionEquations<DIM>::output(std::ostream &outfile,
265  const unsigned &nplot)
266 {
267  //Vector of local coordinates
268  Vector<double> s(DIM);
269 
270 
271  // Tecplot header info
272  outfile << tecplot_zone_string(nplot);
273 
274  // Loop over plot points
275  unsigned num_plot_points=nplot_points(nplot);
276  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
277  {
278  // Get local coordinates of plot point
279  get_s_plot(iplot,nplot,s);
280 
281  // Get Eulerian coordinate of plot point
282  Vector<double> x(DIM);
283  interpolated_x(s,x);
284 
285  for(unsigned i=0;i<DIM;i++)
286  {
287  outfile << x[i] << " ";
288  }
289  outfile << interpolated_u_adv_diff(s) << " ";
290 
291  // Get the wind
292  Vector<double> wind(DIM);
293  // Dummy ipt argument needed... ?
294  unsigned ipt=0;
295  get_wind_adv_diff(ipt,s,x,wind);
296  for(unsigned i=0;i<DIM;i++)
297  {
298  outfile << wind[i] << " ";
299  }
300  outfile << std::endl;
301 
302  }
303 
304  // Write tecplot footer (e.g. FE connectivity lists)
305  write_tecplot_zone_footer(outfile,nplot);
306 
307 }
308 
309 
310 //======================================================================
311 /// C-style output function:
312 ///
313 /// x,y,u or x,y,z,u
314 ///
315 /// nplot points in each coordinate direction
316 //======================================================================
317 template <unsigned DIM>
319  const unsigned &nplot)
320 {
321  //Vector of local coordinates
322  Vector<double> s(DIM);
323 
324  // Tecplot header info
325  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
326 
327  // Loop over plot points
328  unsigned num_plot_points=nplot_points(nplot);
329  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
330  {
331 
332  // Get local coordinates of plot point
333  get_s_plot(iplot,nplot,s);
334 
335  for(unsigned i=0;i<DIM;i++)
336  {
337  fprintf(file_pt,"%g ",interpolated_x(s,i));
338 
339  }
340  fprintf(file_pt,"%g \n",interpolated_u_adv_diff(s));
341  }
342 
343  // Write tecplot footer (e.g. FE connectivity lists)
344  write_tecplot_zone_footer(file_pt,nplot);
345 
346 }
347 
348 
349 
350 //======================================================================
351  /// \short Output exact solution
352  ///
353  /// Solution is provided via function pointer.
354  /// Plot at a given number of plot points.
355  ///
356  /// x,y,u_exact or x,y,z,u_exact
357 //======================================================================
358 template <unsigned DIM>
360  const unsigned &nplot,
362  {
363 
364  //Vector of local coordinates
365  Vector<double> s(DIM);
366 
367  // Vector for coordintes
368  Vector<double> x(DIM);
369 
370  // Tecplot header info
371  outfile << tecplot_zone_string(nplot);
372 
373  // Exact solution Vector (here a scalar)
374  Vector<double> exact_soln(1);
375 
376  // Loop over plot points
377  unsigned num_plot_points=nplot_points(nplot);
378  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
379  {
380 
381  // Get local coordinates of plot point
382  get_s_plot(iplot,nplot,s);
383 
384  // Get x position as Vector
385  interpolated_x(s,x);
386 
387  // Get exact solution at this point
388  (*exact_soln_pt)(x,exact_soln);
389 
390  //Output x,y,...,u_exact
391  for(unsigned i=0;i<DIM;i++)
392  {
393  outfile << x[i] << " ";
394  }
395  outfile << exact_soln[0] << std::endl;
396  }
397 
398  // Write tecplot footer (e.g. FE connectivity lists)
399  write_tecplot_zone_footer(outfile,nplot);
400 
401  }
402 
403 
404 
405 
406 //======================================================================
407  /// \short Validate against exact solution
408  ///
409  /// Solution is provided via function pointer.
410  /// Plot error at a given number of plot points.
411  ///
412 //======================================================================
413 template <unsigned DIM>
416  double& error, double& norm)
417 {
418 
419  // Initialise
420  error=0.0;
421  norm=0.0;
422 
423  //Vector of local coordinates
424  Vector<double> s(DIM);
425 
426  // Vector for coordintes
427  Vector<double> x(DIM);
428 
429  //Find out how many nodes there are in the element
430  unsigned n_node = nnode();
431 
432  Shape psi(n_node);
433 
434  //Set the value of n_intpt
435  unsigned n_intpt = integral_pt()->nweight();
436 
437  // Tecplot header info
438  outfile << "ZONE" << std::endl;
439 
440  // Exact solution Vector (here a scalar)
441  Vector<double> exact_soln(1);
442 
443  //Loop over the integration points
444  for(unsigned ipt=0;ipt<n_intpt;ipt++)
445  {
446 
447  //Assign values of s
448  for(unsigned i=0;i<DIM;i++)
449  {
450  s[i] = integral_pt()->knot(ipt,i);
451  }
452 
453  //Get the integral weight
454  double w = integral_pt()->weight(ipt);
455 
456  // Get jacobian of mapping
457  double J=J_eulerian(s);
458 
459  //Premultiply the weights and the Jacobian
460  double W = w*J;
461 
462  // Get x position as Vector
463  interpolated_x(s,x);
464 
465  // Get FE function value
466  double u_fe=interpolated_u_adv_diff(s);
467 
468  // Get exact solution at this point
469  (*exact_soln_pt)(x,exact_soln);
470 
471  //Output x,y,...,error
472  for(unsigned i=0;i<DIM;i++)
473  {
474  outfile << x[i] << " ";
475  }
476  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
477 
478  // Add to error and norm
479  norm+=exact_soln[0]*exact_soln[0]*W;
480  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
481 
482  }
483 
484 }
485 
486 
487 //======================================================================
488 // Set the data for the number of Variables at each node
489 //======================================================================
490 template<unsigned DIM, unsigned NNODE_1D>
492 
493 //====================================================================
494 // Force build of templates
495 //====================================================================
496 template class QAdvectionDiffusionElement<1,2>;
497 template class QAdvectionDiffusionElement<1,3>;
498 template class QAdvectionDiffusionElement<1,4>;
499 
500 template class QAdvectionDiffusionElement<2,2>;
501 template class QAdvectionDiffusionElement<2,3>;
502 template class QAdvectionDiffusionElement<2,4>;
503 
504 template class QAdvectionDiffusionElement<3,2>;
505 template class QAdvectionDiffusionElement<3,3>;
506 template class QAdvectionDiffusionElement<3,4>;
507 
508 
509 }
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...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
cstr elem_len * i
Definition: cfortran.h:607
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
QAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Advection Diffusion element...
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points.
A class for all elements that solve the Advection Diffusion equations using isoparametric elements...
const double & pe() const
Peclet number.
const double & pe_st() const
Peclet number multiplied by Strouhal number.
static char t char * s
Definition: cfortran.h:572
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
void output(std::ostream &outfile)
Output with default number of plot points.
unsigned self_test()
Self-test: Return 0 for OK.