gen_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 GeneralisedAdvection Diffusion elements
32 
33 namespace oomph
34 {
35 
36 ///2D GeneralisedAdvection 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  const unsigned n_node = nnode();
61 
62  //Get the nodal index at which the unknown is stored
63  const unsigned u_nodal_index = u_index_cons_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  const 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  const double peclet = pe();
77 
78  //Get the Peclet*Strouhal number
79  const 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_cons_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_cons_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_cons_adv_diff(ipt,interpolated_x,source);
145 
146 
147  //Get wind
148  //--------
149  Vector<double> wind(DIM);
150  get_wind_cons_adv_diff(ipt,s,interpolated_x,wind);
151 
152  //Get the conserved wind (non-divergence free)
153  Vector<double> conserved_wind(DIM);
154  get_conserved_wind_cons_adv_diff(ipt,s,interpolated_x,conserved_wind);
155 
156 
157  //Get diffusivity tensor
158  DenseMatrix<double> D(DIM,DIM);
159  get_diff_cons_adv_diff(ipt,s,interpolated_x,D);
160 
161  // Assemble residuals and Jacobian
162  //--------------------------------
163 
164  // Loop over the test functions
165  for(unsigned l=0;l<n_node;l++)
166  {
167  //Set the local equation number
168  local_eqn = nodal_local_eqn(l,u_nodal_index);
169 
170  /*IF it's not a boundary condition*/
171  if(local_eqn >= 0)
172  {
173  // Add body force/source term and time derivative
174  residuals[local_eqn] -= (peclet_st*dudt + source)*test(l)*W;
175 
176  // The Generalised Advection Diffusion bit itself
177  for(unsigned k=0;k<DIM;k++)
178  {
179  //Terms that multiply the test function
180  //divergence-free wind
181  double tmp = peclet*wind[k];
182  //If the mesh is moving need to subtract the mesh velocity
183  if(!ALE_is_disabled) {tmp -= peclet_st*mesh_velocity[k];}
184  tmp *= interpolated_dudx[k];
185 
186  //Terms that multiply the derivative of the test function
187  //Advective term
188  double tmp2 = -conserved_wind[k]*interpolated_u;
189  //Now the diuffusive term
190  for(unsigned j=0;j<DIM;j++)
191  {
192  tmp2 += D(k,j)*interpolated_dudx[j];
193  }
194  //Now construct the contribution to the residuals
195  residuals[local_eqn] -= (tmp*test(l) + tmp2*dtestdx(l,k))*W;
196  }
197 
198  // Calculate the jacobian
199  //-----------------------
200  if(flag)
201  {
202  //Loop over the velocity shape functions again
203  for(unsigned l2=0;l2<n_node;l2++)
204  {
205  //Set the number of the unknown
206  local_unknown = nodal_local_eqn(l2,u_nodal_index);
207 
208  //If at a non-zero degree of freedom add in the entry
209  if(local_unknown >= 0)
210  {
211  //Mass matrix term
212  jacobian(local_eqn,local_unknown)
213  -= peclet_st*test(l)*psi(l2)*
214  node_pt(l2)->time_stepper_pt()->weight(1,0)*W;
215 
216  //Add the mass matrix term
217  if(flag==2)
218  {
219  mass_matrix(local_eqn,local_unknown)
220  += peclet_st*test(l)*psi(l2)*W;
221  }
222 
223  //Add contribution to Elemental Matrix
224  for(unsigned k=0;k<DIM;k++)
225  {
226  //Temporary term used in assembly
227  double tmp = -peclet*wind[k];
228  if(!ALE_is_disabled)
229  {tmp -= peclet_st*mesh_velocity[k];}
230  tmp *= dpsidx(l2,k);
231 
232  double tmp2 = -conserved_wind[k]*psi(l2);
233  //Now the diffusive term
234  for(unsigned j=0;j<DIM;j++)
235  {
236  tmp2 += D(k,j)*dpsidx(l2,j);
237  }
238 
239  //Now assemble Jacobian term
240  jacobian(local_eqn,local_unknown)
241  -= (tmp*test(l) + tmp2*dtestdx(l,k))*W;
242  }
243  }
244  }
245  }
246  }
247  }
248 
249 
250  } // End of loop over integration points
251 }
252 
253 
254 
255 //======================================================================
256 /// Self-test: Return 0 for OK
257 //======================================================================
258 template <unsigned DIM>
260 {
261 
262  bool passed=true;
263 
264  // Check lower-level stuff
265  if (FiniteElement::self_test()!=0)
266  {
267  passed=false;
268  }
269 
270  // Return verdict
271  if (passed)
272  {
273  return 0;
274  }
275  else
276  {
277  return 1;
278  }
279 
280 }
281 
282 
283 
284 //======================================================================
285 /// \short Output function:
286 ///
287 /// x,y,u,w_x,w_y or x,y,z,u,w_x,w_y,w_z
288 ///
289 /// nplot points in each coordinate direction
290 //======================================================================
291 template <unsigned DIM>
293  const unsigned &nplot)
294 {
295  //Vector of local coordinates
296  Vector<double> s(DIM);
297 
298 
299  // Tecplot header info
300  outfile << tecplot_zone_string(nplot);
301 
302  const unsigned n_node = this->nnode();
303  Shape psi(n_node);
304  DShape dpsidx(n_node,DIM);
305 
306  // Loop over plot points
307  unsigned num_plot_points=nplot_points(nplot);
308  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
309  {
310  // Get local coordinates of plot point
311  get_s_plot(iplot,nplot,s);
312 
313  // Get Eulerian coordinate of plot point
314  Vector<double> x(DIM);
315  interpolated_x(s,x);
316 
317  for(unsigned i=0;i<DIM;i++)
318  {
319  outfile << x[i] << " ";
320  }
321  outfile << interpolated_u_cons_adv_diff(s) << " ";
322 
323  //Get the gradients
324  (void)this->dshape_eulerian(s,psi,dpsidx);
325  Vector<double> interpolated_dudx(DIM,0.0);
326  for(unsigned n=0;n<n_node;n++)
327  {
328  const double u_ = this->nodal_value(n,0);
329  for(unsigned i=0;i<DIM;i++)
330  {
331  interpolated_dudx[i] += u_*dpsidx(n,i);
332  }
333  }
334 
335  for(unsigned i=0;i<DIM;i++)
336  {
337  outfile << interpolated_dudx[i] << " ";
338  }
339 
340  // Get the wind
341  Vector<double> wind(DIM);
342  // Dummy integration point variable needed
343  unsigned ipt=0;
344  get_wind_cons_adv_diff(ipt,s,x,wind);
345  for(unsigned i=0;i<DIM;i++)
346  {
347  outfile << wind[i] << " ";
348  }
349  outfile << std::endl;
350 
351  }
352 
353  // Write tecplot footer (e.g. FE connectivity lists)
354  write_tecplot_zone_footer(outfile,nplot);
355 
356 }
357 
358 
359 //======================================================================
360 /// C-style output function:
361 ///
362 /// x,y,u or x,y,z,u
363 ///
364 /// nplot points in each coordinate direction
365 //======================================================================
366 template <unsigned DIM>
368  const unsigned &nplot)
369 {
370  //Vector of local coordinates
371  Vector<double> s(DIM);
372 
373  // Tecplot header info
374  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
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  for(unsigned i=0;i<DIM;i++)
385  {
386  fprintf(file_pt,"%g ",interpolated_x(s,i));
387 
388  }
389  fprintf(file_pt,"%g \n",interpolated_u_cons_adv_diff(s));
390  }
391 
392  // Write tecplot footer (e.g. FE connectivity lists)
393  write_tecplot_zone_footer(file_pt,nplot);
394 
395 }
396 
397 
398 
399 //======================================================================
400  /// \short Output exact solution
401  ///
402  /// Solution is provided via function pointer.
403  /// Plot at a given number of plot points.
404  ///
405  /// x,y,u_exact or x,y,z,u_exact
406 //======================================================================
407 template <unsigned DIM>
409  const unsigned &nplot,
411  {
412 
413  //Vector of local coordinates
414  Vector<double> s(DIM);
415 
416  // Vector for coordintes
417  Vector<double> x(DIM);
418 
419  // Tecplot header info
420  outfile << tecplot_zone_string(nplot);
421 
422  // Exact solution Vector (here a scalar)
423  Vector<double> exact_soln(1);
424 
425  // Loop over plot points
426  unsigned num_plot_points=nplot_points(nplot);
427  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
428  {
429 
430  // Get local coordinates of plot point
431  get_s_plot(iplot,nplot,s);
432 
433  // Get x position as Vector
434  interpolated_x(s,x);
435 
436  // Get exact solution at this point
437  (*exact_soln_pt)(x,exact_soln);
438 
439  //Output x,y,...,u_exact
440  for(unsigned i=0;i<DIM;i++)
441  {
442  outfile << x[i] << " ";
443  }
444  outfile << exact_soln[0] << std::endl;
445  }
446 
447  // Write tecplot footer (e.g. FE connectivity lists)
448  write_tecplot_zone_footer(outfile,nplot);
449 
450  }
451 
452 
453 
454 
455 //======================================================================
456  /// \short Validate against exact solution
457  ///
458  /// Solution is provided via function pointer.
459  /// Plot error at a given number of plot points.
460  ///
461 //======================================================================
462 template <unsigned DIM>
465  double& error, double& norm)
466 {
467 
468  // Initialise
469  error=0.0;
470  norm=0.0;
471 
472  //Vector of local coordinates
473  Vector<double> s(DIM);
474 
475  // Vector for coordintes
476  Vector<double> x(DIM);
477 
478  //Find out how many nodes there are in the element
479  unsigned n_node = nnode();
480 
481  Shape psi(n_node);
482 
483  //Set the value of n_intpt
484  unsigned n_intpt = integral_pt()->nweight();
485 
486  // Tecplot header info
487  outfile << "ZONE" << std::endl;
488 
489  // Exact solution Vector (here a scalar)
490  Vector<double> exact_soln(1);
491 
492  //Loop over the integration points
493  for(unsigned ipt=0;ipt<n_intpt;ipt++)
494  {
495 
496  //Assign values of s
497  for(unsigned i=0;i<DIM;i++)
498  {
499  s[i] = integral_pt()->knot(ipt,i);
500  }
501 
502  //Get the integral weight
503  double w = integral_pt()->weight(ipt);
504 
505  // Get jacobian of mapping
506  double J=J_eulerian(s);
507 
508  //Premultiply the weights and the Jacobian
509  double W = w*J;
510 
511  // Get x position as Vector
512  interpolated_x(s,x);
513 
514  // Get FE function value
515  double u_fe=interpolated_u_cons_adv_diff(s);
516 
517  // Get exact solution at this point
518  (*exact_soln_pt)(x,exact_soln);
519 
520  //Output x,y,...,error
521  for(unsigned i=0;i<DIM;i++)
522  {
523  outfile << x[i] << " ";
524  }
525  outfile << exact_soln[0] << " " << exact_soln[0]-u_fe << std::endl;
526 
527  // Add to error and norm
528  norm+=exact_soln[0]*exact_soln[0]*W;
529  error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*W;
530 
531  }
532 
533 }
534 
535 //======================================================================
536 /// \short Calculate the integrated value of the unknown over the element
537 ///
538 //======================================================================
539 template <unsigned DIM>
541 {
542 
543  // Initialise
544  double sum = 0.0;
545 
546  //Vector of local coordinates
547  Vector<double> s(DIM);
548 
549  //Find out how many nodes there are in the element
550  const unsigned n_node = nnode();
551 
552  //Find the index at which the concentration is stored
553  const unsigned u_nodal_index = this->u_index_cons_adv_diff();
554 
555  //Allocate memory for the shape functions
556  Shape psi(n_node);
557 
558  //Set the value of n_intpt
559  const unsigned n_intpt = integral_pt()->nweight();
560 
561  //Loop over the integration points
562  for(unsigned ipt=0;ipt<n_intpt;ipt++)
563  {
564  //Get the integral weight
565  const double w = integral_pt()->weight(ipt);
566 
567  //Get the shape functions
568  this->shape_at_knot(ipt,psi);
569 
570  //Calculate the concentration
571  double interpolated_u = 0.0;
572  for(unsigned l=0;l<n_node;l++)
573  {interpolated_u += this->nodal_value(l,u_nodal_index)*psi(l);}
574 
575  // Get jacobian of mapping
576  const double J=J_eulerian_at_knot(ipt);
577 
578  //Add the values to the sum
579  sum += interpolated_u*w*J;
580  }
581 
582  //return the sum
583  return sum;
584 }
585 
586 
587 //======================================================================
588 // Set the data for the number of Variables at each node
589 //======================================================================
590 template<unsigned DIM, unsigned NNODE_1D>
592 
593 //====================================================================
594 // Force build of templates
595 //====================================================================
599 
603 
607 
611 
612 
613 }
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
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.
cstr elem_len * i
Definition: cfortran.h:607
virtual void fill_in_generic_residual_contribution_cons_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 self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
Definition: elements.cc:4247
double integrate_u()
Integrate the concentration over the element.
A class for all elements that solve the Advection Diffusion equations in conservative form using isop...
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
QGeneralisedAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Advection Diffus...
void output(std::ostream &outfile)
Output with default number of plot points.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
static double Default_peclet_number
Static default value for the Peclet number.