gen_axisym_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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 2016) $
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 //Header file for Advection Diffusion elements
31 #ifndef OOMPH_GEN_AXISYM_ADV_DIFF_ELEMENTS_HEADER
32 #define OOMPH_GEN_AXISYM_ADV_DIFF_ELEMENTS_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 //OOMPH-LIB headers
41 #include "../generic/nodes.h"
42 #include "../generic/Qelements.h"
43 #include "../generic/oomph_utilities.h"
44 
45 namespace oomph
46 {
47 
48 //=============================================================
49 // DOXYERROR: the maths below is wrong somehow but I have no idea what it's supposed to be doing!
50 /// \short A class for all elements that solve the Advection
51 /// Diffusion equations in conservative form using isoparametric elements
52 /// in a cylindrical polar coordinate system.
53 /// \mbox{\boldmath$\nabla\cdot$} \left(
54 /// Pe \mbox{\boldmath$w$}(\mbox{\boldmath$x$}) u
55 /// - D(\mbox{\boldmath$x$)\mbox{\boldmath$\nabla$} u\right)
56 /// = f(\mbox{\boldmath$x$})
57 /// This contains the generic maths. Shape functions, geometric
58 /// mapping etc. must get implemented in derived class.
59 //=============================================================
60 class GeneralisedAxisymAdvectionDiffusionEquations :
61  public virtual FiniteElement
62 {
63 
64 public:
65 
66  /// \short Function pointer to source function fct(x,f(x)) --
67  /// x is a Vector!
68  typedef void (*GeneralisedAxisymAdvectionDiffusionSourceFctPt)
69  (const Vector<double>& x, double& f);
70 
71  /// \short Function pointer to wind function fct(x,w(x)) --
72  /// x is a Vector!
73  typedef void (*GeneralisedAxisymAdvectionDiffusionWindFctPt)
74  (const Vector<double>& x, Vector<double>& wind);
75 
76 
77  /// \short Function pointer to a diffusivity function
78  typedef void (*GeneralisedAxisymAdvectionDiffusionDiffFctPt)
79  (const Vector<double> &x, DenseMatrix<double> &D);
80 
81  /// \short Constructor: Initialise the Source_fct_pt and Wind_fct_pt
82  /// to null and set (pointer to) Peclet number to default
83  GeneralisedAxisymAdvectionDiffusionEquations() :
86  {
87  //Set Peclet number to default
89  //Set Peclet Strouhal number to default
91  }
92 
93  /// Broken copy constructor
94  GeneralisedAxisymAdvectionDiffusionEquations(
95  const GeneralisedAxisymAdvectionDiffusionEquations& dummy)
96  {
97  BrokenCopy::broken_copy("GeneralisedAxisymAdvectionDiffusionEquations");
98  }
99 
100  /// Broken assignment operator
101 //Commented out broken assignment operator because this can lead to a conflict warning
102 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
103 //realise that two separate implementations of the broken function are the same and so,
104 //quite rightly, it shouts.
105  /*void operator=(const GeneralisedAxisymAdvectionDiffusionEquations&)
106  {
107  BrokenCopy::broken_assign("GeneralisedAxisymAdvectionDiffusionEquations");
108  }*/
109 
110  /// \short Return the index at which the unknown value
111  /// is stored. The default value, 0, is appropriate for single-physics
112  /// problems, when there is only one variable, the value that satisfies
113  /// the advection-diffusion equation.
114  /// In derived multi-physics elements, this function should be overloaded
115  /// to reflect the chosen storage scheme. Note that these equations require
116  /// that the unknown is always stored at the same index at each node.
117  virtual inline unsigned u_index_cons_axisym_adv_diff() const {return 0;}
118 
119 /// \short du/dt at local node n.
120  /// Uses suitably interpolated value for hanging nodes.
121  double du_dt_cons_axisym_adv_diff(const unsigned &n) const
122  {
123  // Get the data's timestepper
124  TimeStepper* time_stepper_pt= this->node_pt(n)->time_stepper_pt();
125 
126  //Initialise dudt
127  double dudt=0.0;
128  //Loop over the timesteps, if there is a non Steady timestepper
129  if (!time_stepper_pt->is_steady())
130  {
131  //Find the index at which the variable is stored
132  const unsigned u_nodal_index = u_index_cons_axisym_adv_diff();
133 
134  // Number of timsteps (past & present)
135  const unsigned n_time = time_stepper_pt->ntstorage();
136 
137  for(unsigned t=0;t<n_time;t++)
138  {
139  dudt += time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
140  }
141  }
142  return dudt;
143  }
144 
145  /// \short Disable ALE, i.e. assert the mesh is not moving -- you do this
146  /// at your own risk!
147  void disable_ALE()
148  {
149  ALE_is_disabled=true;
150  }
151 
152 
153  /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
154  /// when evaluating the time-derivative. Note: By default, ALE is
155  /// enabled, at the expense of possibly creating unnecessary work
156  /// in problems where the mesh is, in fact, stationary.
157  void enable_ALE()
158  {
159  ALE_is_disabled=false;
160  }
161 
162 
163  /// Output with default number of plot points
164  void output(std::ostream &outfile)
165  {
166  unsigned nplot=5;
167  output(outfile,nplot);
168  }
169 
170  /// \short Output FE representation of soln: r,z,u at
171  /// nplot^2 plot points
172  void output(std::ostream &outfile, const unsigned &nplot);
173 
174  /// C_style output with default number of plot points
175  void output(FILE* file_pt)
176  {
177  unsigned n_plot=5;
178  output(file_pt,n_plot);
179  }
180 
181  /// \short C-style output FE representation of soln: r,z,u at
182  /// n_plot^2 plot points
183  void output(FILE* file_pt, const unsigned &n_plot);
184 
185 
186  /// Output exact soln: r,z,u_exact at nplot^2 plot points
187  void output_fct(std::ostream &outfile, const unsigned &nplot,
189  exact_soln_pt);
190 
191  /// \short Output exact soln: r,z,,u_exact at
192  /// nplot^2 plot points (dummy time-dependent version to
193  /// keep intel compiler happy)
194  virtual void output_fct(std::ostream &outfile, const unsigned &nplot,
195  const double& time,
197  {
198  throw OomphLibError(
199  "There is no time-dependent output_fct() for Advection Diffusion elements",
200  OOMPH_CURRENT_FUNCTION,
201  OOMPH_EXCEPTION_LOCATION);
202  }
203 
204 
205  /// Get error against and norm of exact solution
206  void compute_error(std::ostream &outfile,
208  exact_soln_pt, double& error, double& norm);
209 
210 
211  /// Dummy, time dependent error checker
212  void compute_error(std::ostream &outfile,
214  exact_soln_pt,
215  const double& time, double& error, double& norm)
216  {
217  throw OomphLibError(
218  "No time-dependent compute_error() for Advection Diffusion elements",
219  OOMPH_CURRENT_FUNCTION,
220  OOMPH_EXCEPTION_LOCATION);
221  }
222 
223  /// \short Integrate the concentration over the element
224  double integrate_u();
225 
226 
227  /// Access function: Pointer to source function
228  GeneralisedAxisymAdvectionDiffusionSourceFctPt& source_fct_pt()
229  {return Source_fct_pt;}
230 
231 
232  /// Access function: Pointer to source function. Const version
233  GeneralisedAxisymAdvectionDiffusionSourceFctPt source_fct_pt() const
234  {return Source_fct_pt;}
235 
236 
237  /// Access function: Pointer to wind function
238  GeneralisedAxisymAdvectionDiffusionWindFctPt& wind_fct_pt()
239  {return Wind_fct_pt;}
240 
241 
242  /// Access function: Pointer to wind function. Const version
243  GeneralisedAxisymAdvectionDiffusionWindFctPt wind_fct_pt() const
244  {return Wind_fct_pt;}
245 
246 
247  /// Access function: Pointer to additional (conservative) wind function
248  GeneralisedAxisymAdvectionDiffusionWindFctPt& conserved_wind_fct_pt()
249  {return Conserved_wind_fct_pt;}
250 
251 
252  /// \short Access function: Pointer to additional (conservative)
253  /// wind function.
254  /// Const version
255  GeneralisedAxisymAdvectionDiffusionWindFctPt conserved_wind_fct_pt() const
256  {return Conserved_wind_fct_pt;}
257 
258  /// Access function: Pointer to diffusion function
259  GeneralisedAxisymAdvectionDiffusionDiffFctPt& diff_fct_pt()
260  {return Diff_fct_pt;}
261 
262  /// Access function: Pointer to diffusion function. Const version
263  GeneralisedAxisymAdvectionDiffusionDiffFctPt diff_fct_pt() const
264  {return Diff_fct_pt;}
265 
266  /// Peclet number
267  const double &pe() const {return *Pe_pt;}
268 
269  /// Pointer to Peclet number
270  double* &pe_pt() {return Pe_pt;}
271 
272  /// Peclet number multiplied by Strouhal number
273  const double &pe_st() const {return *PeSt_pt;}
274 
275  /// Pointer to Peclet number multipled by Strouha number
276  double* &pe_st_pt() {return PeSt_pt;}
277 
278  /// \short Get source term at (Eulerian) position x. This function is
279  /// virtual to allow overloading in multi-physics problems where
280  /// the strength of the source function might be determined by
281  /// another system of equations
282  inline virtual void get_source_cons_axisym_adv_diff(const unsigned& ipt,
283  const Vector<double>& x,
284  double& source) const
285  {
286  //If no source function has been set, return zero
287  if(Source_fct_pt==0) {source = 0.0;}
288  else
289  {
290  // Get source strength
291  (*Source_fct_pt)(x,source);
292  }
293  }
294 
295  /// \short Get wind at (Eulerian) position x and/or local coordinate s.
296  /// This function is
297  /// virtual to allow overloading in multi-physics problems where
298  /// the wind function might be determined by
299  /// another system of equations
300  inline virtual void get_wind_cons_axisym_adv_diff(const unsigned& ipt,
301  const Vector<double> &s,
302  const Vector<double>& x,
303  Vector<double>& wind) const
304  {
305  //If no wind function has been set, return zero
306  //There are three components of the wind, but only two matter
307  if(Wind_fct_pt==0)
308  {
309  for(unsigned i=0;i<3;i++) {wind[i]= 0.0;}
310  }
311  else
312  {
313  // Get wind
314  (*Wind_fct_pt)(x,wind);
315  }
316  }
317 
318 
319  /// \short Get additional (conservative)
320  /// wind at (Eulerian) position x and/or local coordinate s.
321  /// This function is
322  /// virtual to allow overloading in multi-physics problems where
323  /// the wind function might be determined by
324  /// another system of equations
326  const unsigned& ipt,
327  const Vector<double> &s,
328  const Vector<double>& x,
329  Vector<double>& wind) const
330  {
331  //If no wind function has been set, return zero
332  if(Conserved_wind_fct_pt==0)
333  {
334  for(unsigned i=0;i<3;i++) {wind[i]= 0.0;}
335  }
336  else
337  {
338  // Get wind
339  (*Conserved_wind_fct_pt)(x,wind);
340  }
341  }
342 
343 
344  /// \short Get diffusivity tensor at (Eulerian) position
345  /// x and/or local coordinate s.
346  /// This function is
347  /// virtual to allow overloading in multi-physics problems where
348  /// the wind function might be determined by
349  /// another system of equations
350  inline virtual void get_diff_cons_axisym_adv_diff(const unsigned& ipt,
351  const Vector<double> &s,
352  const Vector<double>& x,
353  DenseMatrix<double>& D) const
354  {
355  //If no wind function has been set, return identity
356  //Again three components, but not all of them matter
357  //Those in the theta direction are ignored
358  if(Diff_fct_pt==0)
359  {
360  for(unsigned i=0;i<3;i++)
361  {
362  for(unsigned j=0;j<3;j++)
363  {
364  if(i==j) {D(i,j) = 1.0;} else {D(i,j) = 0.0;}
365  }
366  }
367  }
368  else
369  {
370  // Get diffusivity tensor
371  (*Diff_fct_pt)(x,D);
372  }
373  }
374 
375 
376  /// Get flux: \f$\mbox{flux}[i] = \mbox{d}u / \mbox{d}x_i \f$
377  void get_flux(const Vector<double>& s, Vector<double>& flux) const
378  {
379  //Find out how many nodes there are in the element
380  const unsigned n_node = this->nnode();
381 
382  //Get the nodal index at which the unknown is stored
383  const unsigned u_nodal_index = this->u_index_cons_axisym_adv_diff();
384 
385  //Set up memory for the shape and test functions
386  Shape psi(n_node);
387  DShape dpsidx(n_node,2);
388 
389  //Call the derivatives of the shape and test functions
390  dshape_eulerian(s,psi,dpsidx);
391 
392  //Initialise to zero
393  for(unsigned j=0;j<2;j++) {flux[j] = 0.0;}
394 
395  // Loop over nodes
396  for(unsigned l=0;l<n_node;l++)
397  {
398  const double u_value = this->nodal_value(l,u_nodal_index);
399  //Add in the derivative directions
400  flux[0] += u_value*dpsidx(l,0);
401  flux[1] += u_value*dpsidx(l,1);
402  }
403  }
404 
405  /// Get flux: \f$\mbox{flux}[i] = \mbox{d}u / \mbox{d}x_i \f$
407  Vector<double>& total_flux) const
408  {
409  //Find out how many nodes there are in the element
410  const unsigned n_node = nnode();
411 
412  //Get the nodal index at which the unknown is stored
413  const unsigned u_nodal_index = u_index_cons_axisym_adv_diff();
414 
415  //Set up memory for the shape and test functions
416  Shape psi(n_node);
417  DShape dpsidx(n_node,2);
418 
419  //Call the derivatives of the shape and test functions
420  dshape_eulerian(s,psi,dpsidx);
421 
422  //Storage for the Eulerian position
423  Vector<double> interpolated_x(2,0.0);
424  //Storage for the concentration
425  double interpolated_u = 0.0;
426  //Storage for the derivatives of the concentration
427  Vector<double> interpolated_dudx(2,0.0);
428 
429  // Loop over nodes
430  for(unsigned l=0;l<n_node;l++)
431  {
432  //Get the value at the node
433  const double u_value = this->nodal_value(l,u_nodal_index);
434  interpolated_u += u_value*psi(l);
435  // Loop over directions
436  for(unsigned j=0;j<2;j++)
437  {
438  interpolated_x[j] += this->nodal_position(l,j)*psi(l);
439  interpolated_dudx[j] += u_value*dpsidx(l,j);
440  }
441  }
442 
443  //Dummy integration point
444  unsigned ipt=0;
445 
446  //Get the conserved wind (non-divergence free)
447  Vector<double> conserved_wind(3);
448  get_conserved_wind_cons_axisym_adv_diff(ipt,s,interpolated_x,conserved_wind);
449 
450  //Get diffusivity tensor
451  DenseMatrix<double> D(3,3);
452  get_diff_cons_axisym_adv_diff(ipt,s,interpolated_x,D);
453 
454  //Calculate the total flux made up of the diffusive flux
455  //and the conserved wind only bother with the
456  //first two components in each case becuase there can be no
457  //variation in the azimuthal direction
458  for(unsigned i=0;i<2;i++)
459  {
460  total_flux[i] = 0.0;
461  for(unsigned j=0;j<2;j++)
462  {
463  total_flux[i] += D(i,j)*interpolated_dudx[j];
464  }
465  total_flux[i] -= conserved_wind[i]*interpolated_u;
466  }
467  }
468 
469 
470  /// Add the element's contribution to its residual vector (wrapper)
472  {
473  //Call the generic residuals function with flag set to 0 and using
474  //a dummy matrix
478  }
479 
480 
481  /// \short Add the element's contribution to its residual vector and
482  /// the element Jacobian matrix (wrapper)
484  DenseMatrix<double> &jacobian)
485  {
486  //Call the generic routine with the flag set to 1
488  residuals,jacobian,GeneralisedElement::Dummy_matrix,1);
489  }
490 
491 
492  /// Add the element's contribution to its residuals vector,
493  /// jacobian matrix and mass matrix
495  Vector<double> &residuals, DenseMatrix<double> &jacobian,
496  DenseMatrix<double> &mass_matrix)
497  {
498  //Call the generic routine with the flag set to 2
500  jacobian,mass_matrix,2);
501  }
502 
503 
504  /// Return FE representation of function value u(s) at local coordinate s
506  const
507  {
508  //Find number of nodes
509  unsigned n_node = nnode();
510 
511  //Get the nodal index at which the unknown is stored
512  unsigned u_nodal_index = u_index_cons_axisym_adv_diff();
513 
514  //Local shape function
515  Shape psi(n_node);
516 
517  //Find values of shape function
518  shape(s,psi);
519 
520  //Initialise value of u
521  double interpolated_u = 0.0;
522 
523  //Loop over the local nodes and sum
524  for(unsigned l=0;l<n_node;l++)
525  {
526  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
527  }
528 
529  return(interpolated_u);
530  }
531 
532 
533  /// \short Self-test: Return 0 for OK
534  unsigned self_test();
535 
536 protected:
537 
538 
539  /// \short Shape/test functions and derivs w.r.t. to global coords at
540  /// local coord. s; return Jacobian of mapping
542  const Vector<double> &s,
543  Shape &psi,
544  DShape &dpsidx,
545  Shape &test,
546  DShape &dtestdx) const=0;
547 
548  /// \short Shape/test functions and derivs w.r.t. to global coords at
549  /// integration point ipt; return Jacobian of mapping
551  const unsigned &ipt,
552  Shape &psi,
553  DShape &dpsidx,
554  Shape &test,
555  DShape &dtestdx)
556  const=0;
557 
558  /// \short Add the element's contribution to its residual vector only
559  /// (if flag=and/or element Jacobian matrix
561  Vector<double> &residuals, DenseMatrix<double> &jacobian,
562  DenseMatrix<double> &mass_matrix, unsigned flag);
563 
564  /// Pointer to global Peclet number
565  double *Pe_pt;
566 
567  /// Pointer to global Peclet number multiplied by Strouhal number
568  double *PeSt_pt;
569 
570  /// Pointer to source function:
571  GeneralisedAxisymAdvectionDiffusionSourceFctPt Source_fct_pt;
572 
573  /// Pointer to wind function:
574  GeneralisedAxisymAdvectionDiffusionWindFctPt Wind_fct_pt;
575 
576  /// Pointer to additional (conservative) wind function:
577  GeneralisedAxisymAdvectionDiffusionWindFctPt Conserved_wind_fct_pt;
578 
579  /// Pointer to diffusivity funciton
580  GeneralisedAxisymAdvectionDiffusionDiffFctPt Diff_fct_pt;
581 
582  /// \short Boolean flag to indicate if ALE formulation is disabled when
583  /// time-derivatives are computed. Only set to false if you're sure
584  /// that the mesh is stationary.
586 
587  private:
588 
589  /// Static default value for the Peclet number
590  static double Default_peclet_number;
591 
592 };
593 
594 
595 ///////////////////////////////////////////////////////////////////////////
596 ///////////////////////////////////////////////////////////////////////////
597 ///////////////////////////////////////////////////////////////////////////
598 
599 
600 
601 //======================================================================
602 /// \short QGeneralisedAxisymAdvectionDiffusionElement elements are
603 /// linear/quadrilateral/brick-shaped Advection Diffusion elements with
604 /// isoparametric interpolation for the function.
605 //======================================================================
606 template <unsigned NNODE_1D>
608  public virtual QElement<2,NNODE_1D>,
609  public virtual GeneralisedAxisymAdvectionDiffusionEquations
610 {
611 
612 private:
613 
614  /// \short Static array of ints to hold number of variables at
615  /// nodes: Initial_Nvalue[n]
616  static const unsigned Initial_Nvalue;
617 
618  public:
619 
620 
621  ///\short Constructor: Call constructors for QElement and
622  /// Advection Diffusion equations
623  QGeneralisedAxisymAdvectionDiffusionElement() : QElement<2,NNODE_1D>(),
624  GeneralisedAxisymAdvectionDiffusionEquations()
625  { }
626 
627  /// Broken copy constructor
630  {
631  BrokenCopy::broken_copy("QGeneralisedAxisymAdvectionDiffusionElement");
632  }
633 
634  /// Broken assignment operator
635  /*void operator=(const QGeneralisedAxisymAdvectionDiffusionElement<NNODE_1D>&)
636  {
637  BrokenCopy::broken_assign("QGeneralisedAxisymAdvectionDiffusionElement");
638  }*/
639 
640  /// \short Required # of `values' (pinned or dofs)
641  /// at node n
642  inline unsigned required_nvalue(const unsigned &n) const
643  {return Initial_Nvalue;}
644 
645  /// \short Output function:
646  /// r,z,u
647  void output(std::ostream &outfile)
649 
650  /// \short Output function:
651  /// r,z,u at n_plot^2 plot points
652  void output(std::ostream &outfile, const unsigned &n_plot)
654 
655 
656  /// \short C-style output function:
657  /// r,z,u
658  void output(FILE* file_pt)
659  {
661  }
662 
663  /// \short C-style output function:
664  /// r,z,u at n_plot^2 plot points
665  void output(FILE* file_pt, const unsigned &n_plot)
666  {
668  }
669 
670  /// \short Output function for an exact solution:
671  /// r,z,u_exact at n_plot^2 plot points
672  void output_fct(std::ostream &outfile, const unsigned &n_plot,
674  exact_soln_pt)
676  ::output_fct(outfile,n_plot,exact_soln_pt);}
677 
678 
679  /// \short Output function for a time-dependent exact solution.
680  /// r,z,u_exact at n_plot^2 plot points
681  /// (Calls the steady version)
682  void output_fct(std::ostream &outfile, const unsigned &n_plot,
683  const double& time,
685  exact_soln_pt)
686  {
688  output_fct(outfile,n_plot,time,exact_soln_pt);
689  }
690 
691 
692 protected:
693 
694  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
696  const Vector<double> &s,
697  Shape &psi,
698  DShape &dpsidx,
699  Shape &test,
700  DShape &dtestdx) const;
701 
702  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
703  /// integration point ipt. Return Jacobian.
705  const unsigned& ipt,
706  Shape &psi,
707  DShape &dpsidx,
708  Shape &test,
709  DShape &dtestdx)
710  const;
711 
712 };
713 
714 //Inline functions:
715 
716 
717 //======================================================================
718 /// \short Define the shape functions and test functions and derivatives
719 /// w.r.t. global coordinates and return Jacobian of mapping.
720 ///
721 /// Galerkin: Test functions = shape functions
722 //======================================================================
723 template<unsigned NNODE_1D>
726  Shape &psi,
727  DShape &dpsidx,
728  Shape &test,
729  DShape &dtestdx) const
730 {
731  //Call the geometrical shape functions and derivatives
732  double J = this->dshape_eulerian(s,psi,dpsidx);
733 
734  //Loop over the test functions and derivatives and set them equal to the
735  //shape functions
736  for(unsigned i=0;i<NNODE_1D;i++)
737  {
738  test[i] = psi[i];
739  for(unsigned j=0;j<2;j++)
740  {
741  dtestdx(i,j) = dpsidx(i,j);
742  }
743  }
744 
745  //Return the jacobian
746  return J;
747 }
748 
749 
750 
751 //======================================================================
752 /// Define the shape functions and test functions and derivatives
753 /// w.r.t. global coordinates and return Jacobian of mapping.
754 ///
755 /// Galerkin: Test functions = shape functions
756 //======================================================================
757 template<unsigned NNODE_1D>
760  const unsigned &ipt,
761  Shape &psi,
762  DShape &dpsidx,
763  Shape &test,
764  DShape &dtestdx) const
765 {
766  //Call the geometrical shape functions and derivatives
767  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
768 
769  //Set the test functions equal to the shape functions (pointer copy)
770  test = psi;
771  dtestdx = dpsidx;
772 
773  //Return the jacobian
774  return J;
775 }
776 
777 
778 ////////////////////////////////////////////////////////////////////////
779 ////////////////////////////////////////////////////////////////////////
780 ////////////////////////////////////////////////////////////////////////
781 
782 
783 
784 //=======================================================================
785 /// \short Face geometry for the QGeneralisedAxisymAdvectionDiffusionElement elements:
786 /// The spatial dimension of the face elements is one lower than that
787 /// of the bulk element but they have the same number of points along
788 /// their 1D edges.
789 //=======================================================================
790 template<unsigned NNODE_1D>
791 class FaceGeometry<QGeneralisedAxisymAdvectionDiffusionElement<NNODE_1D> >:
792  public virtual QElement<1,NNODE_1D>
793 {
794 
795  public:
796 
797  /// \short Constructor: Call the constructor for the
798  /// appropriate lower-dimensional QElement
799  FaceGeometry() : QElement<1,NNODE_1D>() {}
800 
801 };
802 
803 }
804 
805 #endif
virtual void get_wind_cons_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...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void get_total_flux(const Vector< double > &s, Vector< double > &total_flux) const
Get flux: .
GeneralisedAxisymAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,z,u at n_plot^2 plot points.
GeneralisedAxisymAdvectionDiffusionWindFctPt & conserved_wind_fct_pt()
Access function: Pointer to additional (conservative) wind function.
GeneralisedAxisymAdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
GeneralisedAxisymAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
double dshape_and_dtest_eulerian_at_knot_cons_axisym_adv_diff(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
QGeneralisedAxisymAdvectionDiffusionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
cstr elem_len * i
Definition: cfortran.h:607
void output(FILE *file_pt, const unsigned &n_plot)
C-style output FE representation of soln: r,z,u at n_plot^2 plot points.
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.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:557
char t
Definition: cfortran.h:572
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:375
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
double interpolated_u_cons_axisym_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
const double & pe() const
Peclet number.
virtual double dshape_and_dtest_eulerian_at_knot_cons_axisym_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 ...
const double & pe_st() const
Peclet number multiplied by Strouhal number.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative. Note: By default, ALE is enabled, at the expense of possibly creating unnecessary work in problems where the mesh is, in fact, stationary.
unsigned self_test()
Self-test: Return 0 for OK.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and the element Jacobian matrix (wrapper) ...
virtual unsigned u_index_cons_axisym_adv_diff() const
A class for all elements that solve the Advection Diffusion equations in conservative form using isop...
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
virtual void fill_in_generic_residual_contribution_cons_axisym_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...
GeneralisedAxisymAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
GeneralisedAxisymAdvectionDiffusionWindFctPt Conserved_wind_fct_pt
Pointer to additional (conservative) wind function:
static char t char * s
Definition: cfortran.h:572
virtual void get_source_cons_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(std::ostream &outfile, const unsigned &n_plot)
Output function: r,z,u at n_plot^2 plot points.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
virtual void output_fct(std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,z,,u_exact at nplot^2 plot points (dummy time-dependent version to keep intel co...
QGeneralisedAxisymAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Advection ...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
double dshape_and_dtest_eulerian_cons_axisym_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.
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
Definition: shape.h:549
void output(std::ostream &outfile)
Output with default number of plot points.
GeneralisedAxisymAdvectionDiffusionDiffFctPt Diff_fct_pt
Pointer to diffusivity funciton.
virtual void get_conserved_wind_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Get additional (conservative) wind at (Eulerian) position x and/or local coordinate s...
double du_dt_cons_axisym_adv_diff(const unsigned &n) const
du/dt at local node n.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
double integrate_u()
Integrate the concentration over the element.
GeneralisedAxisymAdvectionDiffusionDiffFctPt & diff_fct_pt()
Access function: Pointer to diffusion function.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
double *& pe_pt()
Pointer to Peclet number.
static double Default_peclet_number
Static default value for the Peclet number.
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: r,z,u_exact at n_plot^2 plot points.
void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. r,z,u_exact at n_plot^2 plot points (Calls the s...
void output(std::ostream &outfile)
Output function: r,z,u.
virtual double dshape_and_dtest_eulerian_cons_axisym_adv_diff(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
double * Pe_pt
Pointer to global Peclet number.
QGeneralisedAxisymAdvectionDiffusionElement(const QGeneralisedAxisymAdvectionDiffusionElement< NNODE_1D > &dummy)
Broken copy constructor.
virtual void get_diff_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &D) const
Get diffusivity tensor at (Eulerian) position x and/or local coordinate s. This function is virtual t...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
void output(FILE *file_pt)
C-style output function: r,z,u.
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228