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// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Header file for advection diffusion elements in a spherical polar coordinate
27 // system
28 #ifndef OOMPH_SPHERICAL_ADV_DIFF_ELEMENTS_HEADER
29 #define OOMPH_SPHERICAL_ADV_DIFF_ELEMENTS_HEADER
30 
31 
32 // Config header generated by autoconfig
33 #ifdef HAVE_CONFIG_H
34 #include <oomph-lib-config.h>
35 #endif
36 
37 // OOMPH-LIB headers
38 #include "../generic/nodes.h"
39 #include "../generic/Qelements.h"
40 #include "../generic/refineable_elements.h"
41 #include "../generic/oomph_utilities.h"
42 
43 namespace oomph
44 {
45  //=============================================================
46  /// A class for all elements that solve the
47  /// Advection Diffusion equations in a spherical polar coordinate system
48  /// using isoparametric elements.
49  /// \f[ Pe \mathbf{w}\cdot(\mathbf{x}) \nabla u = \nabla \cdot \left( \nabla u \right) + f(\mathbf{x}) \f]
50  /// This contains the generic maths. Shape functions, geometric
51  /// mapping etc. must get implemented in derived class.
52  //=============================================================
54  {
55  public:
56  /// Function pointer to source function fct(x,f(x)) --
57  /// x is a Vector!
59  const Vector<double>& x, double& f);
60 
61 
62  /// Function pointer to wind function fct(x,w(x)) --
63  /// x is a Vector!
65  const Vector<double>& x, Vector<double>& wind);
66 
67 
68  /// Constructor: Initialise the Source_fct_pt and Wind_fct_pt
69  /// to null and set (pointer to) Peclet number to default
71  {
72  // Set pointer to Peclet number to the default value zero
75  }
76 
77  /// Broken copy constructor
79  const SphericalAdvectionDiffusionEquations& dummy) = delete;
80 
81  /// Broken assignment operator
83 
84  /// Return the index at which the unknown value
85  /// is stored. The default value, 0, is appropriate for single-physics
86  /// problems, when there is only one variable, the value that satisfies
87  /// the spherical advection-diffusion equation.
88  /// In derived multi-physics elements, this function should be overloaded
89  /// to reflect the chosen storage scheme. Note that these equations require
90  /// that the unknown is always stored at the same index at each node.
91  virtual inline unsigned u_index_spherical_adv_diff() const
92  {
93  return 0;
94  }
95 
96 
97  /// du/dt at local node n.
98  /// Uses suitably interpolated value for hanging nodes.
99  double du_dt_spherical_adv_diff(const unsigned& n) const
100  {
101  // Get the data's timestepper
103 
104  // Initialise dudt
105  double dudt = 0.0;
106  // Loop over the timesteps, if there is a non Steady timestepper
107  if (!time_stepper_pt->is_steady())
108  {
109  // Find the index at which the variable is stored
110  const unsigned u_nodal_index = u_index_spherical_adv_diff();
111 
112  // Number of timsteps (past & present)
113  const unsigned n_time = time_stepper_pt->ntstorage();
114 
115  for (unsigned t = 0; t < n_time; t++)
116  {
117  dudt +=
118  time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
119  }
120  }
121  return dudt;
122  }
123 
124 
125  /// Disable ALE -- empty overload to suppress warning.
126  /// ALE isn't implemented anyway
127  void disable_ALE() {}
128 
129  /// Output with default number of plot points
130  void output(std::ostream& outfile)
131  {
132  unsigned nplot = 5;
133  output(outfile, nplot);
134  }
135 
136  /// Output FE representation of soln: r,z,u at
137  /// nplot^2 plot points
138  void output(std::ostream& outfile, const unsigned& nplot);
139 
140 
141  /// C_style output with default number of plot points
142  void output(FILE* file_pt)
143  {
144  unsigned n_plot = 5;
145  output(file_pt, n_plot);
146  }
147 
148  /// C-style output FE representation of soln: r,z,u at
149  /// n_plot^2 plot points
150  void output(FILE* file_pt, const unsigned& n_plot);
151 
152 
153  /// Output exact soln: r,z,u_exact at nplot^2 plot points
154  void output_fct(std::ostream& outfile,
155  const unsigned& nplot,
157 
158  /// Get error against and norm of exact solution
159  void compute_error(std::ostream& outfile,
161  double& error,
162  double& norm);
163 
164  /// Access function: Pointer to source function
166  {
167  return Source_fct_pt;
168  }
169 
170 
171  /// Access function: Pointer to source function. Const version
173  {
174  return Source_fct_pt;
175  }
176 
177 
178  /// Access function: Pointer to wind function
180  {
181  return Wind_fct_pt;
182  }
183 
184 
185  /// Access function: Pointer to wind function. Const version
187  {
188  return Wind_fct_pt;
189  }
190 
191  // Access functions for the physical constants
192 
193  /// Peclet number
194  inline const double& pe() const
195  {
196  return *Pe_pt;
197  }
198 
199  /// Pointer to Peclet number
200  inline double*& pe_pt()
201  {
202  return Pe_pt;
203  }
204 
205  /// Peclet number multiplied by Strouhal number
206  inline const double& pe_st() const
207  {
208  return *PeSt_pt;
209  }
210 
211  /// Pointer to Peclet number multipled by Strouha number
212  inline double*& pe_st_pt()
213  {
214  return PeSt_pt;
215  }
216 
217  /// Get source term at (Eulerian) position x. This function is
218  /// virtual to allow overloading in multi-physics problems where
219  /// the strength of the source function might be determined by
220  /// another system of equations
221  inline virtual void get_source_spherical_adv_diff(const unsigned& ipt,
222  const Vector<double>& x,
223  double& source) const
224  {
225  // If no source function has been set, return zero
226  if (Source_fct_pt == 0)
227  {
228  source = 0.0;
229  }
230  else
231  {
232  // Get source strength
233  (*Source_fct_pt)(x, source);
234  }
235  }
236 
237  /// Get wind at (Eulerian) position x and/or local coordinate s.
238  /// This function is
239  /// virtual to allow overloading in multi-physics problems where
240  /// the wind function might be determined by
241  /// another system of equations
242  inline virtual void get_wind_spherical_adv_diff(const unsigned& ipt,
243  const Vector<double>& s,
244  const Vector<double>& x,
245  Vector<double>& wind) const
246  {
247  // If no wind function has been set, return zero
248  if (Wind_fct_pt == 0)
249  {
250  for (unsigned i = 0; i < 3; i++)
251  {
252  wind[i] = 0.0;
253  }
254  }
255  else
256  {
257  // Get wind
258  (*Wind_fct_pt)(x, wind);
259  }
260  }
261 
262  /// Get flux:
263  /// \f[ \mbox{flux}[i] = \nabla u = \mbox{d}u / \mbox{d} r + 1/r \mbox{d}u / \mbox{d} \theta \f]
264  void get_flux(const Vector<double>& s, Vector<double>& flux) const
265  {
266  // Find out how many nodes there are in the element
267  const unsigned n_node = nnode();
268 
269  // Get the nodal index at which the unknown is stored
270  const unsigned u_nodal_index = u_index_spherical_adv_diff();
271 
272  // Set up memory for the shape and test functions
273  Shape psi(n_node);
274  DShape dpsidx(n_node, 2);
275 
276  // Call the derivatives of the shape and test functions
277  dshape_eulerian(s, psi, dpsidx);
278 
279  // Initialise to zero
280  for (unsigned j = 0; j < 2; j++)
281  {
282  flux[j] = 0.0;
283  }
284 
285  // Loop over nodes
286  for (unsigned l = 0; l < n_node; l++)
287  {
288  const double u_value = this->nodal_value(l, u_nodal_index);
289  const double r = this->nodal_position(l, 0);
290  // Add in the derivative directions
291  flux[0] += u_value * dpsidx(l, 0);
292  flux[1] += u_value * dpsidx(l, 1) / r;
293  }
294  }
295 
296 
297  /// Add the element's contribution to its residual vector (wrapper)
299  {
300  // Call the generic residuals function with flag set to 0 and using
301  // a dummy matrix
303  residuals,
306  0);
307  }
308 
309 
310  /// Add the element's contribution to its residual vector and
311  /// the element Jacobian matrix (wrapper)
313  DenseMatrix<double>& jacobian)
314  {
315  // Call the generic routine with the flag set to 1
317  residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
318  }
319 
320  /// Add the element's contribution to its residual vector and
321  /// the element Jacobian matrix (wrapper) and mass matrix
323  Vector<double>& residuals,
324  DenseMatrix<double>& jacobian,
325  DenseMatrix<double>& mass_matrix)
326  {
327  // Call the generic routine with the flag set to 2
329  residuals, jacobian, mass_matrix, 2);
330  }
331 
332 
333  /// Return FE representation of function value u(s) at local coordinate s
335  const Vector<double>& s) const
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  /// 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  /// Self-test: Return 0 for OK
421  unsigned self_test();
422 
423  protected:
424  /// Shape/test functions and derivs w.r.t. to global coords at
425  /// local coord. s; return Jacobian of mapping
427  const Vector<double>& s,
428  Shape& psi,
429  DShape& dpsidx,
430  Shape& test,
431  DShape& dtestdx) const = 0;
432 
433  /// Shape/test functions and derivs w.r.t. to global coords at
434  /// integration point ipt; return Jacobian of mapping
436  const unsigned& ipt,
437  Shape& psi,
438  DShape& dpsidx,
439  Shape& test,
440  DShape& dtestdx) const = 0;
441 
442  /// Add the element's contribution to its residual vector only
443  /// (if flag=and/or element Jacobian matrix
445  Vector<double>& residuals,
446  DenseMatrix<double>& jacobian,
447  DenseMatrix<double>& mass_matrix,
448  unsigned flag);
449 
450  // Physical constants
451 
452  /// Pointer to global Peclet number
453  double* Pe_pt;
454 
455  /// Pointer to global Peclet number multiplied by Strouhal number
456  double* PeSt_pt;
457 
458  /// Pointer to source function:
460 
461  /// Pointer to wind function:
463 
464  private:
465  /// Static default value for the Peclet number
466  static double Default_peclet_number;
467 
468 
469  }; // End class SphericalAdvectionDiffusionEquations
470 
471 
472  /// ////////////////////////////////////////////////////////////////////////
473  /// ////////////////////////////////////////////////////////////////////////
474  /// ////////////////////////////////////////////////////////////////////////
475 
476 
477  //======================================================================
478  /// QSphericalAdvectionDiffusionElement elements are
479  /// linear/quadrilateral/brick-shaped Axisymmetric Advection Diffusion
480  /// elements with isoparametric interpolation for the function.
481  //======================================================================
482  template<unsigned NNODE_1D>
484  : public virtual QElement<2, NNODE_1D>,
486  {
487  private:
488  /// Static array of ints to hold number of variables at
489  /// nodes: Initial_Nvalue[n]
490  static const unsigned Initial_Nvalue;
491 
492  public:
493  /// Constructor: Call constructors for QElement and
494  /// Advection Diffusion equations
497  {
498  }
499 
500  /// Broken copy constructor
502  const QSphericalAdvectionDiffusionElement<NNODE_1D>& dummy) = delete;
503 
504  /// Broken assignment operator
506  delete;
507 
508  /// Required # of `values' (pinned or dofs)
509  /// at node n
510  inline unsigned required_nvalue(const unsigned& n) const
511  {
512  return Initial_Nvalue;
513  }
514 
515  /// Output function:
516  /// r,z,u
517  void output(std::ostream& outfile)
518  {
520  }
521 
522  /// Output function:
523  /// r,z,u at n_plot^2 plot points
524  void output(std::ostream& outfile, const unsigned& n_plot)
525  {
527  }
528 
529 
530  /// C-style output function:
531  /// r,z,u
532  void output(FILE* file_pt)
533  {
535  }
536 
537  /// C-style output function:
538  /// r,z,u at n_plot^2 plot points
539  void output(FILE* file_pt, const unsigned& n_plot)
540  {
542  }
543 
544  /// Output function for an exact solution:
545  /// r,z,u_exact at n_plot^2 plot points
546  void output_fct(std::ostream& outfile,
547  const unsigned& n_plot,
549  {
551  outfile, n_plot, exact_soln_pt);
552  }
553 
554 
555  protected:
556  /// Shape, test functions & derivs. w.r.t. to global coords. Return
557  /// Jacobian.
559  const Vector<double>& s,
560  Shape& psi,
561  DShape& dpsidx,
562  Shape& test,
563  DShape& dtestdx) const;
564 
565  /// Shape, test functions & derivs. w.r.t. to global coords. at
566  /// integration point ipt. Return Jacobian.
568  const unsigned& ipt,
569  Shape& psi,
570  DShape& dpsidx,
571  Shape& test,
572  DShape& dtestdx) const;
573 
574  }; // End class QSphericalAdvectionDiffusionElement
575 
576  // Inline functions:
577 
578  //======================================================================
579  /// Define the shape functions and test functions and derivatives
580  /// w.r.t. global coordinates and return Jacobian of mapping.
581  ///
582  /// Galerkin: Test functions = shape functions
583  //======================================================================
584 
585  template<unsigned NNODE_1D>
588  Shape& psi,
589  DShape& dpsidx,
590  Shape& test,
591  DShape& dtestdx) const
592  {
593  // Call the geometrical shape functions and derivatives
594  double J = this->dshape_eulerian(s, psi, dpsidx);
595 
596  // Loop over the test functions and derivatives and set them equal to the
597  // shape functions
598  for (unsigned i = 0; i < NNODE_1D; i++)
599  {
600  test[i] = psi[i];
601  for (unsigned j = 0; j < 2; j++)
602  {
603  dtestdx(i, j) = dpsidx(i, j);
604  }
605  }
606 
607  // Return the jacobian
608  return J;
609  }
610 
611 
612  //======================================================================
613  /// Define the shape functions and test functions and derivatives
614  /// w.r.t. global coordinates and return Jacobian of mapping.
615  ///
616  /// Galerkin: Test functions = shape functions
617  //======================================================================
618 
619  template<unsigned NNODE_1D>
622  Shape& psi,
623  DShape& dpsidx,
624  Shape& test,
625  DShape& dtestdx) const
626  {
627  // Call the geometrical shape functions and derivatives
628  double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
629 
630  // Set the test functions equal to the shape functions (pointer copy)
631  test = psi;
632  dtestdx = dpsidx;
633 
634  // Return the jacobian
635  return J;
636  }
637 
638  /// /////////////////////////////////////////////////////////////////////
639  /// /////////////////////////////////////////////////////////////////////
640  /// /////////////////////////////////////////////////////////////////////
641 
642  template<unsigned NNODE_1D>
644  : public virtual QElement<1, NNODE_1D>
645  {
646  public:
647  /// Constructor: Call the constructor for the
648  /// appropriate lower-dimensional QElement
649  FaceGeometry() : QElement<1, NNODE_1D>() {}
650  };
651 
652 
653  /// /////////////////////////////////////////////////////////////////////
654  /// /////////////////////////////////////////////////////////////////////
655  /// /////////////////////////////////////////////////////////////////////
656 
657  //======================================================================
658  /// A class for elements that allow the imposition of an
659  /// applied Robin boundary condition on the boundaries of Steady
660  /// Axisymmnetric Advection Diffusion Flux elements.
661  /// \f[ -\Delta u \cdot \mathbf{n} + \alpha(r,z) u = \beta(r,z) \f]
662  /// The element geometry is obtained from the FaceGeometry<ELEMENT>
663  /// policy class.
664  //======================================================================
665  template<class ELEMENT>
667  : public virtual FaceGeometry<ELEMENT>,
668  public virtual FaceElement
669  {
670  public:
671  /// Function pointer to the prescribed-beta function fct(x,beta(x))
672  /// -- x is a Vector!
674  const Vector<double>& x, double& beta);
675 
676  /// Function pointer to the prescribed-alpha function fct(x,alpha(x))
677  /// -- x is a Vector!
679  const Vector<double>& x, double& alpha);
680 
681 
682  /// Constructor, takes the pointer to the "bulk" element
683  /// and the index of the face to be created
685  const int& face_index);
686 
687 
688  /// Broken empty constructor
690  {
691  throw OomphLibError("Don't call empty constructor for "
692  "SphericalAdvectionDiffusionFluxElement",
693  OOMPH_CURRENT_FUNCTION,
694  OOMPH_EXCEPTION_LOCATION);
695  }
696 
697  /// Broken copy constructor
699  const SphericalAdvectionDiffusionFluxElement& dummy) = delete;
700 
701  /// Broken assignment operator
703 
704  /// Access function for the prescribed-beta function pointer
706  {
707  return Beta_fct_pt;
708  }
709 
710  /// Access function for the prescribed-alpha function pointer
712  {
713  return Alpha_fct_pt;
714  }
715 
716 
717  /// Add the element's contribution to its residual vector
719  {
720  // Call the generic residuals function with flag set to 0
721  // using a dummy matrix
723  residuals, GeneralisedElement::Dummy_matrix, 0);
724  }
725 
726 
727  /// Add the element's contribution to its residual vector and
728  /// its Jacobian matrix
730  DenseMatrix<double>& jacobian)
731  {
732  // Call the generic routine with the flag set to 1
734  residuals, jacobian, 1);
735  }
736 
737  /// Specify the value of nodal zeta from the face geometry
738  /// The "global" intrinsic coordinate of the element when
739  /// viewed as part of a geometric object should be given by
740  /// the FaceElement representation, by default (needed to break
741  /// indeterminacy if bulk element is SolidElement)
742  double zeta_nodal(const unsigned& n,
743  const unsigned& k,
744  const unsigned& i) const
745  {
746  return FaceElement::zeta_nodal(n, k, i);
747  }
748 
749  /// Output function -- forward to broken version in FiniteElement
750  /// until somebody decides what exactly they want to plot here...
751  void output(std::ostream& outfile)
752  {
753  FiniteElement::output(outfile);
754  }
755 
756  /// Output function -- forward to broken version in FiniteElement
757  /// until somebody decides what exactly they want to plot here...
758  void output(std::ostream& outfile, const unsigned& nplot)
759  {
760  FiniteElement::output(outfile, nplot);
761  }
762 
763 
764  protected:
765  /// Function to compute the shape and test functions and to return
766  /// the Jacobian of mapping between local and global (Eulerian)
767  /// coordinates
768  inline double shape_and_test(const Vector<double>& s,
769  Shape& psi,
770  Shape& test) const
771  {
772  // Find number of nodes
773  unsigned n_node = nnode();
774 
775  // Get the shape functions
776  shape(s, psi);
777 
778  // Set the test functions to be the same as the shape functions
779  for (unsigned i = 0; i < n_node; i++)
780  {
781  test[i] = psi[i];
782  }
783 
784  // Return the value of the jacobian
785  return J_eulerian(s);
786  }
787 
788 
789  /// Function to compute the shape and test functions and to return
790  /// the Jacobian of mapping between local and global (Eulerian)
791  /// coordinates
792  inline double shape_and_test_at_knot(const unsigned& ipt,
793  Shape& psi,
794  Shape& test) const
795  {
796  // Find number of nodes
797  unsigned n_node = nnode();
798 
799  // Get the shape functions
800  shape_at_knot(ipt, psi);
801 
802  // Set the test functions to be the same as the shape functions
803  for (unsigned i = 0; i < n_node; i++)
804  {
805  test[i] = psi[i];
806  }
807 
808  // Return the value of the jacobian
809  return J_eulerian_at_knot(ipt);
810  }
811 
812  /// Function to calculate the prescribed beta at a given spatial
813  /// position
814  void get_beta(const Vector<double>& x, double& beta)
815  {
816  // If the function pointer is zero return zero
817  if (Beta_fct_pt == 0)
818  {
819  beta = 0.0;
820  }
821  // Otherwise call the function
822  else
823  {
824  (*Beta_fct_pt)(x, beta);
825  }
826  }
827 
828  /// Function to calculate the prescribed alpha at a given spatial
829  /// position
830  void get_alpha(const Vector<double>& x, double& alpha)
831  {
832  // If the function pointer is zero return zero
833  if (Alpha_fct_pt == 0)
834  {
835  alpha = 0.0;
836  }
837  // Otherwise call the function
838  else
839  {
840  (*Alpha_fct_pt)(x, alpha);
841  }
842  }
843 
844  private:
845  /// Add the element's contribution to its residual vector.
846  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
848  Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
849 
850 
851  /// Function pointer to the (global) prescribed-beta function
853 
854  /// Function pointer to the (global) prescribed-alpha function
856 
857  /// The index at which the unknown is stored at the nodes
859 
860 
861  }; // End class SphericalAdvectionDiffusionFluxElement
862 
863 
864  /// ////////////////////////////////////////////////////////////////////
865  /// ////////////////////////////////////////////////////////////////////
866  /// ////////////////////////////////////////////////////////////////////
867 
868 
869  //===========================================================================
870  /// Constructor, takes the pointer to the "bulk" element and the index
871  /// of the face to be created
872  //===========================================================================
873  template<class ELEMENT>
876  const int& face_index)
877  : FaceGeometry<ELEMENT>(), FaceElement()
878 
879  {
880  // Let the bulk element build the FaceElement, i.e. setup the pointers
881  // to its nodes (by referring to the appropriate nodes in the bulk
882  // element), etc.
883  bulk_el_pt->build_face_element(face_index, this);
884 
885 #ifdef PARANOID
886  {
887  // Check that the element is not a refineable 3d element
888  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
889 
890  // If it's three-d
891  if (elem_pt->dim() == 3)
892  {
893  // Is it refineable
894  RefineableElement* ref_el_pt =
895  dynamic_cast<RefineableElement*>(elem_pt);
896  if (ref_el_pt != 0)
897  {
898  if (this->has_hanging_nodes())
899  {
900  throw OomphLibError("This flux element will not work correctly if "
901  "nodes are hanging\n",
902  OOMPH_CURRENT_FUNCTION,
903  OOMPH_EXCEPTION_LOCATION);
904  }
905  }
906  }
907  }
908 #endif
909 
910  // Initialise the prescribed-beta function pointer to zero
911  Beta_fct_pt = 0;
912 
913  // Set up U_index_adv_diff. Initialise to zero, which probably won't change
914  // in most cases, oh well, the price we pay for generality
915  U_index_adv_diff = 0;
916 
917  // Cast to the appropriate AdvectionDiffusionEquation so that we can
918  // find the index at which the variable is stored
919  // We assume that the dimension of the full problem is the same
920  // as the dimension of the node, if this is not the case you will have
921  // to write custom elements, sorry
923  dynamic_cast<SphericalAdvectionDiffusionEquations*>(bulk_el_pt);
924 
925  // If the cast has failed die
926  if (eqn_pt == 0)
927  {
928  std::string error_string =
929  "Bulk element must inherit from SphericalAdvectionDiffusionEquations.";
930  error_string +=
931  "Nodes are two dimensional, but cannot cast the bulk element to\n";
932  error_string += "SphericalAdvectionDiffusionEquations<2>\n.";
933  error_string +=
934  "If you desire this functionality, you must implement it yourself\n";
935 
936  throw OomphLibError(
937  error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
938  }
939  else
940  {
941  // Read the index from the (cast) bulk element.
943  }
944  }
945 
946 
947  //===========================================================================
948  /// Compute the element's residual vector and the (zero) Jacobian
949  /// matrix for the Robin boundary condition:
950  /// \f[ \Delta u \cdot \mathbf{n} + \alpha (\mathbf{x}) = \beta (\mathbf{x}) \f]
951  //===========================================================================
952  template<class ELEMENT>
955  Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
956  {
957  // Find out how many nodes there are
958  const unsigned n_node = nnode();
959 
960  // Locally cache the index at which the variable is stored
961  const unsigned u_index_spherical_adv_diff = U_index_adv_diff;
962 
963  // Set up memory for the shape and test functions
964  Shape psif(n_node), testf(n_node);
965 
966  // Set the value of n_intpt
967  const unsigned n_intpt = integral_pt()->nweight();
968 
969  // Set the Vector to hold local coordinates
970  Vector<double> s(1);
971 
972  // Integers used to store the local equation number and local unknown
973  // indices for the residuals and jacobians
974  int local_eqn = 0, local_unknown = 0;
975 
976  // Loop over the integration points
977  //--------------------------------
978  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
979  {
980  // Assign values of s
981  for (unsigned i = 0; i < 1; i++)
982  {
983  s[i] = integral_pt()->knot(ipt, i);
984  }
985 
986  // Get the integral weight
987  double w = integral_pt()->weight(ipt);
988 
989  // Find the shape and test functions and return the Jacobian
990  // of the mapping
991  double J = shape_and_test(s, psif, testf);
992 
993  // Premultiply the weights and the Jacobian
994  double W = w * J;
995 
996  // Calculate local values of the solution and its derivatives
997  // Allocate
998  double interpolated_u = 0.0;
999  Vector<double> interpolated_x(2, 0.0);
1000 
1001  // Calculate position
1002  for (unsigned l = 0; l < n_node; l++)
1003  {
1004  // Get the value at the node
1005  double u_value = raw_nodal_value(l, u_index_spherical_adv_diff);
1006  interpolated_u += u_value * psif(l);
1007  // Loop over coordinate direction
1008  for (unsigned i = 0; i < 2; i++)
1009  {
1010  interpolated_x[i] += nodal_position(l, i) * psif(l);
1011  }
1012  }
1013 
1014  // Get the imposed beta (beta=flux when alpha=0.0)
1015  double beta;
1016  get_beta(interpolated_x, beta);
1017 
1018  // Get the imposed alpha
1019  double alpha;
1020  get_alpha(interpolated_x, alpha);
1021 
1022  // calculate the area weighting dS = r^{2} sin theta dr dtheta
1023  // r = x[0] and theta = x[1]
1024  double dS =
1025  interpolated_x[0] * interpolated_x[0] * sin(interpolated_x[1]);
1026 
1027  // Now add to the appropriate equations
1028 
1029  // Loop over the test functions
1030  for (unsigned l = 0; l < n_node; l++)
1031  {
1032  // Set the local equation number
1033  local_eqn = nodal_local_eqn(l, u_index_spherical_adv_diff);
1034  /*IF it's not a boundary condition*/
1035  if (local_eqn >= 0)
1036  {
1037  // Add the prescribed beta terms
1038  residuals[local_eqn] -=
1039  dS * (beta - alpha * interpolated_u) * testf(l) * W;
1040 
1041  // Calculate the Jacobian
1042  //----------------------
1043  if ((flag) && (alpha != 0.0))
1044  {
1045  // Loop over the velocity shape functions again
1046  for (unsigned l2 = 0; l2 < n_node; l2++)
1047  {
1048  // Set the number of the unknown
1049  local_unknown = nodal_local_eqn(l2, u_index_spherical_adv_diff);
1050 
1051  // If at a non-zero degree of freedom add in the entry
1052  if (local_unknown >= 0)
1053  {
1054  jacobian(local_eqn, local_unknown) +=
1055  dS * alpha * psif[l2] * testf[l] * W;
1056  }
1057  }
1058  }
1059  }
1060  } // end loop over test functions
1061 
1062  } // end loop over integration points
1063 
1064  } // end fill_in_generic_residual_contribution_adv_diff_flux
1065 
1066 
1067 } // namespace oomph
1068 
1069 #endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4342
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4630
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,...
Definition: elements.h:4501
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:5272
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:5358
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5002
A general Finite Element class.
Definition: elements.h:1317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2179
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:2597
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:3054
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...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2214
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1763
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:2321
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:3328
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:5162
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:3250
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2474
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:227
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:192
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
void output(FILE *file_pt)
C-style output function: r,z,u.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,z,u at n_plot^2 plot points.
QSphericalAdvectionDiffusionElement(const QSphericalAdvectionDiffusionElement< NNODE_1D > &dummy)=delete
Broken copy constructor.
void operator=(const QSphericalAdvectionDiffusionElement< NNODE_1D > &)=delete
Broken assignment operator.
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.
QSphericalAdvectionDiffusionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,z,u at n_plot^2 plot points.
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....
unsigned required_nvalue(const unsigned &n) const
Required # of ‘values’ (pinned or dofs) at node n.
void output(std::ostream &outfile)
Output function: r,z,u.
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.
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
A class for all elements that solve the Advection Diffusion equations in a spherical polar coordinate...
void disable_ALE()
Disable ALE – empty overload to suppress warning. ALE isn't implemented anyway.
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...
SphericalAdvectionDiffusionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
double interpolated_u_spherical_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
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 ...
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.
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
void output(std::ostream &outfile)
Output with default number of plot points.
double du_dt_spherical_adv_diff(const unsigned &n) const
du/dt at local node n. Uses suitably interpolated value for hanging nodes.
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux:
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
SphericalAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
void output(FILE *file_pt)
C_style output with default number of plot points.
SphericalAdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
SphericalAdvectionDiffusionEquations(const SphericalAdvectionDiffusionEquations &dummy)=delete
Broken copy constructor.
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...
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
const double & pe_st() const
Peclet number multiplied by Strouhal number.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
SphericalAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
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...
SphericalAdvectionDiffusionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
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(* SphericalAdvectionDiffusionSourceFctPt)(const Vector< double > &x, double &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
static double Default_peclet_number
Static default value for the Peclet number.
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.
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(* SphericalAdvectionDiffusionWindFctPt)(const Vector< double > &x, Vector< double > &wind)
Function pointer to wind function fct(x,w(x)) – x is a Vector!
SphericalAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
void operator=(const SphericalAdvectionDiffusionEquations &)=delete
Broken assignment operator.
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....
SphericalAdvectionDiffusionEquations()
Constructor: Initialise the Source_fct_pt and Wind_fct_pt to null and set (pointer to) Peclet number ...
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...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
unsigned U_index_adv_diff
The index at which the unknown is stored at the nodes.
SphericalAdvectionDiffusionPrescribedAlphaFctPt Alpha_fct_pt
Function pointer to the (global) prescribed-alpha function.
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!
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 ...
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.
SphericalAdvectionDiffusionPrescribedBetaFctPt & beta_fct_pt()
Access function for the prescribed-beta function pointer.
void get_beta(const Vector< double > &x, double &beta)
Function to calculate the prescribed beta at a given spatial position.
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...
void output(std::ostream &outfile, const unsigned &nplot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
void output(std::ostream &outfile)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
void(* SphericalAdvectionDiffusionPrescribedAlphaFctPt)(const Vector< double > &x, double &alpha)
Function pointer to the prescribed-alpha function fct(x,alpha(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...
SphericalAdvectionDiffusionPrescribedAlphaFctPt & alpha_fct_pt()
Access function for the prescribed-alpha function pointer.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
void operator=(const SphericalAdvectionDiffusionFluxElement &)=delete
Broken assignment operator.
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 ...
SphericalAdvectionDiffusionFluxElement(const SphericalAdvectionDiffusionFluxElement &dummy)=delete
Broken copy constructor.
SphericalAdvectionDiffusionPrescribedBetaFctPt Beta_fct_pt
Function pointer to the (global) prescribed-beta function.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:389
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...