spherical_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 in a spherical polar coordinate
31 //system
32 #ifndef OOMPH_SPHERICAL_ADV_DIFF_ELEMENTS_HEADER
33 #define OOMPH_SPHERICAL_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 spherical 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  {
82  //Set pointer to Peclet number to the default value zero
85  }
86 
87  /// Broken copy constructor
90  {
91  BrokenCopy::broken_copy("SphericalAdvectionDiffusionEquations");
92  }
93 
94  /// Broken assignment operator
96  {
97  BrokenCopy::broken_assign("SphericalAdvectionDiffusionEquations");
98  }
99 
100  /// \short Return the index at which the unknown value
101  /// is stored. The default value, 0, is appropriate for single-physics
102  /// problems, when there is only one variable, the value that satisfies
103  /// the spherical advection-diffusion equation.
104  /// In derived multi-physics elements, this function should be overloaded
105  /// to reflect the chosen storage scheme. Note that these equations require
106  /// that the unknown is always stored at the same index at each node.
107  virtual inline unsigned u_index_spherical_adv_diff() const
108  {
109  return 0;
110  }
111 
112 
113 /// \short du/dt at local node n.
114  /// Uses suitably interpolated value for hanging nodes.
115  double du_dt_spherical_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_spherical_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 
140  /// \short Disable ALE -- empty overload to suppress warning.
141  /// ALE isn't implemented anyway
142  void disable_ALE() {}
143 
144  /// Output with default number of plot points
145  void output(std::ostream &outfile)
146  {
147  unsigned nplot = 5;
148  output(outfile,nplot);
149  }
150 
151  /// \short Output FE representation of soln: r,z,u at
152  /// nplot^2 plot points
153  void output(std::ostream &outfile, const unsigned &nplot);
154 
155 
156  /// C_style output with default number of plot points
157  void output(FILE* file_pt)
158  {
159  unsigned n_plot = 5;
160  output(file_pt,n_plot);
161  }
162 
163  /// \short C-style output FE representation of soln: r,z,u at
164  /// n_plot^2 plot points
165  void output(FILE* file_pt, const unsigned &n_plot);
166 
167 
168  /// Output exact soln: r,z,u_exact at nplot^2 plot points
169  void output_fct(std::ostream &outfile,
170  const unsigned &nplot,
172  exact_soln_pt);
173 
174  /// Get error against and norm of exact solution
175  void compute_error(std::ostream &outfile,
177  exact_soln_pt, double& error,
178  double& norm);
179 
180  /// Access function: Pointer to source function
182  {return Source_fct_pt;}
183 
184 
185  /// Access function: Pointer to source function. Const version
187  {return Source_fct_pt;}
188 
189 
190  /// Access function: Pointer to wind function
192  {return Wind_fct_pt;}
193 
194 
195  /// Access function: Pointer to wind function. Const version
197  {return Wind_fct_pt;}
198 
199  // Access functions for the physical constants
200 
201  /// Peclet number
202  inline const double &pe() const {return *Pe_pt;}
203 
204  /// Pointer to Peclet number
205  inline double* &pe_pt() {return Pe_pt;}
206 
207  /// Peclet number multiplied by Strouhal number
208  inline const double &pe_st() const {return *PeSt_pt;}
209 
210  /// Pointer to Peclet number multipled by Strouha number
211  inline double* &pe_st_pt() {return PeSt_pt;}
212 
213  /// \short Get source term at (Eulerian) position x. This function is
214  /// virtual to allow overloading in multi-physics problems where
215  /// the strength of the source function might be determined by
216  /// another system of equations
217  inline virtual void get_source_spherical_adv_diff(const unsigned& ipt,
218  const Vector<double>& x,
219  double& source) const
220  {
221  //If no source function has been set, return zero
222  if(Source_fct_pt==0)
223  {
224  source = 0.0;
225  }
226  else
227  {
228  //Get source strength
229  (*Source_fct_pt)(x,source);
230  }
231  }
232 
233  /// \short Get wind at (Eulerian) position x and/or local coordinate s.
234  /// This function is
235  /// virtual to allow overloading in multi-physics problems where
236  /// the wind function might be determined by
237  /// another system of equations
238  inline virtual void get_wind_spherical_adv_diff(const unsigned& ipt,
239  const Vector<double> &s,
240  const Vector<double>& x,
241  Vector<double>& wind) const
242  {
243  //If no wind function has been set, return zero
244  if(Wind_fct_pt==0)
245  {
246  for(unsigned i=0;i<3;i++)
247  {
248  wind[i] = 0.0;
249  }
250  }
251  else
252  {
253  //Get wind
254  (*Wind_fct_pt)(x,wind);
255  }
256  }
257 
258  /// \short Get flux:
259  /// \f[
260  /// \mbox{flux}[i] = \nabla u = \mbox{d}u / \mbox{d} r
261  /// + 1/r \mbox{d}u / \mbox{d} \theta
262  /// \f]
263  void get_flux(const Vector<double>& s, Vector<double>& flux) const
264  {
265  //Find out how many nodes there are in the element
266  const unsigned n_node = nnode();
267 
268  //Get the nodal index at which the unknown is stored
269  const unsigned u_nodal_index = u_index_spherical_adv_diff();
270 
271  //Set up memory for the shape and test functions
272  Shape psi(n_node);
273  DShape dpsidx(n_node,2);
274 
275  //Call the derivatives of the shape and test functions
276  dshape_eulerian(s,psi,dpsidx);
277 
278  //Initialise to zero
279  for(unsigned j=0;j<2;j++) {flux[j] = 0.0;}
280 
281  //Loop over nodes
282  for(unsigned l=0;l<n_node;l++)
283  {
284  const double u_value = this->nodal_value(l,u_nodal_index);
285  const double r = this->nodal_position(l,0);
286  //Add in the derivative directions
287  flux[0] += u_value*dpsidx(l,0);
288  flux[1] += u_value*dpsidx(l,1)/r;
289  }
290  }
291 
292 
293  /// Add the element's contribution to its residual vector (wrapper)
295  {
296  //Call the generic residuals function with flag set to 0 and using
297  //a dummy matrix
299  residuals,
302  0);
303  }
304 
305 
306  /// \short Add the element's contribution to its residual vector and
307  /// the element Jacobian matrix (wrapper)
309  DenseMatrix<double> &jacobian)
310  {
311  //Call the generic routine with the flag set to 1
313  residuals,
314  jacobian,
316  }
317 
318  /// \short Add the element's contribution to its residual vector and
319  /// the element Jacobian matrix (wrapper) and mass matrix
321  Vector<double> &residuals,
322  DenseMatrix<double> &jacobian,
323  DenseMatrix<double> &mass_matrix)
324  {
325  //Call the generic routine with the flag set to 2
327  residuals,
328  jacobian,
329  mass_matrix,2);
330  }
331 
332 
333 
334  /// Return FE representation of function value u(s) at local coordinate s
336  {
337  //Find number of nodes
338  const unsigned n_node = nnode();
339 
340  //Get the nodal index at which the unknown is stored
341  const unsigned u_nodal_index = u_index_spherical_adv_diff();
342 
343  //Local shape function
344  Shape psi(n_node);
345 
346  //Find values of shape function
347  shape(s,psi);
348 
349  //Initialise value of u
350  double interpolated_u = 0.0;
351 
352  //Loop over the local nodes and sum
353  for(unsigned l=0;l<n_node;l++)
354  {
355  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
356  }
357 
358  return(interpolated_u);
359  }
360 
361 
362  ///\short Return derivative of u at point s with respect to all data
363  ///that can affect its value.
364  ///In addition, return the global equation numbers corresponding to the
365  ///data. This is virtual so that it can be overloaded in the
366  ///refineable version
368  const Vector<double> &s,
369  Vector<double> &du_ddata,
370  Vector<unsigned> &global_eqn_number)
371  {
372  //Find number of nodes
373  const unsigned n_node = nnode();
374 
375  //Get the nodal index at which the unknown is stored
376  const unsigned u_nodal_index = u_index_spherical_adv_diff();
377 
378  //Local shape function
379  Shape psi(n_node);
380 
381  //Find values of shape function
382  shape(s,psi);
383 
384  //Find the number of dofs associated with interpolated u
385  unsigned n_u_dof = 0;
386  for(unsigned l=0;l<n_node;l++)
387  {
388  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
389  //If it's positive add to the count
390  if (global_eqn >=0)
391  {
392  ++n_u_dof;
393  }
394  }
395 
396  //Now resize the storage schemes
397  du_ddata.resize(n_u_dof,0.0);
398  global_eqn_number.resize(n_u_dof,0);
399 
400  //Loop over the nodes again and set the derivatives
401  unsigned count = 0;
402  for(unsigned l=0;l<n_node;l++)
403  {
404  //Get the global equation number
405  int global_eqn=this->node_pt(l)->eqn_number(u_nodal_index);
406  //If it's positive
407  if (global_eqn >= 0)
408  {
409  //Set the global equation number
410  global_eqn_number[count] = global_eqn;
411  //Set the derivative with respect to the unknown
412  du_ddata[count] = psi[l];
413  //Increase the counter
414  ++count;
415  }
416  }
417  }
418 
419 
420  /// \short Self-test: Return 0 for OK
421  unsigned self_test();
422 
423 protected:
424 
425  /// \short Shape/test functions and derivs w.r.t. to global coords at
426  /// local coord. s; return Jacobian of mapping
428  const Vector<double> &s,
429  Shape &psi,
430  DShape &dpsidx,
431  Shape &test,
432  DShape &dtestdx)
433  const = 0;
434 
435  /// \short Shape/test functions and derivs w.r.t. to global coords at
436  /// integration point ipt; return Jacobian of mapping
438  const unsigned &ipt,
439  Shape &psi,
440  DShape &dpsidx,
441  Shape &test,
442  DShape &dtestdx)
443  const = 0;
444 
445  /// \short Add the element's contribution to its residual vector only
446  /// (if flag=and/or element Jacobian matrix
448  Vector<double> &residuals, DenseMatrix<double> &jacobian,
449  DenseMatrix<double> &mass_matrix, unsigned flag);
450 
451  // Physical constants
452 
453  /// Pointer to global Peclet number
454  double *Pe_pt;
455 
456  /// Pointer to global Peclet number multiplied by Strouhal number
457  double *PeSt_pt;
458 
459  /// Pointer to source function:
461 
462  /// Pointer to wind function:
464 
465  private:
466 
467  /// Static default value for the Peclet number
468  static double Default_peclet_number;
469 
470 
471 };//End class SphericalAdvectionDiffusionEquations
472 
473 
474 ///////////////////////////////////////////////////////////////////////////
475 ///////////////////////////////////////////////////////////////////////////
476 ///////////////////////////////////////////////////////////////////////////
477 
478 
479 
480 //======================================================================
481 /// \short QSphericalAdvectionDiffusionElement elements are
482 /// linear/quadrilateral/brick-shaped Axisymmetric Advection Diffusion
483 /// elements with isoparametric interpolation for the function.
484 //======================================================================
485 template <unsigned NNODE_1D>
486 class QSphericalAdvectionDiffusionElement : public virtual QElement<2,NNODE_1D>,
488 {
489 
490 private:
491 
492  /// \short Static array of ints to hold number of variables at
493  /// nodes: Initial_Nvalue[n]
494  static const unsigned Initial_Nvalue;
495 
496 public:
497 
498 
499  ///\short Constructor: Call constructors for QElement and
500  /// Advection Diffusion equations
503  { }
504 
505  /// Broken copy constructor
508  {
509  BrokenCopy::broken_copy("QSphericalAdvectionDiffusionElement");
510  }
511 
512  /// Broken assignment operator
514  {
515  BrokenCopy::broken_assign("QSphericalAdvectionDiffusionElement");
516  }
517 
518  /// \short Required # of `values' (pinned or dofs)
519  /// at node n
520  inline unsigned required_nvalue(const unsigned &n) const
521  {
522  return Initial_Nvalue;
523  }
524 
525  /// \short Output function:
526  /// r,z,u
527  void output(std::ostream &outfile)
528  {
530  }
531 
532  /// \short Output function:
533  /// r,z,u at n_plot^2 plot points
534  void output(std::ostream &outfile, const unsigned &n_plot)
535  {
537  }
538 
539 
540  /// \short C-style output function:
541  /// r,z,u
542  void output(FILE* file_pt)
543  {
545  }
546 
547  /// \short C-style output function:
548  /// r,z,u at n_plot^2 plot points
549  void output(FILE* file_pt, const unsigned &n_plot)
550  {
552  }
553 
554  /// \short Output function for an exact solution:
555  /// r,z,u_exact at n_plot^2 plot points
556  void output_fct(std::ostream &outfile,
557  const unsigned &n_plot,
559  {
560  SphericalAdvectionDiffusionEquations::output_fct(outfile,n_plot,exact_soln_pt);
561  }
562 
563 
564 protected:
565 
566  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
568  const Vector<double> &s,
569  Shape &psi,
570  DShape &dpsidx,
571  Shape &test,
572  DShape &dtestdx) const;
573 
574  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
575  /// integration point ipt. Return Jacobian.
577  const unsigned& ipt,
578  Shape &psi,
579  DShape &dpsidx,
580  Shape &test,
581  DShape &dtestdx) const;
582 
583 };//End class QSphericalAdvectionDiffusionElement
584 
585 //Inline functions:
586 
587 //======================================================================
588 /// \short Define the shape functions and test functions and derivatives
589 /// w.r.t. global coordinates and return Jacobian of mapping.
590 ///
591 /// Galerkin: Test functions = shape functions
592 //======================================================================
593 
594 template<unsigned NNODE_1D>
597  Shape &psi,
598  DShape &dpsidx,
599  Shape &test,
600  DShape &dtestdx) const
601 {
602  //Call the geometrical shape functions and derivatives
603  double J = this->dshape_eulerian(s,psi,dpsidx);
604 
605  //Loop over the test functions and derivatives and set them equal to the
606  //shape functions
607  for(unsigned i=0;i<NNODE_1D;i++)
608  {
609  test[i] = psi[i];
610  for(unsigned j=0;j<2;j++)
611  {
612  dtestdx(i,j) = dpsidx(i,j);
613  }
614  }
615 
616  //Return the jacobian
617  return J;
618 }
619 
620 
621 
622 //======================================================================
623 /// \short Define the shape functions and test functions and derivatives
624 /// w.r.t. global coordinates and return Jacobian of mapping.
625 ///
626 /// Galerkin: Test functions = shape functions
627 //======================================================================
628 
629 template<unsigned NNODE_1D>
632  Shape &psi,
633  DShape &dpsidx,
634  Shape &test,
635  DShape &dtestdx) const
636 {
637  //Call the geometrical shape functions and derivatives
638  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
639 
640  //Set the test functions equal to the shape functions (pointer copy)
641  test = psi;
642  dtestdx = dpsidx;
643 
644  //Return the jacobian
645  return J;
646 }
647 
648 ////////////////////////////////////////////////////////////////////////
649 ////////////////////////////////////////////////////////////////////////
650 ////////////////////////////////////////////////////////////////////////
651 
652 template<unsigned NNODE_1D>
654 public virtual QElement<1,NNODE_1D>
655 {
656 
657  public:
658 
659  /// \short Constructor: Call the constructor for the
660  /// appropriate lower-dimensional QElement
661  FaceGeometry() : QElement<1,NNODE_1D>() {}
662 
663 };
664 
665 
666 
667 
668 ////////////////////////////////////////////////////////////////////////
669 ////////////////////////////////////////////////////////////////////////
670 ////////////////////////////////////////////////////////////////////////
671 
672 //======================================================================
673 /// \short A class for elements that allow the imposition of an
674 /// applied Robin boundary condition on the boundaries of Steady
675 /// Axisymmnetric Advection Diffusion Flux elements.
676 /// \f[
677 /// -\Delta u \cdot \mathbf{n} + \alpha(r,z) u = \beta(r,z)
678 /// \f]
679 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
680 /// policy class.
681 //======================================================================
682 template <class ELEMENT>
684 public virtual FaceGeometry<ELEMENT>,
685  public virtual FaceElement
686 {
687 
688 public:
689 
690 
691  /// \short Function pointer to the prescribed-beta function fct(x,beta(x)) --
692  /// x is a Vector!
694  const Vector<double>& x,
695  double& beta);
696 
697  /// \short Function pointer to the prescribed-alpha function fct(x,alpha(x)) --
698  /// x is a Vector!
700  const Vector<double>& x,
701  double& alpha);
702 
703 
704  /// \short Constructor, takes the pointer to the "bulk" element
705  /// and the index of the face to be created
707  const int &face_index);
708 
709 
710  /// Broken empty constructor
712  {
713  throw OomphLibError(
714  "Don't call empty constructor for SphericalAdvectionDiffusionFluxElement",
715  OOMPH_CURRENT_FUNCTION,
716  OOMPH_EXCEPTION_LOCATION);
717  }
718 
719  /// Broken copy constructor
722  {
723  BrokenCopy::broken_copy("SphericalAdvectionDiffusionFluxElement");
724  }
725 
726  /// Broken assignment operator
728  {
729  BrokenCopy::broken_assign("SphericalAdvectionDiffusionFluxElement");
730  }
731 
732  /// Access function for the prescribed-beta function pointer
734  {
735  return Beta_fct_pt;
736  }
737 
738  /// Access function for the prescribed-alpha function pointer
740  {
741  return Alpha_fct_pt;
742  }
743 
744 
745  /// Add the element's contribution to its residual vector
747  {
748  //Call the generic residuals function with flag set to 0
749  //using a dummy matrix
751  residuals,
753  }
754 
755 
756  /// \short Add the element's contribution to its residual vector and
757  /// its Jacobian matrix
759  DenseMatrix<double> &jacobian)
760  {
761  //Call the generic routine with the flag set to 1
763  jacobian,
764  1);
765  }
766 
767  /// Specify the value of nodal zeta from the face geometry
768  /// \short The "global" intrinsic coordinate of the element when
769  /// viewed as part of a geometric object should be given by
770  /// the FaceElement representation, by default (needed to break
771  /// indeterminacy if bulk element is SolidElement)
772  double zeta_nodal(const unsigned &n, const unsigned &k,
773  const unsigned &i) const
774  {return FaceElement::zeta_nodal(n,k,i);}
775 
776  /// \short Output function -- forward to broken version in FiniteElement
777  /// until somebody decides what exactly they want to plot here...
778  void output(std::ostream &outfile)
779  {
780  FiniteElement::output(outfile);
781  }
782 
783  /// \short Output function -- forward to broken version in FiniteElement
784  /// until somebody decides what exactly they want to plot here...
785  void output(std::ostream &outfile, const unsigned &nplot)
786  {
787  FiniteElement::output(outfile,nplot);
788  }
789 
790 
791 protected:
792 
793  /// \short Function to compute the shape and test functions and to return
794  /// the Jacobian of mapping between local and global (Eulerian)
795  /// coordinates
796  inline double shape_and_test(const Vector<double> &s,
797  Shape &psi,
798  Shape &test) const
799  {
800  //Find number of nodes
801  unsigned n_node = nnode();
802 
803  //Get the shape functions
804  shape(s,psi);
805 
806  //Set the test functions to be the same as the shape functions
807  for(unsigned i=0;i<n_node;i++)
808  {
809  test[i] = psi[i];
810  }
811 
812  //Return the value of the jacobian
813  return J_eulerian(s);
814  }
815 
816 
817  /// \short Function to compute the shape and test functions and to return
818  /// the Jacobian of mapping between local and global (Eulerian)
819  /// coordinates
820  inline double shape_and_test_at_knot(const unsigned &ipt,
821  Shape &psi,
822  Shape &test) const
823  {
824  //Find number of nodes
825  unsigned n_node = nnode();
826 
827  //Get the shape functions
828  shape_at_knot(ipt,psi);
829 
830  //Set the test functions to be the same as the shape functions
831  for(unsigned i=0;i<n_node;i++)
832  {
833  test[i] = psi[i];
834  }
835 
836  //Return the value of the jacobian
837  return J_eulerian_at_knot(ipt);
838  }
839 
840  /// \short Function to calculate the prescribed beta at a given spatial
841  /// position
842  void get_beta(const Vector<double>& x, double& beta)
843  {
844  //If the function pointer is zero return zero
845  if(Beta_fct_pt == 0)
846  {
847  beta = 0.0;
848  }
849  //Otherwise call the function
850  else
851  {
852  (*Beta_fct_pt)(x,beta);
853  }
854  }
855 
856  /// \short Function to calculate the prescribed alpha at a given spatial
857  /// position
858  void get_alpha(const Vector<double>& x, double& alpha)
859  {
860  //If the function pointer is zero return zero
861  if(Alpha_fct_pt == 0)
862  {
863  alpha = 0.0;
864  }
865  //Otherwise call the function
866  else
867  {
868  (*Alpha_fct_pt)(x,alpha);
869  }
870  }
871 
872 private:
873 
874 
875  /// \short Add the element's contribution to its residual vector.
876  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
878  Vector<double> &residuals, DenseMatrix<double> &jacobian,
879  unsigned flag);
880 
881 
882  /// Function pointer to the (global) prescribed-beta function
884 
885  /// Function pointer to the (global) prescribed-alpha function
887 
888  /// The index at which the unknown is stored at the nodes
890 
891 
892 };//End class SphericalAdvectionDiffusionFluxElement
893 
894 
895 ///////////////////////////////////////////////////////////////////////
896 ///////////////////////////////////////////////////////////////////////
897 ///////////////////////////////////////////////////////////////////////
898 
899 
900 //===========================================================================
901 /// \short Constructor, takes the pointer to the "bulk" element and the index
902 /// of the face to be created
903 //===========================================================================
904 template<class ELEMENT>
907  const int &face_index) :
908 FaceGeometry<ELEMENT>(), FaceElement()
909 
910 {
911  //Let the bulk element build the FaceElement, i.e. setup the pointers
912  //to its nodes (by referring to the appropriate nodes in the bulk
913  //element), etc.
914  bulk_el_pt->build_face_element(face_index,this);
915 
916 #ifdef PARANOID
917  {
918  //Check that the element is not a refineable 3d element
919  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
920 
921  //If it's three-d
922  if(elem_pt->dim()==3)
923  {
924  //Is it refineable
925  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
926  if(ref_el_pt!=0)
927  {
928  if (this->has_hanging_nodes())
929  {
930  throw OomphLibError(
931  "This flux element will not work correctly if nodes are hanging\n",
932  OOMPH_CURRENT_FUNCTION,
933  OOMPH_EXCEPTION_LOCATION);
934  }
935  }
936  }
937  }
938 #endif
939 
940  //Initialise the prescribed-beta function pointer to zero
941  Beta_fct_pt = 0;
942 
943  //Set up U_index_adv_diff. Initialise to zero, which probably won't change
944  //in most cases, oh well, the price we pay for generality
945  U_index_adv_diff = 0;
946 
947  //Cast to the appropriate AdvectionDiffusionEquation so that we can
948  //find the index at which the variable is stored
949  //We assume that the dimension of the full problem is the same
950  //as the dimension of the node, if this is not the case you will have
951  //to write custom elements, sorry
953  dynamic_cast<SphericalAdvectionDiffusionEquations*>(bulk_el_pt);
954 
955  //If the cast has failed die
956  if(eqn_pt==0)
957  {
958  std::string error_string =
959  "Bulk element must inherit from SphericalAdvectionDiffusionEquations.";
960  error_string +=
961  "Nodes are two dimensional, but cannot cast the bulk element to\n";
962  error_string += "SphericalAdvectionDiffusionEquations<2>\n.";
963  error_string +=
964  "If you desire this functionality, you must implement it yourself\n";
965 
966  throw OomphLibError(
967  error_string,
968  OOMPH_CURRENT_FUNCTION,
969  OOMPH_EXCEPTION_LOCATION);
970  }
971  else
972  {
973  //Read the index from the (cast) bulk element.
975  }
976 
977 
978 }
979 
980 
981 //===========================================================================
982 /// \short Compute the element's residual vector and the (zero) Jacobian
983 /// matrix for the Robin boundary condition:
984 /// \f[
985 /// \Delta u \cdot \mathbf{n} + \alpha (\mathbf{x}) = \beta (\mathbf{x})
986 /// \f]
987 //===========================================================================
988 template<class ELEMENT>
991  Vector<double> &residuals,
992  DenseMatrix<double> &jacobian,
993  unsigned flag)
994 {
995  //Find out how many nodes there are
996  const unsigned n_node = nnode();
997 
998  // Locally cache the index at which the variable is stored
999  const unsigned u_index_spherical_adv_diff = U_index_adv_diff;
1000 
1001  //Set up memory for the shape and test functions
1002  Shape psif(n_node), testf(n_node);
1003 
1004  //Set the value of n_intpt
1005  const unsigned n_intpt = integral_pt()->nweight();
1006 
1007  //Set the Vector to hold local coordinates
1008  Vector<double> s(1);
1009 
1010  //Integers used to store the local equation number and local unknown
1011  //indices for the residuals and jacobians
1012  int local_eqn=0, local_unknown=0;
1013 
1014  //Loop over the integration points
1015  //--------------------------------
1016  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1017  {
1018 
1019  //Assign values of s
1020  for(unsigned i=0;i<1;i++)
1021  {
1022  s[i] = integral_pt()->knot(ipt,i);
1023  }
1024 
1025  //Get the integral weight
1026  double w = integral_pt()->weight(ipt);
1027 
1028  //Find the shape and test functions and return the Jacobian
1029  //of the mapping
1030  double J = shape_and_test(s,psif,testf);
1031 
1032  //Premultiply the weights and the Jacobian
1033  double W = w*J;
1034 
1035  //Calculate local values of the solution and its derivatives
1036  //Allocate
1037  double interpolated_u=0.0;
1038  Vector<double> interpolated_x(2,0.0);
1039 
1040  //Calculate position
1041  for(unsigned l=0;l<n_node;l++)
1042  {
1043  //Get the value at the node
1044  double u_value = raw_nodal_value(l,u_index_spherical_adv_diff);
1045  interpolated_u += u_value*psif(l);
1046  //Loop over coordinate direction
1047  for(unsigned i=0;i<2;i++)
1048  {
1049  interpolated_x[i] += nodal_position(l,i)*psif(l);
1050  }
1051  }
1052 
1053  //Get the imposed beta (beta=flux when alpha=0.0)
1054  double beta;
1055  get_beta(interpolated_x,beta);
1056 
1057  //Get the imposed alpha
1058  double alpha;
1059  get_alpha(interpolated_x,alpha);
1060 
1061  //calculate the area weighting dS = r^{2} sin theta dr dtheta
1062  // r = x[0] and theta = x[1]
1063  double dS = interpolated_x[0]*interpolated_x[0]*sin(interpolated_x[1]);
1064 
1065  //Now add to the appropriate equations
1066 
1067  //Loop over the test functions
1068  for(unsigned l=0;l<n_node;l++)
1069  {
1070  //Set the local equation number
1071  local_eqn = nodal_local_eqn(l,u_index_spherical_adv_diff);
1072  /*IF it's not a boundary condition*/
1073  if(local_eqn >= 0)
1074  {
1075  //Add the prescribed beta terms
1076  residuals[local_eqn] -= dS*(beta - alpha*interpolated_u)*testf(l)*W;
1077 
1078  //Calculate the Jacobian
1079  //----------------------
1080  if ( (flag) && (alpha!=0.0) )
1081  {
1082  //Loop over the velocity shape functions again
1083  for(unsigned l2=0;l2<n_node;l2++)
1084  {
1085  //Set the number of the unknown
1086  local_unknown = nodal_local_eqn(l2,u_index_spherical_adv_diff);
1087 
1088  //If at a non-zero degree of freedom add in the entry
1089  if(local_unknown >= 0)
1090  {
1091  jacobian(local_eqn,local_unknown) +=
1092  dS*alpha*psif[l2]*testf[l]*W;
1093  }
1094  }
1095  }
1096 
1097  }
1098  } //end loop over test functions
1099 
1100  } //end loop over integration points
1101 
1102 }//end fill_in_generic_residual_contribution_adv_diff_flux
1103 
1104 
1105 } //End SphericalAdvectionDiffusionFluxElement class
1106 
1107 #endif
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
QSphericalAdvectionDiffusionElement(const QSphericalAdvectionDiffusionElement< NNODE_1D > &dummy)
Broken copy constructor.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
QSphericalAdvectionDiffusionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
virtual void dinterpolated_u_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_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
void get_alpha(const Vector< double > &x, double &alpha)
Function to calculate the prescribed alpha at a given spatial position.
void(* SphericalAdvectionDiffusionPrescribedBetaFctPt)(const Vector< double > &x, double &beta)
Function pointer to the prescribed-beta function fct(x,beta(x)) – x is a Vector!
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
void output(std::ostream &outfile)
Output with default number of plot points.
const double & pe_st() const
Peclet number multiplied by Strouhal number.
void(* SphericalAdvectionDiffusionPrescribedAlphaFctPt)(const Vector< double > &x, double &alpha)
Function pointer to the prescribed-alpha function fct(x,alpha(x)) – x is a Vector! ...
SphericalAdvectionDiffusionPrescribedAlphaFctPt Alpha_fct_pt
Function pointer to the (global) prescribed-alpha function.
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
SphericalAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
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
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.
double interpolated_u_spherical_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
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) ...
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
void operator=(const SphericalAdvectionDiffusionEquations &)
Broken assignment operator.
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 void fill_in_generic_residual_contribution_spherical_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...
A general Finite Element class.
Definition: elements.h:1271
unsigned required_nvalue(const unsigned &n) const
Required # of `values' (pinned or dofs) at node n.
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(* SphericalAdvectionDiffusionWindFctPt)(const Vector< double > &x, Vector< double > &wind)
Function pointer to wind function fct(x,w(x)) – x is a Vector!
SphericalAdvectionDiffusionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
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
static double Default_peclet_number
Static default value for the Peclet number.
void output(std::ostream &outfile)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
double dshape_and_dtest_eulerian_spherical_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.
SphericalAdvectionDiffusionPrescribedBetaFctPt Beta_fct_pt
Function pointer to the (global) prescribed-beta function.
A class for all elements that solve the Advection Diffusion equations in a spherical polar coordinate...
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Add the element's contribution to its residual vector and the element Jacobian matrix (wrapper) and m...
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
SphericalAdvectionDiffusionEquations()
Constructor: Initialise the Source_fct_pt and Wind_fct_pt to null and set (pointer to) Peclet number ...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
double du_dt_spherical_adv_diff(const unsigned &n) const
du/dt at local node n.
double dshape_and_dtest_eulerian_at_knot_spherical_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.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
static char t char * s
Definition: cfortran.h:572
virtual void get_source_spherical_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(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.
SphericalAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
virtual double dshape_and_dtest_eulerian_at_knot_spherical_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 ...
A class for elements that allow the imposition of an applied Robin boundary condition on the boundari...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
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 output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,z,u at n_plot^2 plot points.
unsigned U_index_adv_diff
The index at which the unknown is stored at the nodes.
void operator=(const SphericalAdvectionDiffusionFluxElement &)
Broken assignment operator.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
void output(std::ostream &outfile, const unsigned &nplot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
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
QSphericalAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Axisymmetric Advec...
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
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
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.
void output(std::ostream &outfile)
Output function: r,z,u.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
SphericalAdvectionDiffusionPrescribedBetaFctPt & beta_fct_pt()
Access function for the prescribed-beta function pointer.
SphericalAdvectionDiffusionEquations(const SphericalAdvectionDiffusionEquations &dummy)
Broken copy constructor.
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.
void fill_in_generic_residual_contribution_spherical_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...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void(* SphericalAdvectionDiffusionSourceFctPt)(const Vector< double > &x, double &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
SphericalAdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2344
SphericalAdvectionDiffusionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
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
void get_beta(const Vector< double > &x, double &beta)
Function to calculate the prescribed beta at a given spatial position.
virtual void get_wind_spherical_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...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,z,u at n_plot^2 plot points.
void operator=(const QSphericalAdvectionDiffusionElement< NNODE_1D > &)
Broken assignment operator.
SphericalAdvectionDiffusionFluxElement(const SphericalAdvectionDiffusionFluxElement &dummy)
Broken copy constructor.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
virtual unsigned u_index_spherical_adv_diff() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
SphericalAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
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 ...
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 ...
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...
SphericalAdvectionDiffusionPrescribedAlphaFctPt & alpha_fct_pt()
Access function for the prescribed-alpha function pointer.
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 disable_ALE()
Disable ALE – empty overload to suppress warning. ALE isn't implemented anyway.
virtual double dshape_and_dtest_eulerian_spherical_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...
void output(FILE *file_pt)
C-style output function: r,z,u.