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