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-2022 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
43namespace 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
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
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
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
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
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
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
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:4497
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:5242
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:5328
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
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:2593
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:3050
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:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
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:2317
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:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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:5132
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:3220
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2470
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.
SphericalAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: 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 ...
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 *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
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:
const double & pe_st() const
Peclet number multiplied by Strouhal number.
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...
SphericalAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
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!
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.
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...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
SphericalAdvectionDiffusionPrescribedBetaFctPt & beta_fct_pt()
Access function for the prescribed-beta function pointer.
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 ...
SphericalAdvectionDiffusionPrescribedAlphaFctPt & alpha_fct_pt()
Access function for the prescribed-alpha function pointer.
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...