steady_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 Steady axisymmetric advection diffusion elements
31 #ifndef OOMPH_STEADY_AXISYM_ADV_DIFF_ELEMENTS_HEADER
32 #define OOMPH_STEADY_AXISYM_ADV_DIFF_ELEMENTS_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
40 //OOMPH-LIB headers
41 #include "../generic/nodes.h"
42 #include "../generic/Qelements.h"
43 #include "../generic/refineable_elements.h"
44 #include "../generic/oomph_utilities.h"
45 
46 namespace oomph
47 {
48 
49 //=============================================================
50 /// \short A class for all elements that solve the Steady Axisymmetric
51 /// Advection Diffusion equations using isoparametric elements.
52 /// \f[
53 /// Pe \mathbf{w}\cdot(\mathbf{x}) \nabla u =
54 /// \nabla \cdot \left( \nabla u \right) + f(\mathbf{x})
55 /// \f]
56 /// This contains the generic maths. Shape functions, geometric
57 /// mapping etc. must get implemented in derived class.
58 //=============================================================
60 {
61 
62 public:
63 
64  /// \short Function pointer to source function fct(x,f(x)) --
65  /// x is a Vector!
67  double& f);
68 
69 
70  /// \short Function pointer to wind function fct(x,w(x)) --
71  /// x is a Vector!
73  Vector<double>& wind);
74 
75 
76  /// \short Constructor: Initialise the Source_fct_pt and Wind_fct_pt
77  /// to null and set (pointer to) Peclet number to default
79  {
80  //Set pointer to Peclet number to the default value zero
82  }
83 
84  /// Broken copy constructor
86  {
87  BrokenCopy::broken_copy("SteadyAxisymAdvectionDiffusionEquations");
88  }
89 
90  /// Broken assignment operator
91 //Commented out broken assignment operator because this can lead to a conflict warning
92 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
93 //realise that two separate implementations of the broken function are the same and so,
94 //quite rightly, it shouts.
95  /*void operator=(const SteadyAxisymAdvectionDiffusionEquations&)
96  {
97  BrokenCopy::broken_assign("SteadyAxisymAdvectionDiffusionEquations");
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 steady axisymmetric 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_axisym_adv_diff() const
108  {
109  return 0;
110  }
111 
112  /// Output with default number of plot points
113  void output(std::ostream &outfile)
114  {
115  unsigned nplot = 5;
116  output(outfile,nplot);
117  }
118 
119  /// \short Output FE representation of soln: r,z,u at
120  /// nplot^2 plot points
121  void output(std::ostream &outfile, const unsigned &nplot);
122 
123 
124  /// C_style output with default number of plot points
125  void output(FILE* file_pt)
126  {
127  unsigned n_plot = 5;
128  output(file_pt,n_plot);
129  }
130 
131  /// \short C-style output FE representation of soln: r,z,u at
132  /// n_plot^2 plot points
133  void output(FILE* file_pt, const unsigned &n_plot);
134 
135 
136  /// Output exact soln: r,z,u_exact at nplot^2 plot points
137  void output_fct(std::ostream &outfile,
138  const unsigned &nplot,
140  exact_soln_pt);
141 
142  /// Get error against and norm of exact solution
143  void compute_error(std::ostream &outfile,
145  exact_soln_pt, double& error,
146  double& norm);
147 
148  /// Access function: Pointer to source function
150  {
151  return Source_fct_pt;
152  }
153 
154 
155  /// Access function: Pointer to source function. Const version
157  {
158  return Source_fct_pt;
159  }
160 
161 
162  /// Access function: Pointer to wind function
164  {
165  return Wind_fct_pt;
166  }
167 
168 
169  /// Access function: Pointer to wind function. Const version
171  {
172  return Wind_fct_pt;
173  }
174 
175  // Access functions for the physical constants
176 
177  /// Peclet number
178  const double &pe() const
179  {
180  return *Pe_pt;
181  }
182 
183  /// Pointer to Peclet number
184  double* &pe_pt()
185  {
186  return Pe_pt;
187  }
188 
189  /// \short Get source term at (Eulerian) position x. This function is
190  /// virtual to allow overloading in multi-physics problems where
191  /// the strength of the source function might be determined by
192  /// another system of equations
193  inline virtual void get_source_axisym_adv_diff(const unsigned& ipt,
194  const Vector<double>& x,
195  double& source) const
196  {
197  //If no source function has been set, return zero
198  if(Source_fct_pt==0)
199  {
200  source = 0.0;
201  }
202  else
203  {
204  //Get source strength
205  (*Source_fct_pt)(x,source);
206  }
207  }
208 
209  /// \short Get wind at (Eulerian) position x and/or local coordinate s.
210  /// This function is
211  /// virtual to allow overloading in multi-physics problems where
212  /// the wind function might be determined by
213  /// another system of equations
214  inline virtual void get_wind_axisym_adv_diff(const unsigned& ipt,
215  const Vector<double> &s,
216  const Vector<double>& x,
217  Vector<double>& wind) const
218  {
219  //If no wind function has been set, return zero
220  if(Wind_fct_pt==0)
221  {
222  for(unsigned i=0;i<2;i++)
223  {
224  wind[i] = 0.0;
225  }
226  }
227  else
228  {
229  //Get wind
230  (*Wind_fct_pt)(x,wind);
231  }
232  }
233 
234  /// Get flux: \f$ \mbox{flux}[i] = \nabla u = \mbox{d}u / \mbox{d}x_i \f$
235  void get_flux(const Vector<double>& s, Vector<double>& flux) const
236  {
237  //Find out how many nodes there are in the element
238  unsigned n_node = nnode();
239 
240  //Get the nodal index at which the unknown is stored
241  unsigned u_nodal_index = u_index_axisym_adv_diff();
242 
243  //Set up memory for the shape and test functions
244  Shape psi(n_node);
245  DShape dpsidx(n_node,2);
246 
247  //Call the derivatives of the shape and test functions
248  dshape_eulerian(s,psi,dpsidx);
249 
250  //Initialise to zero
251  for(unsigned j=0;j<2;j++)
252  {
253  flux[j] = 0.0;
254  }
255 
256  //Loop over nodes
257  for(unsigned l=0;l<n_node;l++)
258  {
259  //Loop over derivative directions
260  for(unsigned j=0;j<2;j++)
261  {
262  flux[j] += nodal_value(l,u_nodal_index)*dpsidx(l,j);
263  }
264  }
265  }
266 
267 
268  /// Add the element's contribution to its residual vector (wrapper)
270  {
271  //Call the generic residuals function with flag set to 0 and using
272  //a dummy matrix
276  0);
277  }
278 
279 
280  /// \short Add the element's contribution to its residual vector and
281  /// the element Jacobian matrix (wrapper)
283  DenseMatrix<double> &jacobian)
284  {
285  //Call the generic routine with the flag set to 1
287  jacobian,
289  1);
290  }
291 
292 
293  /// Return FE representation of function value u(s) at local coordinate s
294  inline double interpolated_u_adv_diff(const Vector<double> &s) const
295  {
296  //Find number of nodes
297  unsigned n_node = nnode();
298 
299  //Get the nodal index at which the unknown is stored
300  unsigned u_nodal_index = u_index_axisym_adv_diff();
301 
302  //Local shape function
303  Shape psi(n_node);
304 
305  //Find values of shape function
306  shape(s,psi);
307 
308  //Initialise value of u
309  double interpolated_u = 0.0;
310 
311  //Loop over the local nodes and sum
312  for(unsigned l=0;l<n_node;l++)
313  {
314  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
315  }
316 
317  return(interpolated_u);
318  }
319 
320 
321  ///\short Return derivative of u at point s with respect to all data
322  ///that can affect its value.
323  ///In addition, return the global equation numbers corresponding to the
324  ///data. This is virtual so that it can be overloaded in the
325  ///refineable version
327  Vector<double> &du_ddata,
328  Vector<unsigned> &global_eqn_number)
329  {
330  //Find number of nodes
331  const unsigned n_node = nnode();
332 
333  //Get the nodal index at which the unknown is stored
334  const unsigned u_nodal_index = u_index_axisym_adv_diff();
335 
336  //Local shape function
337  Shape psi(n_node);
338 
339  //Find values of shape function
340  shape(s,psi);
341 
342  //Find the number of dofs associated with interpolated u
343  unsigned n_u_dof = 0;
344  for(unsigned l=0;l<n_node;l++)
345  {
346  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
347  //If it's positive add to the count
348  if (global_eqn >=0)
349  {
350  ++n_u_dof;
351  }
352  }
353 
354  //Now resize the storage schemes
355  du_ddata.resize(n_u_dof,0.0);
356  global_eqn_number.resize(n_u_dof,0);
357 
358  //Loop over the nodes again and set the derivatives
359  unsigned count = 0;
360  for(unsigned l=0;l<n_node;l++)
361  {
362  //Get the global equation number
363  int global_eqn=this->node_pt(l)->eqn_number(u_nodal_index);
364  //If it's positive
365  if (global_eqn >= 0)
366  {
367  //Set the global equation number
368  global_eqn_number[count] = global_eqn;
369  //Set the derivative with respect to the unknown
370  du_ddata[count] = psi[l];
371  //Increase the counter
372  ++count;
373  }
374  }
375  }
376 
377 
378  /// \short Self-test: Return 0 for OK
379  unsigned self_test();
380 
381 protected:
382 
383  /// \short Shape/test functions and derivs w.r.t. to global coords at
384  /// local coord. s; return Jacobian of mapping
385  virtual double dshape_and_dtest_eulerian_adv_diff(const Vector<double> &s,
386  Shape &psi,
387  DShape &dpsidx,
388  Shape &test,
389  DShape &dtestdx)
390  const = 0;
391 
392  /// \short Shape/test functions and derivs w.r.t. to global coords at
393  /// integration point ipt; return Jacobian of mapping
394  virtual double dshape_and_dtest_eulerian_at_knot_adv_diff(const unsigned &ipt,
395  Shape &psi,
396  DShape &dpsidx,
397  Shape &test,
398  DShape &dtestdx)
399  const = 0;
400 
401  /// \short Add the element's contribution to its residual vector only
402  /// (if flag=and/or element Jacobian matrix
404  DenseMatrix<double> &jacobian,
405  DenseMatrix<double> &mass_matrix,
406  unsigned flag);
407 
408  // Physical constants
409 
410  /// Pointer to global Peclet number
411  double *Pe_pt;
412 
413  /// Pointer to source function:
415 
416  /// Pointer to wind function:
418 
419  private:
420 
421  /// Static default value for the Peclet number
422  static double Default_peclet_number;
423 
424 
425 };//End class SteadyAxisymAdvectionDiffusionEquations
426 
427 
428 ///////////////////////////////////////////////////////////////////////////
429 ///////////////////////////////////////////////////////////////////////////
430 ///////////////////////////////////////////////////////////////////////////
431 
432 
433 
434 //======================================================================
435 /// \short QSteadyAxisymAdvectionDiffusionElement elements are
436 /// linear/quadrilateral/brick-shaped Axisymmetric Advection Diffusion
437 /// elements with isoparametric interpolation for the function.
438 //======================================================================
439 template <unsigned NNODE_1D>
440 class QSteadyAxisymAdvectionDiffusionElement : public virtual QElement<2,NNODE_1D>,
442 {
443 
444 private:
445 
446  /// \short Static array of ints to hold number of variables at
447  /// nodes: Initial_Nvalue[n]
448  static const unsigned Initial_Nvalue;
449 
450 public:
451 
452 
453  ///\short Constructor: Call constructors for QElement and
454  /// Advection Diffusion equations
457  { }
458 
459  /// Broken copy constructor
461  {
462  BrokenCopy::broken_copy("QSteadyAxisymAdvectionDiffusionElement");
463  }
464 
465  /// Broken assignment operator
466  /*void operator=(const QSteadyAxisymAdvectionDiffusionElement<NNODE_1D>&)
467  {
468  BrokenCopy::broken_assign("QSteadyAxisymAdvectionDiffusionElement");
469  }*/
470 
471  /// \short Required # of `values' (pinned or dofs)
472  /// at node n
473  inline unsigned required_nvalue(const unsigned &n) const
474  {
475  return Initial_Nvalue;
476  }
477 
478  /// \short Output function:
479  /// r,z,u
480  void output(std::ostream &outfile)
481  {
483  }
484 
485  /// \short Output function:
486  /// r,z,u at n_plot^2 plot points
487  void output(std::ostream &outfile, const unsigned &n_plot)
488  {
490  }
491 
492 
493  /// \short C-style output function:
494  /// r,z,u
495  void output(FILE* file_pt)
496  {
498  }
499 
500  /// \short C-style output function:
501  /// r,z,u at n_plot^2 plot points
502  void output(FILE* file_pt, const unsigned &n_plot)
503  {
505  }
506 
507  /// \short Output function for an exact solution:
508  /// r,z,u_exact at n_plot^2 plot points
509  void output_fct(std::ostream &outfile,
510  const unsigned &n_plot,
512  {
513  SteadyAxisymAdvectionDiffusionEquations::output_fct(outfile,n_plot,exact_soln_pt);
514  }
515 
516 
517 protected:
518 
519  /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
521  Shape &psi,
522  DShape &dpsidx,
523  Shape &test,
524  DShape &dtestdx) const;
525 
526  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
527  /// integration point ipt. Return Jacobian.
528  inline double dshape_and_dtest_eulerian_at_knot_adv_diff(const unsigned& ipt,
529  Shape &psi,
530  DShape &dpsidx,
531  Shape &test,
532  DShape &dtestdx) const;
533 
534 };//End class QSteadyAxisymAdvectionDiffusionElement
535 
536 //Inline functions:
537 
538 //======================================================================
539 /// \short Define the shape functions and test functions and derivatives
540 /// w.r.t. global coordinates and return Jacobian of mapping.
541 ///
542 /// Galerkin: Test functions = shape functions
543 //======================================================================
544 
545 template<unsigned NNODE_1D>
548  Shape &psi,
549  DShape &dpsidx,
550  Shape &test,
551  DShape &dtestdx) const
552 {
553  //Call the geometrical shape functions and derivatives
554  double J = this->dshape_eulerian(s,psi,dpsidx);
555 
556  //Loop over the test functions and derivatives and set them equal to the
557  //shape functions
558  for(unsigned i=0;i<NNODE_1D;i++)
559  {
560  test[i] = psi[i];
561  for(unsigned j=0;j<2;j++)
562  {
563  dtestdx(i,j) = dpsidx(i,j);
564  }
565  }
566 
567  //Return the jacobian
568  return J;
569 }
570 
571 
572 
573 //======================================================================
574 /// \short Define the shape functions and test functions and derivatives
575 /// w.r.t. global coordinates and return Jacobian of mapping.
576 ///
577 /// Galerkin: Test functions = shape functions
578 //======================================================================
579 
580 template<unsigned NNODE_1D>
583  Shape &psi,
584  DShape &dpsidx,
585  Shape &test,
586  DShape &dtestdx) const
587 {
588  //Call the geometrical shape functions and derivatives
589  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
590 
591  //Set the test functions equal to the shape functions (pointer copy)
592  test = psi;
593  dtestdx = dpsidx;
594 
595  //Return the jacobian
596  return J;
597 }
598 
599 ////////////////////////////////////////////////////////////////////////
600 ////////////////////////////////////////////////////////////////////////
601 ////////////////////////////////////////////////////////////////////////
602 
603 template<unsigned NNODE_1D>
605 public virtual QElement<1,NNODE_1D>
606 {
607 
608  public:
609 
610  /// \short Constructor: Call the constructor for the
611  /// appropriate lower-dimensional QElement
612  FaceGeometry() : QElement<1,NNODE_1D>() {}
613 
614 };
615 
616 
617 
618 
619 ////////////////////////////////////////////////////////////////////////
620 ////////////////////////////////////////////////////////////////////////
621 ////////////////////////////////////////////////////////////////////////
622 
623 //======================================================================
624 /// \short A class for elements that allow the imposition of an
625 /// applied Robin boundary condition on the boundaries of Steady
626 /// Axisymmnetric Advection Diffusion Flux elements.
627 /// \f[
628 /// -\Delta u \cdot \mathbf{n} + \alpha(r,z) u = \beta(r,z)
629 /// \f]
630 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
631 /// policy class.
632 //======================================================================
633 template <class ELEMENT>
635 public virtual FaceGeometry<ELEMENT>,
636  public virtual FaceElement
637 {
638 
639 public:
640 
641 
642  /// \short Function pointer to the prescribed-beta function fct(x,beta(x)) --
643  /// x is a Vector!
645  double& beta);
646 
647  /// \short Function pointer to the prescribed-alpha function fct(x,alpha(x)) --
648  /// x is a Vector!
650  double& alpha);
651 
652 
653  /// \short Constructor, takes the pointer to the "bulk" element
654  /// and the index of the face to be created
656  const int &face_index);
657 
658 
659  /// Broken empty constructor
661  {
662  throw OomphLibError(
663  "Don't call empty constructor for SteadyAxisymAdvectionDiffusionFluxElement",
664  OOMPH_CURRENT_FUNCTION,
665  OOMPH_EXCEPTION_LOCATION);
666  }
667 
668  /// Broken copy constructor
670  {
671  BrokenCopy::broken_copy("SteadyAxisymAdvectionDiffusionFluxElement");
672  }
673 
674  /// Broken assignment operator
675  /*void operator=(const SteadyAxisymAdvectionDiffusionFluxElement&)
676  {
677  BrokenCopy::broken_assign("SteadyAxisymAdvectionDiffusionFluxElement");
678  }*/
679 
680  /// Access function for the prescribed-beta function pointer
682  {
683  return Beta_fct_pt;
684  }
685 
686  /// Access function for the prescribed-alpha function pointer
688  {
689  return Alpha_fct_pt;
690  }
691 
692 
693  /// Add the element's contribution to its residual vector
695  {
696  //Call the generic residuals function with flag set to 0
697  //using a dummy matrix
700  0);
701  }
702 
703 
704  /// \short Add the element's contribution to its residual vector and
705  /// its Jacobian matrix
707  DenseMatrix<double> &jacobian)
708  {
709  //Call the generic routine with the flag set to 1
711  jacobian,
712  1);
713  }
714 
715 
716  /// Specify the value of nodal zeta from the face geometry
717  /// \short The "global" intrinsic coordinate of the element when
718  /// viewed as part of a geometric object should be given by
719  /// the FaceElement representation, by default (needed to break
720  /// indeterminacy if bulk element is SolidElement)
721  double zeta_nodal(const unsigned &n, const unsigned &k,
722  const unsigned &i) const
723  {return FaceElement::zeta_nodal(n,k,i);}
724 
725 
726  /// \short Output function -- forward to broken version in FiniteElement
727  /// until somebody decides what exactly they want to plot here...
728  void output(std::ostream &outfile)
729  {
730  FiniteElement::output(outfile);
731  }
732 
733  /// \short Output function -- forward to broken version in FiniteElement
734  /// until somebody decides what exactly they want to plot here...
735  void output(std::ostream &outfile, const unsigned &nplot)
736  {
737  FiniteElement::output(outfile,nplot);
738  }
739 
740 
741 protected:
742 
743  /// \short Function to compute the shape and test functions and to return
744  /// the Jacobian of mapping between local and global (Eulerian)
745  /// coordinates
746  inline double shape_and_test(const Vector<double> &s,
747  Shape &psi,
748  Shape &test) const
749  {
750  //Find number of nodes
751  unsigned n_node = nnode();
752 
753  //Get the shape functions
754  shape(s,psi);
755 
756  //Set the test functions to be the same as the shape functions
757  for(unsigned i=0;i<n_node;i++)
758  {
759  test[i] = psi[i];
760  }
761 
762  //Return the value of the jacobian
763  return J_eulerian(s);
764  }
765 
766 
767  /// \short Function to compute the shape and test functions and to return
768  /// the Jacobian of mapping between local and global (Eulerian)
769  /// coordinates
770  inline double shape_and_test_at_knot(const unsigned &ipt,
771  Shape &psi,
772  Shape &test) const
773  {
774  //Find number of nodes
775  unsigned n_node = nnode();
776 
777  //Get the shape functions
778  shape_at_knot(ipt,psi);
779 
780  //Set the test functions to be the same as the shape functions
781  for(unsigned i=0;i<n_node;i++)
782  {
783  test[i] = psi[i];
784  }
785 
786  //Return the value of the jacobian
787  return J_eulerian_at_knot(ipt);
788  }
789 
790  /// \short Function to calculate the prescribed beta at a given spatial
791  /// position
792  void get_beta(const Vector<double>& x, double& beta)
793  {
794  //If the function pointer is zero return zero
795  if(Beta_fct_pt == 0)
796  {
797  beta = 0.0;
798  }
799  //Otherwise call the function
800  else
801  {
802  (*Beta_fct_pt)(x,beta);
803  }
804  }
805 
806  /// \short Function to calculate the prescribed alpha at a given spatial
807  /// position
808  void get_alpha(const Vector<double>& x, double& alpha)
809  {
810  //If the function pointer is zero return zero
811  if(Alpha_fct_pt == 0)
812  {
813  alpha = 0.0;
814  }
815  //Otherwise call the function
816  else
817  {
818  (*Alpha_fct_pt)(x,alpha);
819  }
820  }
821 
822 private:
823 
824 
825  /// \short Add the element's contribution to its residual vector.
826  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
828  DenseMatrix<double> &jacobian,
829  unsigned flag);
830 
831 
832  /// Function pointer to the (global) prescribed-beta function
834 
835  /// Function pointer to the (global) prescribed-alpha function
837 
838  /// The index at which the unknown is stored at the nodes
840 
841 
842 };//End class SteadyAxisymAdvectionDiffusionFluxElement
843 
844 
845 ///////////////////////////////////////////////////////////////////////
846 ///////////////////////////////////////////////////////////////////////
847 ///////////////////////////////////////////////////////////////////////
848 
849 
850 //===========================================================================
851 /// \short Constructor, takes the pointer to the "bulk" element and the index
852 /// of the face to be created
853 //===========================================================================
854 template<class ELEMENT>
857  const int &face_index) :
858 FaceGeometry<ELEMENT>(), FaceElement()
859 
860 {
861  //Let the bulk element build the FaceElement, i.e. setup the pointers
862  //to its nodes (by referring to the appropriate nodes in the bulk
863  //element), etc.
864  bulk_el_pt->build_face_element(face_index,this);
865 
866 
867 #ifdef PARANOID
868  {
869  //Check that the element is not a refineable 3d element
870  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
871  //If it's three-d
872  if(elem_pt->dim()==3)
873  {
874  //Is it refineable
875  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
876  if(ref_el_pt!=0)
877  {
878  if (this->has_hanging_nodes())
879  {
880  throw OomphLibError(
881  "This flux element will not work correctly if nodes are hanging\n",
882  OOMPH_CURRENT_FUNCTION,
883  OOMPH_EXCEPTION_LOCATION);
884  }
885  }
886  }
887  }
888 #endif
889 
890  //Initialise the prescribed-beta function pointer to zero
891  Beta_fct_pt = 0;
892 
893  //Set up U_index_adv_diff. Initialise to zero, which probably won't change
894  //in most cases, oh well, the price we pay for generality
895  U_index_adv_diff = 0;
896 
897  //Cast to the appropriate AdvectionDiffusionEquation so that we can
898  //find the index at which the variable is stored
899  //We assume that the dimension of the full problem is the same
900  //as the dimension of the node, if this is not the case you will have
901  //to write custom elements, sorry
903  dynamic_cast<SteadyAxisymAdvectionDiffusionEquations*>(bulk_el_pt);
904 
905  //If the cast has failed die
906  if(eqn_pt==0)
907  {
908  std::string error_string =
909  "Bulk element must inherit from SteadyAxisymAdvectionDiffusionEquations.";
910  error_string +=
911  "Nodes are two dimensional, but cannot cast the bulk element to\n";
912  error_string += "SteadyAxisymAdvectionDiffusionEquations<2>\n.";
913  error_string +=
914  "If you desire this functionality, you must implement it yourself\n";
915 
916  throw OomphLibError(
917  error_string,
918  OOMPH_CURRENT_FUNCTION,
919  OOMPH_EXCEPTION_LOCATION);
920  }
921  else
922  {
923  //Read the index from the (cast) bulk element.
925  }
926 
927 
928 }
929 
930 
931 //===========================================================================
932 /// \short Compute the element's residual vector and the (zero) Jacobian
933 /// matrix for the Robin boundary condition:
934 /// \f[
935 /// \Delta u \cdot \mathbf{n} + \alpha (\mathbf{x}) = \beta (\mathbf{x})
936 /// \f]
937 //===========================================================================
938 template<class ELEMENT>
941  DenseMatrix<double> &jacobian,
942  unsigned flag)
943 {
944  //Find out how many nodes there are
945  const unsigned n_node = nnode();
946 
947  // Locally cache the index at which the variable is stored
948  const unsigned u_index_axisym_adv_diff = U_index_adv_diff;
949 
950  //Set up memory for the shape and test functions
951  Shape psif(n_node), testf(n_node);
952 
953  //Set the value of n_intpt
954  const unsigned n_intpt = integral_pt()->nweight();
955 
956  //Set the Vector to hold local coordinates
957  Vector<double> s(1);
958 
959  //Integers used to store the local equation number and local unknown
960  //indices for the residuals and jacobians
961  int local_eqn=0, local_unknown=0;
962 
963  //Loop over the integration points
964  //--------------------------------
965  for(unsigned ipt=0;ipt<n_intpt;ipt++)
966  {
967 
968  //Assign values of s
969  for(unsigned i=0;i<1;i++)
970  {
971  s[i] = integral_pt()->knot(ipt,i);
972  }
973 
974  //Get the integral weight
975  double w = integral_pt()->weight(ipt);
976 
977  //Find the shape and test functions and return the Jacobian
978  //of the mapping
979  double J = shape_and_test(s,psif,testf);
980 
981  //Premultiply the weights and the Jacobian
982  double W = w*J;
983 
984  //Calculate local values of the solution and its derivatives
985  //Allocate
986  double interpolated_u=0.0;
987  Vector<double> interpolated_x(2,0.0);
988 
989  //Calculate position
990  for(unsigned l=0;l<n_node;l++)
991  {
992  //Get the value at the node
993  double u_value = raw_nodal_value(l,u_index_axisym_adv_diff);
994  interpolated_u += u_value*psif(l);
995  //Loop over coordinate direction
996  for(unsigned i=0;i<2;i++)
997  {
998  interpolated_x[i] += nodal_position(l,i)*psif(l);
999  }
1000  }
1001 
1002  //Get the imposed beta (beta=flux when alpha=0.0)
1003  double beta;
1004  get_beta(interpolated_x,beta);
1005 
1006  //Get the imposed alpha
1007  double alpha;
1008  get_alpha(interpolated_x,alpha);
1009 
1010  //r is the first position component
1011  double r = interpolated_x[0];
1012 
1013  //Now add to the appropriate equations
1014 
1015  //Loop over the test functions
1016  for(unsigned l=0;l<n_node;l++)
1017  {
1018  //Set the local equation number
1019  local_eqn = nodal_local_eqn(l,u_index_axisym_adv_diff);
1020  /*IF it's not a boundary condition*/
1021  if(local_eqn >= 0)
1022  {
1023  //Add the prescribed beta terms
1024  residuals[local_eqn] -= r*(beta - alpha*interpolated_u)*testf(l)*W;
1025 
1026  //Calculate the Jacobian
1027  //----------------------
1028  if ( (flag) && (alpha!=0.0) )
1029  {
1030  //Loop over the velocity shape functions again
1031  for(unsigned l2=0;l2<n_node;l2++)
1032  {
1033  //Set the number of the unknown
1034  local_unknown = nodal_local_eqn(l2,u_index_axisym_adv_diff);
1035 
1036  //If at a non-zero degree of freedom add in the entry
1037  if(local_unknown >= 0)
1038  {
1039  jacobian(local_eqn,local_unknown) +=
1040  r*alpha*psif[l2]*testf[l]*W;
1041  }
1042  }
1043  }
1044 
1045  }
1046  } //end loop over test functions
1047 
1048  } //end loop over integration points
1049 
1050 }//end fill_in_generic_residual_contribution_adv_diff_flux
1051 
1052 
1053 } //End SteadyAxisymAdvectionDiffusionFluxElement class
1054 
1055 #endif
virtual void get_wind_axisym_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 broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
SteadyAxisymAdvectionDiffusionEquations()
Constructor: Initialise the Source_fct_pt and Wind_fct_pt to null and set (pointer to) Peclet number ...
virtual double dshape_and_dtest_eulerian_at_knot_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 ...
SteadyAxisymAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
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 ...
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
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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
void output(FILE *file_pt)
C_style output with default number of plot points.
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
void output(std::ostream &outfile, const unsigned &nplot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
A class for all elements that solve the Steady Axisymmetric Advection Diffusion equations using isopa...
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 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 ...
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
A general Finite Element class.
Definition: elements.h:1271
QSteadyAxisymAdvectionDiffusionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
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.
void output(std::ostream &outfile)
Output function: r,z,u.
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
void(* SteadyAxisymAdvectionDiffusionSourceFctPt)(const Vector< double > &x, double &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
SteadyAxisymAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
QSteadyAxisymAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Axisymmetric Ad...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,z,u at n_plot^2 plot points.
void output(std::ostream &outfile)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
SteadyAxisymAdvectionDiffusionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
void(* SteadyAxisymAdvectionDiffusionPrescribedAlphaFctPt)(const Vector< double > &x, double &alpha)
Function pointer to the prescribed-alpha function fct(x,alpha(x)) – x is a Vector! ...
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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
SteadyAxisymAdvectionDiffusionEquations(const SteadyAxisymAdvectionDiffusionEquations &dummy)
Broken copy constructor.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,z,u at n_plot^2 plot points.
double dshape_and_dtest_eulerian_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.
SteadyAxisymAdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
static char t char * s
Definition: cfortran.h:572
void(* SteadyAxisymAdvectionDiffusionPrescribedBetaFctPt)(const Vector< double > &x, double &beta)
Function pointer to the prescribed-beta function fct(x,beta(x)) – x is a Vector!
SteadyAxisymAdvectionDiffusionFluxElement(const SteadyAxisymAdvectionDiffusionFluxElement &dummy)
Broken copy constructor.
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.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
double dshape_and_dtest_eulerian_at_knot_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.
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
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
SteadyAxisymAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
SteadyAxisymAdvectionDiffusionPrescribedAlphaFctPt Alpha_fct_pt
Function pointer to the (global) prescribed-alpha function.
void fill_in_generic_residual_contribution_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...
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 get_beta(const Vector< double > &x, double &beta)
Function to calculate the prescribed beta at a given spatial position.
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
QSteadyAxisymAdvectionDiffusionElement(const QSteadyAxisymAdvectionDiffusionElement< NNODE_1D > &dummy)
Broken copy constructor.
SteadyAxisymAdvectionDiffusionPrescribedBetaFctPt Beta_fct_pt
Function pointer to the (global) prescribed-beta function.
virtual double dshape_and_dtest_eulerian_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...
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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) ...
static double Default_peclet_number
Static default value for the Peclet number.
virtual unsigned u_index_axisym_adv_diff() const
Broken assignment operator.
void(* SteadyAxisymAdvectionDiffusionWindFctPt)(const Vector< double > &x, Vector< double > &wind)
Function pointer to wind function fct(x,w(x)) – x is a Vector!
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2344
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
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_alpha(const Vector< double > &x, double &alpha)
Function to calculate the prescribed alpha at a given spatial position.
A class for elements that allow the imposition of an applied Robin boundary condition on the boundari...
virtual void get_source_axisym_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...
SteadyAxisymAdvectionDiffusionPrescribedBetaFctPt & beta_fct_pt()
Broken assignment operator.
void output(FILE *file_pt)
C-style output function: r,z,u.
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 ...
SteadyAxisymAdvectionDiffusionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
SteadyAxisymAdvectionDiffusionPrescribedAlphaFctPt & alpha_fct_pt()
Access function for the prescribed-alpha function pointer.
double interpolated_u_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
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...
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.
void output(std::ostream &outfile)
Output with default number of plot points.
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
virtual void fill_in_generic_residual_contribution_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...
unsigned U_index_adv_diff
The index at which the unknown is stored at the nodes.