gen_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: 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 //Header file for Advection Diffusion elements
31 #ifndef OOMPH_GEN_ADV_DIFF_ELEMENTS_HEADER
32 #define OOMPH_GEN_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 /// \short A class for all elements that solve the Advection
50 /// Diffusion equations in conservative form using isoparametric elements.
51 /// \f[
52 /// \frac{\partial}{\partial x_{i}}\left(
53 /// Pe w_{i}(x_{k}) u - D_{ij}(x_{k})\frac{\partial u}{\partial x_{j}}\right)
54 /// = f(x_{j})
55 /// \f]
56 /// This contains the generic maths. Shape functions, geometric
57 /// mapping etc. must get implemented in derived class.
58 //=============================================================
59 template <unsigned DIM>
61 {
62 
63 public:
64 
65  /// \short Function pointer to source function fct(x,f(x)) --
66  /// x is a Vector!
68  (const Vector<double>& x, double& f);
69 
70  /// \short Function pointer to wind function fct(x,w(x)) --
71  /// x is a Vector!
73  (const Vector<double>& x, Vector<double>& wind);
74 
75 
76  /// \short Funciton pointer to a diffusivity function
79 
80  /// \short Constructor: Initialise the Source_fct_pt and Wind_fct_pt
81  /// to null and set (pointer to) Peclet number to default
84  {
85  //Set Peclet number to default
87  //Set Peclet Strouhal number to default
89  }
90 
91  /// Broken copy constructor
94  {
95  BrokenCopy::broken_copy("GeneralisedAdvectionDiffusionEquations");
96  }
97 
98  /// Broken assignment operator
100  {
101  BrokenCopy::broken_assign("GeneralisedAdvectionDiffusionEquations");
102  }
103 
104  /// \short Return the index at which the unknown value
105  /// is stored. The default value, 0, is appropriate for single-physics
106  /// problems, when there is only one variable, the value that satisfies
107  /// the advection-diffusion equation.
108  /// In derived multi-physics elements, this function should be overloaded
109  /// to reflect the chosen storage scheme. Note that these equations require
110  /// that the unknown is always stored at the same index at each node.
111  virtual inline unsigned u_index_cons_adv_diff() const {return 0;}
112 
113 /// \short du/dt at local node n.
114  /// Uses suitably interpolated value for hanging nodes.
115  double du_dt_cons_adv_diff(const unsigned &n) const
116  {
117  // Get the data's timestepper
119 
120  //Initialise dudt
121  double dudt=0.0;
122  //Loop over the timesteps, if there is a non Steady timestepper
123  if (!time_stepper_pt->is_steady())
124  {
125  //Find the index at which the variable is stored
126  const unsigned u_nodal_index = u_index_cons_adv_diff();
127 
128  // Number of timsteps (past & present)
129  const unsigned n_time = time_stepper_pt->ntstorage();
130 
131  for(unsigned t=0;t<n_time;t++)
132  {
133  dudt += time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
134  }
135  }
136  return dudt;
137  }
138 
139  /// \short Disable ALE, i.e. assert the mesh is not moving -- you do this
140  /// at your own risk!
141  void disable_ALE()
142  {
143  ALE_is_disabled=true;
144  }
145 
146 
147  /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
148  /// when evaluating the time-derivative. Note: By default, ALE is
149  /// enabled, at the expense of possibly creating unnecessary work
150  /// in problems where the mesh is, in fact, stationary.
151  void enable_ALE()
152  {
153  ALE_is_disabled=false;
154  }
155 
156 
157  /// Output with default number of plot points
158  void output(std::ostream &outfile)
159  {
160  unsigned nplot=5;
161  output(outfile,nplot);
162  }
163 
164  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
165  /// nplot^DIM plot points
166  void output(std::ostream &outfile, const unsigned &nplot);
167 
168 
169  /// C_style output with default number of plot points
170  void output(FILE* file_pt)
171  {
172  unsigned n_plot=5;
173  output(file_pt,n_plot);
174  }
175 
176  /// \short C-style output FE representation of soln: x,y,u or x,y,z,u at
177  /// n_plot^DIM plot points
178  void output(FILE* file_pt, const unsigned &n_plot);
179 
180 
181  /// Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points
182  void output_fct(std::ostream &outfile, const unsigned &nplot,
184  exact_soln_pt);
185 
186  /// \short Output exact soln: x,y,u_exact or x,y,z,u_exact at
187  /// nplot^DIM plot points (dummy time-dependent version to
188  /// keep intel compiler happy)
189  virtual void output_fct(std::ostream &outfile, const unsigned &nplot,
190  const double& time,
192  {
193  throw OomphLibError(
194  "There is no time-dependent output_fct() for Advection Diffusion elements",
195  OOMPH_CURRENT_FUNCTION,
196  OOMPH_EXCEPTION_LOCATION);
197  }
198 
199 
200  /// Get error against and norm of exact solution
201  void compute_error(std::ostream &outfile,
203  exact_soln_pt, double& error, double& norm);
204 
205 
206  /// Dummy, time dependent error checker
207  void compute_error(std::ostream &outfile,
209  exact_soln_pt,
210  const double& time, double& error, double& norm)
211  {
212  throw OomphLibError(
213  "No time-dependent compute_error() for Advection Diffusion elements",
214  OOMPH_CURRENT_FUNCTION,
215  OOMPH_EXCEPTION_LOCATION);
216  }
217 
218  /// \short Integrate the concentration over the element
219  double integrate_u();
220 
221 
222  /// Access function: Pointer to source function
224  {return Source_fct_pt;}
225 
226 
227  /// Access function: Pointer to source function. Const version
229  {return Source_fct_pt;}
230 
231 
232  /// Access function: Pointer to wind function
234 
235 
236  /// Access function: Pointer to wind function. Const version
238  {return Wind_fct_pt;}
239 
240 
241  /// Access function: Pointer to additional (conservative) wind function
243  {return Conserved_wind_fct_pt;}
244 
245 
246  /// \short Access function: Pointer to additoinal (conservative)
247  /// wind function.
248  /// Const version
250  {return Conserved_wind_fct_pt;}
251 
252  /// Access function: Pointer to diffusion function
254  {return Diff_fct_pt;}
255 
256  /// Access function: Pointer to diffusion function. Const version
258  {return Diff_fct_pt;}
259 
260  /// Peclet number
261  const double &pe() const {return *Pe_pt;}
262 
263  /// Pointer to Peclet number
264  double* &pe_pt() {return Pe_pt;}
265 
266  /// Peclet number multiplied by Strouhal number
267  const double &pe_st() const {return *PeSt_pt;}
268 
269  /// Pointer to Peclet number multipled by Strouha number
270  double* &pe_st_pt() {return PeSt_pt;}
271 
272  /// \short Get source term at (Eulerian) position x. This function is
273  /// virtual to allow overloading in multi-physics problems where
274  /// the strength of the source function might be determined by
275  /// another system of equations
276  inline virtual void get_source_cons_adv_diff(const unsigned& ipt,
277  const Vector<double>& x,
278  double& source) const
279  {
280  //If no source function has been set, return zero
281  if(Source_fct_pt==0) {source = 0.0;}
282  else
283  {
284  // Get source strength
285  (*Source_fct_pt)(x,source);
286  }
287  }
288 
289  /// \short Get wind at (Eulerian) position x and/or local coordinate s.
290  /// This function is
291  /// virtual to allow overloading in multi-physics problems where
292  /// the wind function might be determined by
293  /// another system of equations
294  inline virtual void get_wind_cons_adv_diff(const unsigned& ipt,
295  const Vector<double> &s,
296  const Vector<double>& x,
297  Vector<double>& wind) const
298  {
299  //If no wind function has been set, return zero
300  if(Wind_fct_pt==0)
301  {
302  for(unsigned i=0;i<DIM;i++) {wind[i]= 0.0;}
303  }
304  else
305  {
306  // Get wind
307  (*Wind_fct_pt)(x,wind);
308  }
309  }
310 
311 
312  /// \short Get additional (conservative)
313  /// wind at (Eulerian) position x and/or local coordinate s.
314  /// This function is
315  /// virtual to allow overloading in multi-physics problems where
316  /// the wind function might be determined by
317  /// another system of equations
318  inline virtual void get_conserved_wind_cons_adv_diff(
319  const unsigned& ipt,
320  const Vector<double> &s,
321  const Vector<double>& x,
322  Vector<double>& wind) const
323  {
324  //If no wind function has been set, return zero
325  if(Conserved_wind_fct_pt==0)
326  {
327  for(unsigned i=0;i<DIM;i++) {wind[i]= 0.0;}
328  }
329  else
330  {
331  // Get wind
332  (*Conserved_wind_fct_pt)(x,wind);
333  }
334  }
335 
336 
337  /// \short Get diffusivity tensor at (Eulerian) position
338  /// x and/or local coordinate s.
339  /// This function is
340  /// virtual to allow overloading in multi-physics problems where
341  /// the wind function might be determined by
342  /// another system of equations
343  inline virtual void get_diff_cons_adv_diff(const unsigned& ipt,
344  const Vector<double> &s,
345  const Vector<double>& x,
346  DenseMatrix<double>& D) const
347  {
348  //If no wind function has been set, return identity
349  if(Diff_fct_pt==0)
350  {
351  for(unsigned i=0;i<DIM;i++)
352  {
353  for(unsigned j=0;j<DIM;j++)
354  {
355  if(i==j) {D(i,j) = 1.0;} else {D(i,j) = 0.0;}
356  }
357  }
358  }
359  else
360  {
361  // Get diffusivity tensor
362  (*Diff_fct_pt)(x,D);
363  }
364  }
365 
366 
367  /// Get flux: \f$\mbox{flux}[i] = \mbox{d}u / \mbox{d}x_i \f$
368  void get_flux(const Vector<double>& s, Vector<double>& flux) const
369  {
370  //Find out how many nodes there are in the element
371  unsigned n_node = nnode();
372 
373  //Get the nodal index at which the unknown is stored
374  unsigned u_nodal_index = u_index_cons_adv_diff();
375 
376  //Set up memory for the shape and test functions
377  Shape psi(n_node);
378  DShape dpsidx(n_node,DIM);
379 
380  //Call the derivatives of the shape and test functions
381  dshape_eulerian(s,psi,dpsidx);
382 
383  //Initialise to zero
384  for(unsigned j=0;j<DIM;j++) {flux[j] = 0.0;}
385 
386  // Loop over nodes
387  for(unsigned l=0;l<n_node;l++)
388  {
389  //Loop over derivative directions
390  for(unsigned j=0;j<DIM;j++)
391  {
392  flux[j] += nodal_value(l,u_nodal_index)*dpsidx(l,j);
393  }
394  }
395  }
396 
397  /// Get flux: \f$\mbox{flux}[i] = \mbox{d}u / \mbox{d}x_i \f$
399  Vector<double>& total_flux) const
400  {
401  //Find out how many nodes there are in the element
402  const unsigned n_node = nnode();
403 
404  //Get the nodal index at which the unknown is stored
405  const unsigned u_nodal_index = u_index_cons_adv_diff();
406 
407  //Set up memory for the shape and test functions
408  Shape psi(n_node);
409  DShape dpsidx(n_node,DIM);
410 
411  //Call the derivatives of the shape and test functions
412  dshape_eulerian(s,psi,dpsidx);
413 
414  //Storage for the Eulerian position
416  //Storage for the concentration
417  double interpolated_u = 0.0;
418  //Storage for the derivatives of the concentration
419  Vector<double> interpolated_dudx(DIM,0.0);
420 
421  // Loop over nodes
422  for(unsigned l=0;l<n_node;l++)
423  {
424  //Get the value at the node
425  const double u_value = this->nodal_value(l,u_nodal_index);
426  interpolated_u += u_value*psi(l);
427  // Loop over directions
428  for(unsigned j=0;j<DIM;j++)
429  {
430  interpolated_x[j] += this->nodal_position(l,j)*psi(l);
431  interpolated_dudx[j] += u_value*dpsidx(l,j);
432  }
433  }
434 
435  //Dummy integration point
436  unsigned ipt=0;
437 
438  //Get the conserved wind (non-divergence free)
439  Vector<double> conserved_wind(DIM);
440  get_conserved_wind_cons_adv_diff(ipt,s,interpolated_x,conserved_wind);
441 
442  //Get diffusivity tensor
443  DenseMatrix<double> D(DIM,DIM);
444  get_diff_cons_adv_diff(ipt,s,interpolated_x,D);
445 
446  //Calculate the total flux made up of the diffusive flux
447  //and the conserved wind
448  for(unsigned i=0;i<DIM;i++)
449  {
450  total_flux[i] = 0.0;
451  for(unsigned j=0;j<DIM;j++)
452  {
453  total_flux[i] += D(i,j)*interpolated_dudx[j];
454  }
455  total_flux[i] -= conserved_wind[i]*interpolated_u;
456  }
457  }
458 
459 
460  /// Add the element's contribution to its residual vector (wrapper)
462  {
463  //Call the generic residuals function with flag set to 0 and using
464  //a dummy matrix
468  }
469 
470 
471  /// \short Add the element's contribution to its residual vector and
472  /// the element Jacobian matrix (wrapper)
474  DenseMatrix<double> &jacobian)
475  {
476  //Call the generic routine with the flag set to 1
478  residuals,jacobian,GeneralisedElement::Dummy_matrix,1);
479  }
480 
481 
482  /// Add the element's contribution to its residuals vector,
483  /// jacobian matrix and mass matrix
485  Vector<double> &residuals, DenseMatrix<double> &jacobian,
486  DenseMatrix<double> &mass_matrix)
487  {
488  //Call the generic routine with the flag set to 2
490  jacobian,mass_matrix,2);
491  }
492 
493 
494  /// Return FE representation of function value u(s) at local coordinate s
495  inline double interpolated_u_cons_adv_diff(const Vector<double> &s) const
496  {
497  //Find number of nodes
498  unsigned n_node = nnode();
499 
500  //Get the nodal index at which the unknown is stored
501  unsigned u_nodal_index = u_index_cons_adv_diff();
502 
503  //Local shape function
504  Shape psi(n_node);
505 
506  //Find values of shape function
507  shape(s,psi);
508 
509  //Initialise value of u
510  double interpolated_u = 0.0;
511 
512  //Loop over the local nodes and sum
513  for(unsigned l=0;l<n_node;l++)
514  {
515  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
516  }
517 
518  return(interpolated_u);
519  }
520 
521 
522  /// \short Self-test: Return 0 for OK
523  unsigned self_test();
524 
525 protected:
526 
527 
528  /// \short Shape/test functions and derivs w.r.t. to global coords at
529  /// local coord. s; return Jacobian of mapping
531  Shape &psi,
532  DShape &dpsidx,
533  Shape &test,
534  DShape &dtestdx) const=0;
535 
536  /// \short Shape/test functions and derivs w.r.t. to global coords at
537  /// integration point ipt; return Jacobian of mapping
539  const unsigned &ipt,
540  Shape &psi,
541  DShape &dpsidx,
542  Shape &test,
543  DShape &dtestdx)
544  const=0;
545 
546  /// \short Add the element's contribution to its residual vector only
547  /// (if flag=and/or element Jacobian matrix
549  Vector<double> &residuals, DenseMatrix<double> &jacobian,
550  DenseMatrix<double> &mass_matrix, unsigned flag);
551 
552  /// Pointer to global Peclet number
553  double *Pe_pt;
554 
555  /// Pointer to global Peclet number multiplied by Strouhal number
556  double *PeSt_pt;
557 
558  /// Pointer to source function:
560 
561  /// Pointer to wind function:
563 
564  /// Pointer to additional (conservative) wind function:
566 
567  /// Pointer to diffusivity funciton
569 
570  /// \short Boolean flag to indicate if ALE formulation is disabled when
571  /// time-derivatives are computed. Only set to false if you're sure
572  /// that the mesh is stationary.
574 
575  private:
576 
577  /// Static default value for the Peclet number
578  static double Default_peclet_number;
579 
580 
581 };
582 
583 
584 ///////////////////////////////////////////////////////////////////////////
585 ///////////////////////////////////////////////////////////////////////////
586 ///////////////////////////////////////////////////////////////////////////
587 
588 
589 
590 //======================================================================
591 /// \short QGeneralisedAdvectionDiffusionElement elements are
592 /// linear/quadrilateral/brick-shaped Advection Diffusion elements with
593 /// isoparametric interpolation for the function.
594 //======================================================================
595 template <unsigned DIM, unsigned NNODE_1D>
597  public virtual QElement<DIM,NNODE_1D>,
598  public virtual GeneralisedAdvectionDiffusionEquations<DIM>
599 {
600 
601 private:
602 
603  /// \short Static array of ints to hold number of variables at
604  /// nodes: Initial_Nvalue[n]
605  static const unsigned Initial_Nvalue;
606 
607  public:
608 
609 
610  ///\short Constructor: Call constructors for QElement and
611  /// Advection Diffusion equations
614  { }
615 
616  /// Broken copy constructor
619  {
620  BrokenCopy::broken_copy("QGeneralisedAdvectionDiffusionElement");
621  }
622 
623  /// Broken assignment operator
625  {
626  BrokenCopy::broken_assign("QGeneralisedAdvectionDiffusionElement");
627  }
628 
629  /// \short Required # of `values' (pinned or dofs)
630  /// at node n
631  inline unsigned required_nvalue(const unsigned &n) const
632  {return Initial_Nvalue;}
633 
634  /// \short Output function:
635  /// x,y,u or x,y,z,u
636  void output(std::ostream &outfile)
638 
639  /// \short Output function:
640  /// x,y,u or x,y,z,u at n_plot^DIM plot points
641  void output(std::ostream &outfile, const unsigned &n_plot)
643 
644 
645  /// \short C-style output function:
646  /// x,y,u or x,y,z,u
647  void output(FILE* file_pt)
648  {
650  }
651 
652  /// \short C-style output function:
653  /// x,y,u or x,y,z,u at n_plot^DIM plot points
654  void output(FILE* file_pt, const unsigned &n_plot)
655  {
657  }
658 
659  /// \short Output function for an exact solution:
660  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
661  void output_fct(std::ostream &outfile, const unsigned &n_plot,
663  exact_soln_pt)
664  {GeneralisedAdvectionDiffusionEquations<DIM>::output_fct(outfile,n_plot,exact_soln_pt);}
665 
666 
667  /// \short Output function for a time-dependent exact solution.
668  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
669  /// (Calls the steady version)
670  void output_fct(std::ostream &outfile, const unsigned &n_plot,
671  const double& time,
673  exact_soln_pt)
674  {
676  output_fct(outfile,n_plot,time,exact_soln_pt);
677  }
678 
679 
680 protected:
681 
682  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
684  const Vector<double> &s,
685  Shape &psi,
686  DShape &dpsidx,
687  Shape &test,
688  DShape &dtestdx) const;
689 
690  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
691  /// integration point ipt. Return Jacobian.
693  const unsigned& ipt,
694  Shape &psi,
695  DShape &dpsidx,
696  Shape &test,
697  DShape &dtestdx)
698  const;
699 
700 };
701 
702 //Inline functions:
703 
704 
705 //======================================================================
706 /// \short Define the shape functions and test functions and derivatives
707 /// w.r.t. global coordinates and return Jacobian of mapping.
708 ///
709 /// Galerkin: Test functions = shape functions
710 //======================================================================
711 template<unsigned DIM, unsigned NNODE_1D>
714  Shape &psi,
715  DShape &dpsidx,
716  Shape &test,
717  DShape &dtestdx) const
718 {
719  //Call the geometrical shape functions and derivatives
720  double J = this->dshape_eulerian(s,psi,dpsidx);
721 
722  //Loop over the test functions and derivatives and set them equal to the
723  //shape functions
724  for(unsigned i=0;i<NNODE_1D;i++)
725  {
726  test[i] = psi[i];
727  for(unsigned j=0;j<DIM;j++)
728  {
729  dtestdx(i,j) = dpsidx(i,j);
730  }
731  }
732 
733  //Return the jacobian
734  return J;
735 }
736 
737 
738 
739 //======================================================================
740 /// Define the shape functions and test functions and derivatives
741 /// w.r.t. global coordinates and return Jacobian of mapping.
742 ///
743 /// Galerkin: Test functions = shape functions
744 //======================================================================
745 template<unsigned DIM, unsigned NNODE_1D>
748  const unsigned &ipt,
749  Shape &psi,
750  DShape &dpsidx,
751  Shape &test,
752  DShape &dtestdx) const
753 {
754  //Call the geometrical shape functions and derivatives
755  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
756 
757  //Set the test functions equal to the shape functions (pointer copy)
758  test = psi;
759  dtestdx = dpsidx;
760 
761  //Return the jacobian
762  return J;
763 }
764 
765 
766 ////////////////////////////////////////////////////////////////////////
767 ////////////////////////////////////////////////////////////////////////
768 ////////////////////////////////////////////////////////////////////////
769 
770 
771 
772 //=======================================================================
773 /// \short Face geometry for the QGeneralisedAdvectionDiffusionElement elements:
774 /// The spatial dimension of the face elements is one lower than that
775 /// of the bulk element but they have the same number of points along
776 /// their 1D edges.
777 //=======================================================================
778 template<unsigned DIM, unsigned NNODE_1D>
780  public virtual QElement<DIM-1,NNODE_1D>
781 {
782 
783  public:
784 
785  /// \short Constructor: Call the constructor for the
786  /// appropriate lower-dimensional QElement
787  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
788 
789 };
790 
791 
792 
793 ////////////////////////////////////////////////////////////////////////
794 ////////////////////////////////////////////////////////////////////////
795 ////////////////////////////////////////////////////////////////////////
796 
797 
798 //=======================================================================
799 /// Face geometry for the 1D QGeneralisedAdvectionDiffusion elements: Point elements
800 //=======================================================================
801 template<unsigned NNODE_1D>
803  public virtual PointElement
804 {
805 
806  public:
807 
808  /// \short Constructor: Call the constructor for the
809  /// appropriate lower-dimensional QElement
811 
812 };
813 
814 }
815 
816 #endif
GeneralisedAdvectionDiffusionEquations()
Constructor: Initialise the Source_fct_pt and Wind_fct_pt to null and set (pointer to) Peclet number ...
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
void output(FILE *file_pt)
C_style output with default number of plot points.
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.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
virtual double dshape_and_dtest_eulerian_cons_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 interpolated_u_cons_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
virtual void get_conserved_wind_cons_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 dshape_and_dtest_eulerian_cons_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.
GeneralisedAdvectionDiffusionDiffFctPt diff_fct_pt() const
Access function: Pointer to diffusion function. Const version.
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
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.
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...
void get_total_flux(const Vector< double > &s, Vector< double > &total_flux) const
Get flux: .
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
double integrate_u()
Integrate the concentration over the element.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
A general Finite Element class.
Definition: elements.h:1271
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
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
GeneralisedAdvectionDiffusionDiffFctPt & diff_fct_pt()
Access function: Pointer to diffusion function.
QGeneralisedAdvectionDiffusionElement(const QGeneralisedAdvectionDiffusionElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.
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
virtual void get_source_cons_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 &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points...
void operator=(const QGeneralisedAdvectionDiffusionElement< DIM, NNODE_1D > &)
Broken assignment operator.
A class for all elements that solve the Advection Diffusion equations in conservative form using isop...
void operator=(const GeneralisedAdvectionDiffusionEquations &)
Broken assignment operator.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
GeneralisedAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
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
GeneralisedAdvectionDiffusionWindFctPt conserved_wind_fct_pt() const
Access function: Pointer to additoinal (conservative) wind function. Const version.
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
void(* GeneralisedAdvectionDiffusionSourceFctPt)(const Vector< double > &x, double &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
virtual void output_fct(std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points (dummy time-dependent versio...
void(* GeneralisedAdvectionDiffusionWindFctPt)(const Vector< double > &x, Vector< double > &wind)
Function pointer to wind function fct(x,w(x)) – x is a Vector!
GeneralisedAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
static char t char * s
Definition: cfortran.h:572
virtual void get_wind_cons_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...
double dshape_and_dtest_eulerian_at_knot_cons_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.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
unsigned required_nvalue(const unsigned &n) const
Required # of `values' (pinned or dofs) at node n.
virtual double dshape_and_dtest_eulerian_at_knot_cons_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 ...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
GeneralisedAdvectionDiffusionWindFctPt Conserved_wind_fct_pt
Pointer to additional (conservative) wind function:
static double Default_peclet_number
Static default value for the Peclet number.
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
QGeneralisedAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Advection Diffus...
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. x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot ...
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
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 const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
GeneralisedAdvectionDiffusionDiffFctPt Diff_fct_pt
Pointer to diffusivity funciton.
GeneralisedAdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
GeneralisedAdvectionDiffusionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
GeneralisedAdvectionDiffusionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
GeneralisedAdvectionDiffusionEquations(const GeneralisedAdvectionDiffusionEquations &dummy)
Broken copy constructor.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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
double du_dt_cons_adv_diff(const unsigned &n) const
du/dt at local node n.
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
GeneralisedAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
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) ...
QGeneralisedAdvectionDiffusionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
virtual void get_diff_cons_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(* GeneralisedAdvectionDiffusionDiffFctPt)(const Vector< double > &x, DenseMatrix< double > &D)
Funciton pointer to a diffusivity function.
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
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
const double & pe_st() const
Peclet number multiplied by Strouhal number.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
GeneralisedAdvectionDiffusionWindFctPt & conserved_wind_fct_pt()
Access function: Pointer to additional (conservative) wind function.
virtual unsigned u_index_cons_adv_diff() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...