advection_diffusion_reaction_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_ADV_DIFF_REACT_ELEMENTS_HEADER
32 #define OOMPH_ADV_DIFF_REACT_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 Reaction equations using isoparametric elements.
51 /// \f[
52 /// \tau_{i} \frac{\partial C_{i}}{\partial t}
53 /// + w_{j} \frac{\partial C_{i}}{\partial x_{j}} =
54 /// D_{i}\frac{\partial^2 C_{i}}{\partial x_j^2} =
55 /// - R_{i}(C_{i}) + f_{i}
56 /// \f]
57 /// This contains the generic maths. Shape functions, geometric
58 /// mapping etc. must get implemented in derived class.
59 //=============================================================
60 template <unsigned NREAGENT, unsigned DIM>
62 {
63 
64 public:
65 
66  /// \short Function pointer to source function fct(x,f(x)) --
67  /// x is a Vector!
69  Vector<double> &f);
70 
71  /// \short Function pointer to reaction terms
73  const Vector<double> &c, Vector<double> &R);
74 
75  /// \short Function pointer to derivative of reaction terms
77  const Vector<double> &c, DenseMatrix<double> &dRdC);
78 
79 
80  /// \short Function pointer to wind function fct(x,w(x)) --
81  /// x is a Vector!
82  typedef void (*AdvectionDiffusionReactionWindFctPt)(const double &time,
83  const Vector<double>& x,
84  Vector<double>& wind);
85 
86 
87  /// \short Constructor: Initialise the Source_fct_pt, Wind_fct_pt,
88  /// Reaction_fct_pt to null and initialise the dimensionless
89  /// timescale and diffusion ratios
92  ALE_is_disabled(false)
93  {
94  //Set diffusion coefficients to default
96  //Set timescales to default
98  }
99 
100  /// Broken copy constructor
103  {
104  BrokenCopy::broken_copy("AdvectionDiffusionReactionEquations");
105  }
106 
107  /// Broken assignment operator
109  {
110  BrokenCopy::broken_assign("AdvectionDiffusionReactionEquations");
111  }
112 
113  /// \short Return the index at which the unknown i-th reagent
114  /// is stored. The default value, i, is appropriate for single-physics
115  /// problems, when there are only i variables, the values that satisfies
116  /// the advection-diffusion-reaction equation.
117  /// In derived multi-physics elements, this function should be overloaded
118  /// to reflect the chosen storage scheme. Note that these equations require
119  /// that the unknown is always stored at the same index at each node.
120  virtual inline unsigned c_index_adv_diff_react(const unsigned &i)
121  const {return i;}
122 
123 /// \short dc_r/dt at local node n.
124  /// Uses suitably interpolated value for hanging nodes.
125  double dc_dt_adv_diff_react(const unsigned &n, const unsigned &r) const
126  {
127  // Get the data's timestepper
129 
130  //Initialise dudt
131  double dudt=0.0;
132  //Loop over the timesteps, if there is a non Steady timestepper
133  if (!time_stepper_pt->is_steady())
134  {
135  //Find the index at which the variable is stored
136  const unsigned c_nodal_index = c_index_adv_diff_react(r);
137 
138  // Number of timsteps (past & present)
139  const unsigned n_time = time_stepper_pt->ntstorage();
140 
141  for(unsigned t=0;t<n_time;t++)
142  {
143  dudt += time_stepper_pt->weight(1,t)*nodal_value(t,n,c_nodal_index);
144  }
145  }
146  return dudt;
147  }
148 
149 
150 /// \short dc/dt at local node n.
151  /// Uses suitably interpolated value for hanging nodes.
152  void dc_dt_adv_diff_react(const unsigned &n, Vector<double> &dc_dt) const
153  {
154  // Get the data's timestepper
156 
157  //Initialise to zero
158  for(unsigned r=0;r<NREAGENT;r++) {dc_dt[r] = 0.0;}
159 
160  //Loop over the timesteps, if there is a non Steady timestepper
161  if (!time_stepper_pt->is_steady())
162  {
163  // Number of timsteps (past & present)
164  const unsigned n_time = time_stepper_pt->ntstorage();
165  //Local storage (cache) for the weights
166  double weight[n_time];
167  for(unsigned t=0;t<n_time;t++)
168  {weight[t] = time_stepper_pt->weight(1,t);}
169 
170  //Loop over the reagents
171  for(unsigned r=0;r<NREAGENT;r++)
172  {
173  //Find the index at which the variable is stored
174  const unsigned c_nodal_index = c_index_adv_diff_react(r);
175 
176  for(unsigned t=0;t<n_time;t++)
177  {
178  dc_dt[r] += weight[t]*nodal_value(t,n,c_nodal_index);
179  }
180  }
181  }
182  }
183 
184  /// \short Disable ALE, i.e. assert the mesh is not moving -- you do this
185  /// at your own risk!
186  void disable_ALE()
187  {
188  ALE_is_disabled=true;
189  }
190 
191  /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
192  /// when evaluating the time-derivative. Note: By default, ALE is
193  /// enabled, at the expense of possibly creating unnecessary work
194  /// in problems where the mesh is, in fact, stationary.
195  void enable_ALE()
196  {
197  ALE_is_disabled=false;
198  }
199 
200 
201  /// Output with default number of plot points
202  void output(std::ostream &outfile)
203  {
204  unsigned nplot=5;
205  output(outfile,nplot);
206  }
207 
208  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
209  /// nplot^DIM plot points
210  void output(std::ostream &outfile, const unsigned &nplot);
211 
212 
213  /// C_style output with default number of plot points
214  void output(FILE* file_pt)
215  {
216  unsigned n_plot=5;
217  output(file_pt,n_plot);
218  }
219 
220  /// \short C-style output FE representation of soln: x,y,u or x,y,z,u at
221  /// n_plot^DIM plot points
222  void output(FILE* file_pt, const unsigned &n_plot);
223 
224 
225  /// Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points
226  void output_fct(std::ostream &outfile, const unsigned &nplot,
228  exact_soln_pt);
229 
230  /// \short Output exact soln: x,y,u_exact or x,y,z,u_exact at
231  /// nplot^DIM plot points (dummy time-dependent version to
232  /// keep intel compiler happy)
233  virtual void output_fct(std::ostream &outfile, const unsigned &nplot,
234  const double& time,
236  {
237  throw OomphLibError(
238  "There is no time-dependent output_fct() for Advection Diffusion elements",
239  OOMPH_CURRENT_FUNCTION,
240  OOMPH_EXCEPTION_LOCATION);
241  }
242 
243 
244  /// Get error against and norm of exact solution
245  void compute_error(std::ostream &outfile,
247  exact_soln_pt, double& error, double& norm);
248 
249 
250  /// Dummy, time dependent error checker
251  void compute_error(std::ostream &outfile,
253  exact_soln_pt,
254  const double& time, double& error, double& norm)
255  {
256  throw OomphLibError(
257  "No time-dependent compute_error() for Advection Diffusion elements",
258  OOMPH_CURRENT_FUNCTION,
259  OOMPH_EXCEPTION_LOCATION);
260  }
261 
262 
263  /// Access function: Pointer to source function
265  {return Source_fct_pt;}
266 
267 
268  /// Access function: Pointer to source function. Const version
270  {return Source_fct_pt;}
271 
272 
273  /// Access function: Pointer to wind function
275  {return Wind_fct_pt;}
276 
277  /// Access function: Pointer to reaction function. Const version
279  {return Wind_fct_pt;}
280 
281  /// Access function: Pointer to reaction function
283  {return Reaction_fct_pt;}
284 
285  /// Access function: Pointer to reaction function. Const version
287  {return Reaction_fct_pt;}
288 
289  /// Access function: Pointer to reaction derivatives function
291  {return Reaction_deriv_fct_pt;}
292 
293  /// Access function: Pointer to reaction function. Const version
295  {return Reaction_deriv_fct_pt;}
296 
297  /// Vector of diffusion coefficients
298  const Vector<double> &diff() const {return *Diff_pt;}
299 
300  /// Pointer to vector of diffusion coefficients
302 
303  /// Vector of dimensionless timescales
304  const Vector<double> &tau() const {return *Tau_pt;}
305 
306  /// Pointer to vector of dimensionless timescales
308 
309  /// \short Get source term at (Eulerian) position x. This function is
310  /// virtual to allow overloading in multi-physics problems where
311  /// the strength of the source function might be determined by
312  /// another system of equations
313  inline virtual void get_source_adv_diff_react(const unsigned& ipt,
314  const Vector<double>& x,
315  Vector<double> &source) const
316  {
317  //If no source function has been set, return zero
318  if(Source_fct_pt==0)
319  {
320  for(unsigned r=0;r<NREAGENT;r++){source[r] = 0.0;}
321  }
322  else
323  {
324  // Get source strength
325  (*Source_fct_pt)(x,source);
326  }
327  }
328 
329  /// \short Get wind at (Eulerian) position x and/or local coordinate s.
330  /// This function is
331  /// virtual to allow overloading in multi-physics problems where
332  /// the wind function might be determined by
333  /// another system of equations
334  inline virtual void get_wind_adv_diff_react(const unsigned& ipt,
335  const Vector<double> &s,
336  const Vector<double>& x,
337  Vector<double>& wind) const
338  {
339  //If no wind function has been set, return zero
340  if(Wind_fct_pt==0)
341  {
342  for(unsigned i=0;i<DIM;i++) {wind[i]= 0.0;}
343  }
344  else
345  {
346  // Get continuous time from timestepper of first node
347  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
348 
349  // Get wind
350  (*Wind_fct_pt)(time,x,wind);
351  }
352  }
353 
354 
355  /// \short Get reaction as a function of the given reagent concentrations
356  /// This function is
357  /// virtual to allow overloading in multi-physics problems where
358  /// the reaction function might be determined by
359  /// another system of equations
360  inline virtual void get_reaction_adv_diff_react(const unsigned& ipt,
361  const Vector<double> &C,
362  Vector<double>& R) const
363  {
364  //If no wind function has been set, return zero
365  if(Reaction_fct_pt==0)
366  {
367  for(unsigned r=0;r<NREAGENT;r++) {R[r]= 0.0;}
368  }
369  else
370  {
371  // Get reaction terms
372  (*Reaction_fct_pt)(C,R);
373  }
374  }
375 
376  /// \short Get the derivatives of the reaction terms with respect to the
377  /// concentration variables. If no explicit function pointer is set,
378  /// these will be calculated by finite differences
379  virtual void get_reaction_deriv_adv_diff_react(const unsigned& ipt,
380  const Vector<double> &C,
381  DenseMatrix<double> &dRdC)
382  const
383  {
384  //If no reaction pointer set, return zero
385  if(Reaction_fct_pt==0)
386  {
387  for(unsigned r=0;r<NREAGENT;r++)
388  {
389  for(unsigned p=0;p<NREAGENT;p++)
390  {
391  dRdC(r,p) = 0.0;
392  }
393  }
394  }
395  else
396  {
397  //If no function pointer get finite differences
398  if(Reaction_deriv_fct_pt==0)
399  {
400  //Local copy of the unknowns
401  Vector<double> C_local = C;
402  //Finite differences
403  Vector<double> R(NREAGENT), R_plus(NREAGENT), R_minus(NREAGENT);
404  //Get the initial reaction terms
405  //(*Reaction_fct_pt)(C,R);
406  const double fd_step = GeneralisedElement::Default_fd_jacobian_step;
407  //Now loop over all the reagents
408  for(unsigned p=0;p<NREAGENT;p++)
409  {
410  //Store the old value
411  double old_var = C_local[p];
412  //Increment the value
413  C_local[p] += fd_step;
414  //Get the new values
415  (*Reaction_fct_pt)(C_local,R_plus);
416 
417  //Reset the value
418  C_local[p] = old_var;
419  //Decrement the value
420  C_local[p] -= fd_step;
421  //Get the new values
422  (*Reaction_fct_pt)(C_local,R_minus);
423 
424  //Assemble the column of the jacobian
425  for(unsigned r=0;r<NREAGENT;r++)
426  {
427  dRdC(r,p) = (R_plus[r] - R_minus[r])/(2.0*fd_step);
428  }
429 
430  //Reset the value
431  C_local[p] = old_var;
432  }
433  }
434  //Otherwise get the terms from the function
435  else
436  {
437  (*Reaction_deriv_fct_pt)(C,dRdC);
438  }
439  }
440  }
441 
442  /// Get flux: \f$\mbox{flux}[DIM r + i] = \mbox{d}C_{r} / \mbox{d}x_i \f$
443  void get_flux(const Vector<double>& s, Vector<double>& flux) const
444  {
445  //Find out how many nodes there are in the element
446  const unsigned n_node = nnode();
447 
448  //Set up memory for the shape and test functions
449  Shape psi(n_node);
450  DShape dpsidx(n_node,DIM);
451 
452  //Call the derivatives of the shape and test functions
453  dshape_eulerian(s,psi,dpsidx);
454 
455  //Initialise to zero
456  for(unsigned j=0;j<DIM*NREAGENT;j++) {flux[j] = 0.0;}
457 
458  //Loop over the reagent terms
459  for(unsigned r=0;r<NREAGENT;r++)
460  {
461  unsigned c_nodal_index = c_index_adv_diff_react(r);
462 
463  //Loop over derivative directions
464  for(unsigned j=0;j<DIM;j++)
465  {
466  unsigned index = r*DIM + j;
467  //Loop over the nodes
468  for(unsigned l=0;l<n_node;l++)
469  {
470  flux[index] += nodal_value(l,c_nodal_index)*dpsidx(l,j);
471  }
472  }
473  }
474  }
475 
476 
477  /// Add the element's contribution to its residual vector (wrapper)
479  {
480  //Call the generic residuals function with flag set to 0 and using
481  //a dummy matrix
485  }
486 
487 
488  /// \short Add the element's contribution to its residual vector and
489  /// the element Jacobian matrix (wrapper)
491  DenseMatrix<double> &jacobian)
492  {
493  //Call the generic routine with the flag set to 1
495  residuals,jacobian,GeneralisedElement::Dummy_matrix,1);
496  }
497 
498 
499  /// Add the element's contribution to its residuals vector,
500  /// jacobian matrix and mass matrix
502  Vector<double> &residuals, DenseMatrix<double> &jacobian,
503  DenseMatrix<double> &mass_matrix)
504  {
505  //Call the generic routine with the flag set to 2
507  jacobian,
508  mass_matrix,2);
509  }
510 
511 
512  /// Return FE representation of function value c_i(s) at local coordinate s
514  const unsigned &i) const
515  {
516  //Find number of nodes
517  unsigned n_node = nnode();
518 
519  //Get the nodal index at which the unknown is stored
520  unsigned c_nodal_index = c_index_adv_diff_react(i);
521 
522  //Local shape function
523  Shape psi(n_node);
524 
525  //Find values of shape function
526  shape(s,psi);
527 
528  //Initialise value of u
529  double interpolated_u = 0.0;
530 
531  //Loop over the local nodes and sum
532  for(unsigned l=0;l<n_node;l++)
533  {
534  interpolated_u += nodal_value(l,c_nodal_index)*psi[l];
535  }
536 
537  return(interpolated_u);
538  }
539 
540 
541  /// \short Self-test: Return 0 for OK
542  unsigned self_test();
543 
544  /// \short Return the integrated reagent concentrations
545  void integrate_reagents(Vector<double> &C) const;
546 
547 protected:
548 
549 
550  /// \short Shape/test functions and derivs w.r.t. to global coords at
551  /// local coord. s; return Jacobian of mapping
553  const Vector<double> &s,
554  Shape &psi,
555  DShape &dpsidx,
556  Shape &test,
557  DShape &dtestdx) const=0;
558 
559  /// \short Shape/test functions and derivs w.r.t. to global coords at
560  /// integration point ipt; return Jacobian of mapping
562  const unsigned &ipt,
563  Shape &psi,
564  DShape &dpsidx,
565  Shape &test,
566  DShape &dtestdx)
567  const=0;
568 
569  /// \short Add the element's contribution to its residual vector only
570  /// (if flag=and/or element Jacobian matrix
572  Vector<double> &residuals, DenseMatrix<double> &jacobian,
573  DenseMatrix<double> &mass_matrix, unsigned flag);
574 
575  /// Pointer to global diffusion coefficients
577 
578  /// Pointer to global timescales
580 
581  /// Pointer to source function:
583 
584  /// Pointer to wind function:
586 
587  /// Pointer to reaction function
589 
590  /// Pointer to reaction derivatives
592 
593  /// \short Boolean flag to indicate if ALE formulation is disabled when
594  /// time-derivatives are computed. Only set to true if you're sure
595  /// that the mesh is stationary.
597 
598  private:
599 
600  /// Static default value for the dimensionless numbers
602 
603 
604 };
605 
606 
607 ///////////////////////////////////////////////////////////////////////////
608 ///////////////////////////////////////////////////////////////////////////
609 ///////////////////////////////////////////////////////////////////////////
610 
611 
612 
613 //======================================================================
614 /// \short QAdvectionDiffusionReactionElement elements are
615 /// linear/quadrilateral/brick-shaped Advection Diffusion elements with
616 /// isoparametric interpolation for the function.
617 //======================================================================
618 template <unsigned NREAGENT, unsigned DIM, unsigned NNODE_1D>
620  public virtual QElement<DIM,NNODE_1D>,
621  public virtual AdvectionDiffusionReactionEquations<NREAGENT,DIM>
622  {
623  public:
624 
625 
626  ///\short Constructor: Call constructors for QElement and
627  /// Advection Diffusion equations
630  { }
631 
632  /// Broken copy constructor
635  {
636  BrokenCopy::broken_copy("QAdvectionDiffusionReactionElement");
637  }
638 
639  /// Broken assignment operator
640  void operator=(
642  {
643  BrokenCopy::broken_assign("QAdvectionDiffusionReactionElement");
644  }
645 
646  /// \short Required # of `values' (pinned or dofs)
647  /// at node n
648  inline unsigned required_nvalue(const unsigned &n) const {return NREAGENT;}
649 
650  /// \short Output function:
651  /// x,y,u or x,y,z,u
652  void output(std::ostream &outfile)
654 
655  /// \short Output function:
656  /// x,y,u or x,y,z,u at n_plot^DIM plot points
657  void output(std::ostream &outfile, const unsigned &n_plot)
659 
660 
661  /// \short C-style output function:
662  /// x,y,u or x,y,z,u
663  void output(FILE* file_pt)
664  {
666  }
667 
668  /// \short C-style output function:
669  /// x,y,u or x,y,z,u at n_plot^DIM plot points
670  void output(FILE* file_pt, const unsigned &n_plot)
671  {
673  }
674 
675  /// \short Output function for an exact solution:
676  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
677  void output_fct(std::ostream &outfile, const unsigned &n_plot,
679  exact_soln_pt)
680  {
682  output_fct(outfile,n_plot,exact_soln_pt);}
683 
684 
685  /// \short Output function for a time-dependent exact solution.
686  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
687  /// (Calls the steady version)
688  void output_fct(std::ostream &outfile, const unsigned &n_plot,
689  const double& time,
691  exact_soln_pt)
692  {
694  output_fct(outfile,n_plot,time,exact_soln_pt);
695  }
696 
697 
698 protected:
699 
700  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
702  const Vector<double> &s,
703  Shape &psi,
704  DShape &dpsidx,
705  Shape &test,
706  DShape &dtestdx) const;
707 
708  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
709  /// integration point ipt. Return Jacobian.
711  const unsigned& ipt,
712  Shape &psi,
713  DShape &dpsidx,
714  Shape &test,
715  DShape &dtestdx)
716  const;
717 
718 };
719 
720 //Inline functions:
721 
722 
723 //======================================================================
724 /// \short Define the shape functions and test functions and derivatives
725 /// w.r.t. global coordinates and return Jacobian of mapping.
726 ///
727 /// Galerkin: Test functions = shape functions
728 //======================================================================
729 template<unsigned NREAGENT, unsigned DIM, unsigned NNODE_1D>
732  Shape &psi,
733  DShape &dpsidx,
734  Shape &test,
735  DShape &dtestdx) const
736 {
737  //Call the geometrical shape functions and derivatives
738  double J = this->dshape_eulerian(s,psi,dpsidx);
739 
740  //Loop over the test functions and derivatives and set them equal to the
741  //shape functions
742  for(unsigned i=0;i<NNODE_1D;i++)
743  {
744  test[i] = psi[i];
745  for(unsigned j=0;j<DIM;j++)
746  {
747  dtestdx(i,j) = dpsidx(i,j);
748  }
749  }
750 
751  //Return the jacobian
752  return J;
753 }
754 
755 
756 
757 //======================================================================
758 /// Define the shape functions and test functions and derivatives
759 /// w.r.t. global coordinates and return Jacobian of mapping.
760 ///
761 /// Galerkin: Test functions = shape functions
762 //======================================================================
763 template<unsigned NREAGENT, unsigned DIM, unsigned NNODE_1D>
766  const unsigned &ipt,
767  Shape &psi,
768  DShape &dpsidx,
769  Shape &test,
770  DShape &dtestdx) const
771 {
772  //Call the geometrical shape functions and derivatives
773  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
774 
775  //Set the test functions equal to the shape functions (pointer copy)
776  test = psi;
777  dtestdx = dpsidx;
778 
779  //Return the jacobian
780  return J;
781 }
782 
783 
784 ////////////////////////////////////////////////////////////////////////
785 ////////////////////////////////////////////////////////////////////////
786 ////////////////////////////////////////////////////////////////////////
787 
788 
789 
790 //=======================================================================
791 /// \short Face geometry for the QAdvectionDiffusionReactionElement elements:
792 /// The spatial dimension of the face elements is one lower than that
793 /// of the bulk element but they have the same number of points along
794 /// their 1D edges.
795 //=======================================================================
796 template<unsigned NREAGENT, unsigned DIM, unsigned NNODE_1D>
798  QAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> >:
799 public virtual QElement<DIM-1,NNODE_1D>
800 {
801 
802  public:
803 
804  /// \short Constructor: Call the constructor for the
805  /// appropriate lower-dimensional QElement
806  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
807 
808 };
809 
810 
811 
812 ////////////////////////////////////////////////////////////////////////
813 ////////////////////////////////////////////////////////////////////////
814 ////////////////////////////////////////////////////////////////////////
815 
816 
817 //=======================================================================
818 /// Face geometry for the 1D QAdvectionDiffusionReaction elements: Point elements
819 //=======================================================================
820 template<unsigned NREAGENT,unsigned NNODE_1D>
821 class FaceGeometry<QAdvectionDiffusionReactionElement<NREAGENT,1,NNODE_1D> >:
822  public virtual PointElement
823 {
824 
825  public:
826 
827  /// \short Constructor: Call the constructor for the
828  /// appropriate lower-dimensional QElement
830 
831 };
832 
833 }
834 
835 #endif
virtual void get_wind_adv_diff_react(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...
const Vector< double > & tau() const
Vector of dimensionless timescales.
virtual void fill_in_generic_residual_contribution_adv_diff_react(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...
double dshape_and_dtest_eulerian_at_knot_adv_diff_react(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 broken_copy(const std::string &class_name)
Issue error message and terminate execution.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
void dc_dt_adv_diff_react(const unsigned &n, Vector< double > &dc_dt) const
dc/dt at local node n.
virtual double dshape_and_dtest_eulerian_adv_diff_react(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 dshape_and_dtest_eulerian_adv_diff_react(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 output(FILE *file_pt)
C_style output with default number of plot points.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
double interpolated_c_adv_diff_react(const Vector< double > &s, const unsigned &i) const
Return FE representation of function value c_i(s) at local coordinate s.
AdvectionDiffusionReactionEquations(const AdvectionDiffusionReactionEquations &dummy)
Broken copy constructor.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
virtual void get_reaction_adv_diff_react(const unsigned &ipt, const Vector< double > &C, Vector< double > &R) const
Get reaction as a function of the given reagent concentrations This function is virtual to allow over...
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 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) ...
cstr elem_len * i
Definition: cfortran.h:607
double dc_dt_adv_diff_react(const unsigned &n, const unsigned &r) const
dc_r/dt at local node n.
QAdvectionDiffusionReactionElement elements are linear/quadrilateral/brick-shaped Advection Diffusion...
AdvectionDiffusionReactionReactionDerivFctPt Reaction_deriv_fct_pt
Pointer to reaction derivatives.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
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 &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...
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
virtual unsigned c_index_adv_diff_react(const unsigned &i) const
Return the index at which the unknown i-th reagent is stored. The default value, i, is appropriate for single-physics problems, when there are only i variables, the values that satisfies the advection-diffusion-reaction equation. In derived multi-physics elements, this function should be overloaded to reflect the chosen storage scheme. Note that these equations require that the unknown is always stored at the same index at each node.
void(* AdvectionDiffusionReactionReactionDerivFctPt)(const Vector< double > &c, DenseMatrix< double > &dRdC)
Function pointer to derivative of reaction terms.
A general Finite Element class.
Definition: elements.h:1271
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 output(std::ostream &outfile)
Output with default number of plot points.
AdvectionDiffusionReactionReactionFctPt & reaction_fct_pt()
Access function: Pointer to reaction function.
void operator=(const AdvectionDiffusionReactionEquations &)
Broken assignment operator.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
void integrate_reagents(Vector< double > &C) const
Return the integrated reagent concentrations.
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
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
const Vector< double > & diff() const
Vector of diffusion coefficients.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
AdvectionDiffusionReactionReactionDerivFctPt & reaction_deriv_fct_pt()
Access function: Pointer to reaction derivatives function.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
AdvectionDiffusionReactionEquations()
Constructor: Initialise the Source_fct_pt, Wind_fct_pt, Reaction_fct_pt to null and initialise the di...
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
AdvectionDiffusionReactionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
void(* AdvectionDiffusionReactionWindFctPt)(const double &time, const Vector< double > &x, Vector< double > &wind)
Function pointer to wind function fct(x,w(x)) – x is a Vector!
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 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.
AdvectionDiffusionReactionWindFctPt wind_fct_pt() const
Access function: Pointer to reaction function. Const version.
A class for all elements that solve the Advection Diffusion Reaction equations using isoparametric el...
AdvectionDiffusionReactionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
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.
AdvectionDiffusionReactionWindFctPt Wind_fct_pt
Pointer to wind function:
static char t char * s
Definition: cfortran.h:572
AdvectionDiffusionReactionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
AdvectionDiffusionReactionReactionDerivFctPt reaction_deriv_fct_pt() const
Access function: Pointer to reaction function. Const version.
Vector< double > *& tau_pt()
Pointer to vector of dimensionless timescales.
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
static Vector< double > Default_dimensionless_number
Static default value for the dimensionless numbers.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
Vector< double > * Diff_pt
Pointer to global diffusion coefficients.
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
void(* AdvectionDiffusionReactionReactionFctPt)(const Vector< double > &c, Vector< double > &R)
Function pointer to reaction terms.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
AdvectionDiffusionReactionSourceFctPt Source_fct_pt
Pointer to source function:
QAdvectionDiffusionReactionElement(const QAdvectionDiffusionReactionElement< NREAGENT, DIM, NNODE_1D > &dummy)
Broken copy constructor.
QAdvectionDiffusionReactionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
virtual double dshape_and_dtest_eulerian_at_knot_adv_diff_react(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 ...
Vector< double > * Tau_pt
Pointer to global timescales.
Vector< double > *& diff_pt()
Pointer to vector of diffusion coefficients.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void operator=(const QAdvectionDiffusionReactionElement< NREAGENT, DIM, NNODE_1D > &)
Broken assignment operator.
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(* AdvectionDiffusionReactionSourceFctPt)(const Vector< double > &x, Vector< double > &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
unsigned required_nvalue(const unsigned &n) const
Required # of `values' (pinned or dofs) at node n.
AdvectionDiffusionReactionReactionFctPt Reaction_fct_pt
Pointer to reaction function.
virtual void get_source_adv_diff_react(const unsigned &ipt, const Vector< double > &x, Vector< 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, 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 ...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164
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.
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...
AdvectionDiffusionReactionReactionFctPt reaction_fct_pt() const
Access function: Pointer to reaction function. Const version.
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
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.
virtual void get_reaction_deriv_adv_diff_react(const unsigned &ipt, const Vector< double > &C, DenseMatrix< double > &dRdC) const
Get the derivatives of the reaction terms with respect to the concentration variables. If no explicit function pointer is set, these will be calculated by finite differences.