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 in a spherical polar coordinate
31 //system
32 #ifndef OOMPH_AXISYM_ADV_DIFF_ELEMENTS_HEADER
33 #define OOMPH_AXISYM_ADV_DIFF_ELEMENTS_HEADER
34 
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38 #include <oomph-lib-config.h>
39 #endif
40 
41 //OOMPH-LIB headers
42 #include "../generic/nodes.h"
43 #include "../generic/Qelements.h"
44 #include "../generic/refineable_elements.h"
45 #include "../generic/oomph_utilities.h"
46 
47 namespace oomph
48 {
49 
50 //=============================================================
51 /// \short A class for all elements that solve the
52 /// Advection Diffusion equations in a cylindrical polar coordinate system
53 /// using isoparametric elements.
54 /// \f[
55 /// Pe \mathbf{w}\cdot(\mathbf{x}) \nabla u =
56 /// \nabla \cdot \left( \nabla u \right) + f(\mathbf{x})
57 /// \f]
58 /// This contains the generic maths. Shape functions, geometric
59 /// mapping etc. must get implemented in derived class.
60 //=============================================================
62 {
63 
64 public:
65 
66  /// \short Function pointer to source function fct(x,f(x)) --
67  /// x is a Vector!
69  (const Vector<double>& x, double& f);
70 
71 
72  /// \short Function pointer to wind function fct(x,w(x)) --
73  /// x is a Vector!
75  const Vector<double>& x, Vector<double>& wind);
76 
77 
78  /// \short Constructor: Initialise the Source_fct_pt and Wind_fct_pt
79  /// to null and set (pointer to) Peclet number to default
81  ALE_is_disabled(false)
82  {
83  //Set pointer to Peclet number to the default value zero
87  }
88 
89  /// Broken copy constructor
92  {
93  BrokenCopy::broken_copy("AxisymAdvectionDiffusionEquations");
94  }
95 
96  /// Broken assignment operator
97 //Commented out broken assignment operator because this can lead to a conflict warning
98 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
99 //realise that two separate implementations of the broken function are the same and so,
100 //quite rightly, it shouts.
101  /*void operator=(const AxisymAdvectionDiffusionEquations&)
102  {
103  BrokenCopy::broken_assign("AxisymAdvectionDiffusionEquations");
104  }*/
105 
106  /// \short Return the index at which the unknown value
107  /// is stored. The default value, 0, is appropriate for single-physics
108  /// problems, when there is only one variable, the value that satisfies
109  /// the spherical advection-diffusion equation.
110  /// In derived multi-physics elements, this function should be overloaded
111  /// to reflect the chosen storage scheme. Note that these equations require
112  /// that the unknown is always stored at the same index at each node.
113  virtual inline unsigned u_index_axi_adv_diff() const
114  {
115  return 0;
116  }
117 
118 
119 /// \short du/dt at local node n.
120  /// Uses suitably interpolated value for hanging nodes.
121  double du_dt_axi_adv_diff(const unsigned &n) const
122  {
123  // Get the data's timestepper
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_axi_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  /// \short Number of scalars/fields output by this element. Reimplements
164  /// broken virtual function in base class.
165  unsigned nscalar_paraview() const
166  {
167  return 4;
168  }
169 
170  /// \short Write values of the i-th scalar field at the plot points. Needs
171  /// to be implemented for each new specific element type.
172  void scalar_value_paraview(std::ofstream& file_out,
173  const unsigned& i,
174  const unsigned& nplot) const
175  {
176  // Vector of local coordinates
177  Vector<double> s(2);
178 
179  // Loop over plot points
180  unsigned num_plot_points=nplot_points_paraview(nplot);
181  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
182  {
183 
184  // Get local coordinates of plot point
185  get_s_plot(iplot,nplot,s);
186 
187  // Get Eulerian coordinate of plot point
188  Vector<double> x(2);
189  interpolated_x(s,x);
190 
191  // Winds
192  if(i<3)
193  {
194  //Get the wind
195  Vector<double> wind(3);
196 
197  //Dummy ipt argument needed... ?
198  unsigned ipt = 0;
199  get_wind_axi_adv_diff(ipt,s,x,wind);
200 
201  file_out << wind[i] << std::endl;
202  }
203  // Advection Diffusion
204  else if(i==3)
205  {
206  file_out << this->interpolated_u_axi_adv_diff(s) << std::endl;
207  }
208  // Never get here
209  else
210  {
211  std::stringstream error_stream;
212  error_stream
213  << "Advection Diffusion Elements only store 4 fields " << std::endl;
214  throw OomphLibError(
215  error_stream.str(),
216  OOMPH_CURRENT_FUNCTION,
217  OOMPH_EXCEPTION_LOCATION);
218  }
219  }
220  }
221 
222  /// \short Name of the i-th scalar field. Default implementation
223  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
224  /// overloaded with more meaningful names in specific elements.
225  std::string scalar_name_paraview(const unsigned& i) const
226  {
227  // Winds
228  if(i<3)
229  {
230  return "Wind "+StringConversion::to_string(i);
231  }
232  // Advection Diffusion field
233  else if(i==3)
234  {
235  return "Advection Diffusion";
236  }
237  // Never get here
238  else
239  {
240  std::stringstream error_stream;
241  error_stream
242  << "Advection Diffusion Elements only store 4 fields "
243  << std::endl;
244  throw OomphLibError(
245  error_stream.str(),
246  OOMPH_CURRENT_FUNCTION,
247  OOMPH_EXCEPTION_LOCATION);
248  // Dummy return
249  return " ";
250  }
251  }
252 
253  /// Output with default number of plot points
254  void output(std::ostream &outfile)
255  {
256  unsigned nplot = 5;
257  output(outfile,nplot);
258  }
259 
260  /// \short Output FE representation of soln: r,z,u at
261  /// nplot^2 plot points
262  void output(std::ostream &outfile, const unsigned &nplot);
263 
264 
265  /// C_style output with default number of plot points
266  void output(FILE* file_pt)
267  {
268  unsigned n_plot = 5;
269  output(file_pt,n_plot);
270  }
271 
272  /// \short C-style output FE representation of soln: r,z,u at
273  /// n_plot^2 plot points
274  void output(FILE* file_pt, const unsigned &n_plot);
275 
276 
277  /// Output exact soln: r,z,u_exact at nplot^2 plot points
278  void output_fct(std::ostream &outfile,
279  const unsigned &nplot,
281  exact_soln_pt);
282 
283  /// Get error against and norm of exact solution
284  void compute_error(std::ostream &outfile,
286  exact_soln_pt, double& error,
287  double& norm);
288 
289  /// Access function: Pointer to source function
291  {return Source_fct_pt;}
292 
293 
294  /// Access function: Pointer to source function. Const version
296  {return Source_fct_pt;}
297 
298 
299  /// Access function: Pointer to wind function
301  {return Wind_fct_pt;}
302 
303 
304  /// Access function: Pointer to wind function. Const version
306  {return Wind_fct_pt;}
307 
308  // Access functions for the physical constants
309 
310  /// Peclet number
311  inline const double &pe() const {return *Pe_pt;}
312 
313  /// Pointer to Peclet number
314  inline double* &pe_pt() {return Pe_pt;}
315 
316  /// Peclet number multiplied by Strouhal number
317  inline const double &pe_st() const {return *PeSt_pt;}
318 
319  /// Pointer to Peclet number multipled by Strouha number
320  inline double* &pe_st_pt() {return PeSt_pt;}
321 
322  /// Peclet number multiplied by Strouhal number
323  inline const double &d() const {return *D_pt;}
324 
325  /// Pointer to Peclet number multipled by Strouha number
326  inline double* &d_pt() {return D_pt;}
327 
328  /// \short Get source term at (Eulerian) position x. This function is
329  /// virtual to allow overloading in multi-physics problems where
330  /// the strength of the source function might be determined by
331  /// another system of equations
332  inline virtual void get_source_axi_adv_diff(const unsigned& ipt,
333  const Vector<double>& x,
334  double& source) const
335  {
336  //If no source function has been set, return zero
337  if(Source_fct_pt==0)
338  {
339  source = 0.0;
340  }
341  else
342  {
343  //Get source strength
344  (*Source_fct_pt)(x,source);
345  }
346  }
347 
348  /// \short Get wind at (Eulerian) position x and/or local coordinate s.
349  /// This function is
350  /// virtual to allow overloading in multi-physics problems where
351  /// the wind function might be determined by
352  /// another system of equations
353  inline virtual void get_wind_axi_adv_diff(const unsigned& ipt,
354  const Vector<double> &s,
355  const Vector<double>& x,
356  Vector<double>& wind) const
357  {
358  //If no wind function has been set, return zero
359  if(Wind_fct_pt==0)
360  {
361  for(unsigned i=0;i<3;i++)
362  {
363  wind[i] = 0.0;
364  }
365  }
366  else
367  {
368  //Get wind
369  (*Wind_fct_pt)(x,wind);
370  }
371  }
372 
373  /// \short Get flux:
374  // \f$ \mbox{flux}[i] = \nabla u = \mbox{d}u / \mbox{d} r
375  // + 1/r \mbox{d}u / \mbox{d} \theta \f$
376  void get_flux(const Vector<double>& s, Vector<double>& flux) const
377  {
378  //Find out how many nodes there are in the element
379  const unsigned n_node = nnode();
380 
381  //Get the nodal index at which the unknown is stored
382  const unsigned u_nodal_index = u_index_axi_adv_diff();
383 
384  //Set up memory for the shape and test functions
385  Shape psi(n_node);
386  DShape dpsidx(n_node,2);
387 
388  //Call the derivatives of the shape and test functions
389  dshape_eulerian(s,psi,dpsidx);
390 
391  //Initialise to zero
392  for(unsigned j=0;j<2;j++) {flux[j] = 0.0;}
393 
394  //Loop over nodes
395  for(unsigned l=0;l<n_node;l++)
396  {
397  const double u_value = this->nodal_value(l,u_nodal_index);
398  //Add in the derivative directions
399  flux[0] += u_value*dpsidx(l,0);
400  flux[1] += u_value*dpsidx(l,1);
401  }
402  }
403 
404 
405  /// Add the element's contribution to its residual vector (wrapper)
407  {
408  //Call the generic residuals function with flag set to 0 and using
409  //a dummy matrix
411  residuals,
414  0);
415  }
416 
417 
418  /// \short Add the element's contribution to its residual vector and
419  /// the element Jacobian matrix (wrapper)
421  DenseMatrix<double> &jacobian)
422  {
423  //Call the generic routine with the flag set to 1
425  residuals,
426  jacobian,
428  }
429 
430 
431  /// Return FE representation of function value u(s) at local coordinate s
432  inline double interpolated_u_axi_adv_diff(const Vector<double> &s) const
433  {
434  //Find number of nodes
435  const unsigned n_node = nnode();
436 
437  //Get the nodal index at which the unknown is stored
438  const unsigned u_nodal_index = u_index_axi_adv_diff();
439 
440  //Local shape function
441  Shape psi(n_node);
442 
443  //Find values of shape function
444  shape(s,psi);
445 
446  //Initialise value of u
447  double interpolated_u = 0.0;
448 
449  //Loop over the local nodes and sum
450  for(unsigned l=0;l<n_node;l++)
451  {
452  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
453  }
454 
455  return(interpolated_u);
456  }
457 
458 
459  ///\short Return derivative of u at point s with respect to all data
460  ///that can affect its value.
461  ///In addition, return the global equation numbers corresponding to the
462  ///data. This is virtual so that it can be overloaded in the
463  ///refineable version
465  const Vector<double> &s,
466  Vector<double> &du_ddata,
467  Vector<unsigned> &global_eqn_number)
468  {
469  //Find number of nodes
470  const unsigned n_node = nnode();
471 
472  //Get the nodal index at which the unknown is stored
473  const unsigned u_nodal_index = u_index_axi_adv_diff();
474 
475  //Local shape function
476  Shape psi(n_node);
477 
478  //Find values of shape function
479  shape(s,psi);
480 
481  //Find the number of dofs associated with interpolated u
482  unsigned n_u_dof = 0;
483  for(unsigned l=0;l<n_node;l++)
484  {
485  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
486  //If it's positive add to the count
487  if (global_eqn >=0)
488  {
489  ++n_u_dof;
490  }
491  }
492 
493  //Now resize the storage schemes
494  du_ddata.resize(n_u_dof,0.0);
495  global_eqn_number.resize(n_u_dof,0);
496 
497  //Loop over the nodes again and set the derivatives
498  unsigned count = 0;
499  for(unsigned l=0;l<n_node;l++)
500  {
501  //Get the global equation number
502  int global_eqn=this->node_pt(l)->eqn_number(u_nodal_index);
503  //If it's positive
504  if (global_eqn >= 0)
505  {
506  //Set the global equation number
507  global_eqn_number[count] = global_eqn;
508  //Set the derivative with respect to the unknown
509  du_ddata[count] = psi[l];
510  //Increase the counter
511  ++count;
512  }
513  }
514  }
515 
516 
517  /// \short Self-test: Return 0 for OK
518  unsigned self_test();
519 
520 protected:
521 
522  /// \short Shape/test functions and derivs w.r.t. to global coords at
523  /// local coord. s; return Jacobian of mapping
525  const Vector<double> &s,
526  Shape &psi,
527  DShape &dpsidx,
528  Shape &test,
529  DShape &dtestdx)
530  const = 0;
531 
532  /// \short Shape/test functions and derivs w.r.t. to global coords at
533  /// integration point ipt; return Jacobian of mapping
535  const unsigned &ipt,
536  Shape &psi,
537  DShape &dpsidx,
538  Shape &test,
539  DShape &dtestdx)
540  const = 0;
541 
542  /// \short Add the element's contribution to its residual vector only
543  /// (if flag=and/or element Jacobian matrix
545  Vector<double> &residuals, DenseMatrix<double> &jacobian,
546  DenseMatrix<double> &mass_matrix, unsigned flag);
547 
548  // Physical constants
549 
550  /// Pointer to global Peclet number
551  double *Pe_pt;
552 
553  /// Pointer to global Peclet number multiplied by Strouhal number
554  double *PeSt_pt;
555 
556  /// Pointer to global Diffusion parameter
557  double *D_pt;
558 
559  /// Pointer to source function:
561 
562  /// Pointer to wind function:
564 
565  /// Boolean flag to indicate whether AlE formulation is disable
567 
568  private:
569 
570  /// Static default value for the Peclet number
571  static double Default_peclet_number;
572 
573  /// Static default value for the Peclet number
575 
576 
577 };//End class AxisymAdvectionDiffusionEquations
578 
579 
580 ///////////////////////////////////////////////////////////////////////////
581 ///////////////////////////////////////////////////////////////////////////
582 ///////////////////////////////////////////////////////////////////////////
583 
584 
585 
586 //======================================================================
587 /// \short QAxisymAdvectionDiffusionElement elements are
588 /// linear/quadrilateral/brick-shaped Axisymmetric Advection Diffusion
589 /// elements with isoparametric interpolation for the function.
590 //======================================================================
591 template <unsigned NNODE_1D>
592 class QAxisymAdvectionDiffusionElement : public virtual QElement<2,NNODE_1D>,
594 {
595 
596 private:
597 
598  /// \short Static array of ints to hold number of variables at
599  /// nodes: Initial_Nvalue[n]
600  static const unsigned Initial_Nvalue;
601 
602 public:
603 
604 
605  ///\short Constructor: Call constructors for QElement and
606  /// Advection Diffusion equations
609  { }
610 
611  /// Broken copy constructor
614  {
615  BrokenCopy::broken_copy("QAxisymAdvectionDiffusionElement");
616  }
617 
618  /// Broken assignment operator
619  /*void operator=(const QAxisymAdvectionDiffusionElement<NNODE_1D>&)
620  {
621  BrokenCopy::broken_assign("QAxisymAdvectionDiffusionElement");
622  }*/
623 
624  /// \short Required # of `values' (pinned or dofs)
625  /// at node n
626  inline unsigned required_nvalue(const unsigned &n) const
627  {
628  return Initial_Nvalue;
629  }
630 
631  /// \short Output function:
632  /// r,z,u
633  void output(std::ostream &outfile)
634  {
636  }
637 
638  /// \short Output function:
639  /// r,z,u at n_plot^2 plot points
640  void output(std::ostream &outfile, const unsigned &n_plot)
641  {
643  }
644 
645 
646  /// \short C-style output function:
647  /// r,z,u
648  void output(FILE* file_pt)
649  {
651  }
652 
653  /// \short C-style output function:
654  /// r,z,u at n_plot^2 plot points
655  void output(FILE* file_pt, const unsigned &n_plot)
656  {
658  }
659 
660  /// \short Output function for an exact solution:
661  /// r,z,u_exact at n_plot^2 plot points
662  void output_fct(std::ostream &outfile,
663  const unsigned &n_plot,
665  {
666  AxisymAdvectionDiffusionEquations::output_fct(outfile,n_plot,exact_soln_pt);
667  }
668 
669 
670 protected:
671 
672  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
674  const Vector<double> &s,
675  Shape &psi,
676  DShape &dpsidx,
677  Shape &test,
678  DShape &dtestdx) const;
679 
680  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
681  /// integration point ipt. Return Jacobian.
683  const unsigned& ipt,
684  Shape &psi,
685  DShape &dpsidx,
686  Shape &test,
687  DShape &dtestdx) const;
688 
689 };//End class QAxisymAdvectionDiffusionElement
690 
691 //Inline functions:
692 
693 //======================================================================
694 /// \short Define the shape functions and test functions and derivatives
695 /// w.r.t. global coordinates and return Jacobian of mapping.
696 ///
697 /// Galerkin: Test functions = shape functions
698 //======================================================================
699 
700 template<unsigned NNODE_1D>
703  Shape &psi,
704  DShape &dpsidx,
705  Shape &test,
706  DShape &dtestdx) const
707 {
708  //Call the geometrical shape functions and derivatives
709  double J = this->dshape_eulerian(s,psi,dpsidx);
710 
711  //Loop over the test functions and derivatives and set them equal to the
712  //shape functions
713  for(unsigned i=0;i<NNODE_1D;i++)
714  {
715  test[i] = psi[i];
716  for(unsigned j=0;j<2;j++)
717  {
718  dtestdx(i,j) = dpsidx(i,j);
719  }
720  }
721 
722  //Return the jacobian
723  return J;
724 }
725 
726 
727 
728 //======================================================================
729 /// \short Define the shape functions and test functions and derivatives
730 /// w.r.t. global coordinates and return Jacobian of mapping.
731 ///
732 /// Galerkin: Test functions = shape functions
733 //======================================================================
734 
735 template<unsigned NNODE_1D>
738  Shape &psi,
739  DShape &dpsidx,
740  Shape &test,
741  DShape &dtestdx) const
742 {
743  //Call the geometrical shape functions and derivatives
744  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
745 
746  //Set the test functions equal to the shape functions (pointer copy)
747  test = psi;
748  dtestdx = dpsidx;
749 
750  //Return the jacobian
751  return J;
752 }
753 
754 ////////////////////////////////////////////////////////////////////////
755 ////////////////////////////////////////////////////////////////////////
756 ////////////////////////////////////////////////////////////////////////
757 
758 template<unsigned NNODE_1D>
760 public virtual QElement<1,NNODE_1D>
761 {
762 
763  public:
764 
765  /// \short Constructor: Call the constructor for the
766  /// appropriate lower-dimensional QElement
767  FaceGeometry() : QElement<1,NNODE_1D>() {}
768 
769 };
770 
771 
772 
773 
774 ////////////////////////////////////////////////////////////////////////
775 ////////////////////////////////////////////////////////////////////////
776 ////////////////////////////////////////////////////////////////////////
777 
778 //======================================================================
779 /// \short A class for elements that allow the imposition of an
780 /// applied Robin boundary condition on the boundaries of Steady
781 /// Axisymmnetric Advection Diffusion Flux elements.
782 /// \f[
783 /// -\Delta u \cdot \mathbf{n} + \alpha(r,z) u = \beta(r,z)
784 /// \f]
785 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
786 /// policy class.
787 //======================================================================
788 template <class ELEMENT>
790 public virtual FaceGeometry<ELEMENT>,
791  public virtual FaceElement
792 {
793 
794 public:
795 
796 
797  /// \short Function pointer to the prescribed-beta function fct(x,beta(x)) --
798  /// x is a Vector!
800  const Vector<double>& x,
801  double& beta);
802 
803  /// \short Function pointer to the prescribed-alpha function fct(x,alpha(x)) --
804  /// x is a Vector!
806  const Vector<double>& x,
807  double& alpha);
808 
809 
810  /// \short Constructor, takes the pointer to the "bulk" element
811  /// and the index of the face to be created
813  const int &face_index);
814 
815 
816  /// Broken empty constructor
818  {
819  throw OomphLibError(
820  "Don't call empty constructor for AxisymAdvectionDiffusionFluxElement",
821  OOMPH_CURRENT_FUNCTION,
822  OOMPH_EXCEPTION_LOCATION);
823  }
824 
825  /// Broken copy constructor
828  {
829  BrokenCopy::broken_copy("AxisymAdvectionDiffusionFluxElement");
830  }
831 
832  /// Broken assignment operator
833  /*void operator=(const AxisymAdvectionDiffusionFluxElement&)
834  {
835  BrokenCopy::broken_assign("AxisymAdvectionDiffusionFluxElement");
836  }*/
837 
838  /// Access function for the prescribed-beta function pointer
840  {
841  return Beta_fct_pt;
842  }
843 
844  /// Access function for the prescribed-alpha function pointer
846  {
847  return Alpha_fct_pt;
848  }
849 
850 
851  /// Add the element's contribution to its residual vector
853  {
854  //Call the generic residuals function with flag set to 0
855  //using a dummy matrix
857  residuals,
859  }
860 
861 
862  /// \short Add the element's contribution to its residual vector and
863  /// its Jacobian matrix
865  DenseMatrix<double> &jacobian)
866  {
867  //Call the generic routine with the flag set to 1
869  jacobian,
870  1);
871  }
872 
873  /// \short Specify the value of nodal zeta from the face geometry
874  /// The "global" intrinsic coordinate of the element when
875  /// viewed as part of a geometric object should be given by
876  /// the FaceElement representation, by default (needed to break
877  /// indeterminacy if bulk element is SolidElement)
878  double zeta_nodal(const unsigned &n, const unsigned &k,
879  const unsigned &i) const
880  {return FaceElement::zeta_nodal(n,k,i);}
881 
882 
883 
884  /// \short Output function -- forward to broken version in FiniteElement
885  /// until somebody decides what exactly they want to plot here...
886  void output(std::ostream &outfile)
887  {
888  FiniteElement::output(outfile);
889  }
890 
891  /// \short Output function -- forward to broken version in FiniteElement
892  /// until somebody decides what exactly they want to plot here...
893  void output(std::ostream &outfile, const unsigned &nplot)
894  {
895  FiniteElement::output(outfile,nplot);
896  }
897 
898 
899 protected:
900 
901  /// \short Function to compute the shape and test functions and to return
902  /// the Jacobian of mapping between local and global (Eulerian)
903  /// coordinates
904  inline double shape_and_test(const Vector<double> &s,
905  Shape &psi,
906  Shape &test) const
907  {
908  //Find number of nodes
909  unsigned n_node = nnode();
910 
911  //Get the shape functions
912  shape(s,psi);
913 
914  //Set the test functions to be the same as the shape functions
915  for(unsigned i=0;i<n_node;i++)
916  {
917  test[i] = psi[i];
918  }
919 
920  //Return the value of the jacobian
921  return J_eulerian(s);
922  }
923 
924 
925  /// \short Function to compute the shape and test functions and to return
926  /// the Jacobian of mapping between local and global (Eulerian)
927  /// coordinates
928  inline double shape_and_test_at_knot(const unsigned &ipt,
929  Shape &psi,
930  Shape &test) const
931  {
932  //Find number of nodes
933  unsigned n_node = nnode();
934 
935  //Get the shape functions
936  shape_at_knot(ipt,psi);
937 
938  //Set the test functions to be the same as the shape functions
939  for(unsigned i=0;i<n_node;i++)
940  {
941  test[i] = psi[i];
942  }
943 
944  //Return the value of the jacobian
945  return J_eulerian_at_knot(ipt);
946  }
947 
948  /// \short Function to calculate the prescribed beta at a given spatial
949  /// position
950  void get_beta(const Vector<double>& x, double& beta)
951  {
952  //If the function pointer is zero return zero
953  if(Beta_fct_pt == 0)
954  {
955  beta = 0.0;
956  }
957  //Otherwise call the function
958  else
959  {
960  (*Beta_fct_pt)(x,beta);
961  }
962  }
963 
964  /// \short Function to calculate the prescribed alpha at a given spatial
965  /// position
966  void get_alpha(const Vector<double>& x, double& alpha)
967  {
968  //If the function pointer is zero return zero
969  if(Alpha_fct_pt == 0)
970  {
971  alpha = 0.0;
972  }
973  //Otherwise call the function
974  else
975  {
976  (*Alpha_fct_pt)(x,alpha);
977  }
978  }
979 
980 private:
981 
982 
983  /// \short Add the element's contribution to its residual vector.
984  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
986  Vector<double> &residuals, DenseMatrix<double> &jacobian,
987  unsigned flag);
988 
989 
990  /// Function pointer to the (global) prescribed-beta function
992 
993  /// Function pointer to the (global) prescribed-alpha function
995 
996  /// The index at which the unknown is stored at the nodes
998 
999 
1000 };//End class AxisymAdvectionDiffusionFluxElement
1001 
1002 
1003 ///////////////////////////////////////////////////////////////////////
1004 ///////////////////////////////////////////////////////////////////////
1005 ///////////////////////////////////////////////////////////////////////
1006 
1007 
1008 //===========================================================================
1009 /// \short Constructor, takes the pointer to the "bulk" element and the index
1010 /// of the face to be created
1011 //===========================================================================
1012 template<class ELEMENT>
1015  const int &face_index) :
1016 FaceGeometry<ELEMENT>(), FaceElement()
1017 
1018 {
1019  //Let the bulk element build the FaceElement, i.e. setup the pointers
1020  //to its nodes (by referring to the appropriate nodes in the bulk
1021  //element), etc.
1022  bulk_el_pt->build_face_element(face_index,this);
1023 
1024 #ifdef PARANOID
1025  {
1026  //Check that the element is not a refineable 3d element
1027  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
1028  //If it's three-d
1029  if(elem_pt->dim()==3)
1030  {
1031  //Is it refineable
1032  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
1033  if(ref_el_pt!=0)
1034  {
1035  if (this->has_hanging_nodes())
1036  {
1037  throw OomphLibError(
1038  "This flux element will not work correctly if nodes are hanging\n",
1039  OOMPH_CURRENT_FUNCTION,
1040  OOMPH_EXCEPTION_LOCATION);
1041  }
1042  }
1043  }
1044  }
1045 #endif
1046 
1047  //Initialise the prescribed-beta function pointer to zero
1048  Beta_fct_pt = 0;
1049 
1050  //Set up U_index_adv_diff. Initialise to zero, which probably won't change
1051  //in most cases, oh well, the price we pay for generality
1052  U_index_adv_diff = 0;
1053 
1054  //Cast to the appropriate AdvectionDiffusionEquation so that we can
1055  //find the index at which the variable is stored
1056  //We assume that the dimension of the full problem is the same
1057  //as the dimension of the node, if this is not the case you will have
1058  //to write custom elements, sorry
1060  dynamic_cast<AxisymAdvectionDiffusionEquations*>(bulk_el_pt);
1061 
1062  //If the cast has failed die
1063  if(eqn_pt==0)
1064  {
1065  std::string error_string =
1066  "Bulk element must inherit from AxisymAdvectionDiffusionEquations.";
1067  error_string +=
1068  "Nodes are two dimensional, but cannot cast the bulk element to\n";
1069  error_string += "AxisymAdvectionDiffusionEquations<2>\n.";
1070  error_string +=
1071  "If you desire this functionality, you must implement it yourself\n";
1072 
1073  throw OomphLibError(
1074  error_string,
1075  OOMPH_CURRENT_FUNCTION,
1076  OOMPH_EXCEPTION_LOCATION);
1077  }
1078  else
1079  {
1080  //Read the index from the (cast) bulk element.
1082  }
1083 
1084 
1085 }
1086 
1087 
1088 //===========================================================================
1089 /// \short Compute the element's residual vector and the (zero) Jacobian
1090 /// matrix for the Robin boundary condition:
1091 /// \f[
1092 /// \Delta u \cdot \mathbf{n} + \alpha (\mathbf{x}) = \beta (\mathbf{x})
1093 /// \f]
1094 //===========================================================================
1095 template<class ELEMENT>
1098  Vector<double> &residuals,
1099  DenseMatrix<double> &jacobian,
1100  unsigned flag)
1101 {
1102  //Find out how many nodes there are
1103  const unsigned n_node = nnode();
1104 
1105  // Locally cache the index at which the variable is stored
1106  const unsigned u_index_axi_adv_diff = U_index_adv_diff;
1107 
1108  //Set up memory for the shape and test functions
1109  Shape psif(n_node), testf(n_node);
1110 
1111  //Set the value of n_intpt
1112  const unsigned n_intpt = integral_pt()->nweight();
1113 
1114  //Set the Vector to hold local coordinates
1115  Vector<double> s(1);
1116 
1117  //Integers used to store the local equation number and local unknown
1118  //indices for the residuals and jacobians
1119  int local_eqn=0, local_unknown=0;
1120 
1121  //Loop over the integration points
1122  //--------------------------------
1123  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1124  {
1125 
1126  //Assign values of s
1127  for(unsigned i=0;i<1;i++)
1128  {
1129  s[i] = integral_pt()->knot(ipt,i);
1130  }
1131 
1132  //Get the integral weight
1133  double w = integral_pt()->weight(ipt);
1134 
1135  //Find the shape and test functions and return the Jacobian
1136  //of the mapping
1137  double J = shape_and_test(s,psif,testf);
1138 
1139  //Premultiply the weights and the Jacobian
1140  double W = w*J;
1141 
1142  //Calculate local values of the solution and its derivatives
1143  //Allocate
1144  double interpolated_u=0.0;
1145  Vector<double> interpolated_x(2,0.0);
1146 
1147  //Calculate position
1148  for(unsigned l=0;l<n_node;l++)
1149  {
1150  //Get the value at the node
1151  double u_value = raw_nodal_value(l,u_index_axi_adv_diff);
1152  interpolated_u += u_value*psif(l);
1153  //Loop over coordinate direction
1154  for(unsigned i=0;i<2;i++)
1155  {
1156  interpolated_x[i] += nodal_position(l,i)*psif(l);
1157  }
1158  }
1159 
1160  //Get the imposed beta (beta=flux when alpha=0.0)
1161  double beta;
1162  get_beta(interpolated_x,beta);
1163 
1164  //Get the imposed alpha
1165  double alpha;
1166  get_alpha(interpolated_x,alpha);
1167 
1168  //calculate the area weighting dS = r^{2} sin theta dr dtheta
1169  // r = x[0] and theta = x[1]
1170  double dS = interpolated_x[0]*interpolated_x[0]*sin(interpolated_x[1]);
1171 
1172  //Now add to the appropriate equations
1173 
1174  //Loop over the test functions
1175  for(unsigned l=0;l<n_node;l++)
1176  {
1177  //Set the local equation number
1178  local_eqn = nodal_local_eqn(l,u_index_axi_adv_diff);
1179  /*IF it's not a boundary condition*/
1180  if(local_eqn >= 0)
1181  {
1182  //Add the prescribed beta terms
1183  residuals[local_eqn] -= dS*(beta - alpha*interpolated_u)*testf(l)*W;
1184 
1185  //Calculate the Jacobian
1186  //----------------------
1187  if ( (flag) && (alpha!=0.0) )
1188  {
1189  //Loop over the velocity shape functions again
1190  for(unsigned l2=0;l2<n_node;l2++)
1191  {
1192  //Set the number of the unknown
1193  local_unknown = nodal_local_eqn(l2,u_index_axi_adv_diff);
1194 
1195  //If at a non-zero degree of freedom add in the entry
1196  if(local_unknown >= 0)
1197  {
1198  jacobian(local_eqn,local_unknown) +=
1199  dS*alpha*psif[l2]*testf[l]*W;
1200  }
1201  }
1202  }
1203 
1204  }
1205  } //end loop over test functions
1206 
1207  } //end loop over integration points
1208 
1209 }//end fill_in_generic_residual_contribution_adv_diff_flux
1210 
1211 
1212 } //End AxisymAdvectionDiffusionFluxElement class
1213 
1214 #endif
QAxisymAdvectionDiffusionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
double *& d_pt()
Pointer to Peclet number multipled by Strouha number.
void(* AxisymAdvectionDiffusionPrescribedBetaFctPt)(const Vector< double > &x, double &beta)
Function pointer to the prescribed-beta function fct(x,beta(x)) – x is a Vector!
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
double dshape_and_dtest_eulerian_at_knot_axi_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.
AxisymAdvectionDiffusionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
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.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
AxisymAdvectionDiffusionFluxElement(const AxisymAdvectionDiffusionFluxElement &dummy)
Broken copy constructor.
void get_alpha(const Vector< double > &x, double &alpha)
Function to calculate the prescribed alpha at a given spatial position.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5121
double dshape_and_dtest_eulerian_axi_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 output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,z,u at n_plot^2 plot points.
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
virtual void get_source_axi_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...
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class...
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
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
bool ALE_is_disabled
Boolean flag to indicate whether AlE formulation is disable.
QAxisymAdvectionDiffusionElement(const QAxisymAdvectionDiffusionElement< NNODE_1D > &dummy)
Broken copy constructor.
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
AxisymAdvectionDiffusionPrescribedAlphaFctPt Alpha_fct_pt
Function pointer to the (global) prescribed-alpha function.
AxisymAdvectionDiffusionEquations(const AxisymAdvectionDiffusionEquations &dummy)
Broken copy constructor.
const double & pe_st() const
Peclet number multiplied by Strouhal number.
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
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
void output(std::ostream &outfile)
Output function: r,z,u.
A general Finite Element class.
Definition: elements.h:1271
void output(std::ostream &outfile, const unsigned &nplot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
virtual double dshape_and_dtest_eulerian_axi_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...
virtual double dshape_and_dtest_eulerian_at_knot_axi_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 ...
virtual void get_wind_axi_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 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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
void(* AxisymAdvectionDiffusionPrescribedAlphaFctPt)(const Vector< double > &x, double &alpha)
Function pointer to the prescribed-alpha function fct(x,alpha(x)) – x is a Vector! ...
double * Pe_pt
Pointer to global Peclet number.
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2700
AxisymAdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
virtual void fill_in_generic_residual_contribution_axi_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 output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,z,u at n_plot^2 plot points.
const double & d() const
Peclet number multiplied by Strouhal number.
AxisymAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
unsigned U_index_adv_diff
The index at which the unknown is stored at the nodes.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
QAxisymAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Axisymmetric Advectio...
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) ...
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.
virtual unsigned u_index_axi_adv_diff() const
Broken assignment operator.
AxisymAdvectionDiffusionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.
void output(std::ostream &outfile)
Output with default number of plot points.
void(* AxisymAdvectionDiffusionWindFctPt)(const Vector< double > &x, Vector< double > &wind)
Function pointer to wind function fct(x,w(x)) – x is a Vector!
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
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
double * D_pt
Pointer to global Diffusion parameter.
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
AxisymAdvectionDiffusionPrescribedAlphaFctPt & alpha_fct_pt()
Access function for the prescribed-alpha function pointer.
static char t char * s
Definition: cfortran.h:572
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
A class for elements that allow the imposition of an applied Robin boundary condition on the boundari...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
AxisymAdvectionDiffusionPrescribedBetaFctPt Beta_fct_pt
Function pointer to the (global) prescribed-beta function.
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3158
void fill_in_generic_residual_contribution_axi_adv_diff_flux(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Add the element's contribution to its residual vector. flag=1(or 0): do (or don't) compute the Jacobi...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
virtual void dinterpolated_u_axi_adv_diff_ddata(const Vector< double > &s, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
Return derivative of u at point s with respect to all data that can affect its value. In addition, return the global equation numbers corresponding to the data. This is virtual so that it can be overloaded in the refineable version.
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:5033
void(* AxisymAdvectionDiffusionSourceFctPt)(const Vector< double > &x, double &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
void output(FILE *file_pt)
C-style output function: r,z,u.
AxisymAdvectionDiffusionPrescribedBetaFctPt & beta_fct_pt()
Broken assignment operator.
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
AxisymAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
static double Default_peclet_number
Static default value for the Peclet number.
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux:
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
AxisymAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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: r,z,u_exact at nplot^2 plot points.
AxisymAdvectionDiffusionEquations()
Constructor: Initialise the Source_fct_pt and Wind_fct_pt to null and set (pointer to) Peclet number ...
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2344
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
static double Default_diffusion_parameter
Static default value for the Peclet number.
double interpolated_u_axi_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
void get_beta(const Vector< double > &x, double &beta)
Function to calculate the prescribed beta at a given spatial position.
void output(std::ostream &outfile)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
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...
A class for all elements that solve the Advection Diffusion equations in a cylindrical polar coordina...
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
double du_dt_axi_adv_diff(const unsigned &n) const
du/dt at local node n.
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types)...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and its Jacobian matrix.