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 advection diffusion elements in a spherical polar coordinate
27// system
28#ifndef OOMPH_AXISYM_ADV_DIFF_ELEMENTS_HEADER
29#define OOMPH_AXISYM_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 cylindrical 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 double& f);
60
61
62 /// Function pointer to wind function fct(x,w(x)) --
63 /// x is a Vector!
65 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
72 {
73 // Set pointer to Peclet number to the default value zero
77 }
78
79 /// Broken copy constructor
81 const AxisymAdvectionDiffusionEquations& dummy) = delete;
82
83 /// Broken assignment operator
84 // Commented out broken assignment operator because this can lead to a
85 // conflict warning when used in the virtual inheritence hierarchy.
86 // Essentially the compiler doesn't realise that two separate
87 // implementations of the broken function are the same and so, quite
88 // rightly, it shouts.
89 /*void operator=(const AxisymAdvectionDiffusionEquations&) = delete;*/
90
91 /// Return the index at which the unknown value
92 /// is stored. The default value, 0, is appropriate for single-physics
93 /// problems, when there is only one variable, the value that satisfies
94 /// the spherical advection-diffusion equation.
95 /// In derived multi-physics elements, this function should be overloaded
96 /// to reflect the chosen storage scheme. Note that these equations require
97 /// that the unknown is always stored at the same index at each node.
98 virtual inline unsigned u_index_axi_adv_diff() const
99 {
100 return 0;
101 }
102
103
104 /// du/dt at local node n.
105 /// Uses suitably interpolated value for hanging nodes.
106 double du_dt_axi_adv_diff(const unsigned& n) const
107 {
108 // Get the data's timestepper
110
111 // Initialise dudt
112 double dudt = 0.0;
113 // Loop over the timesteps, if there is a non Steady timestepper
115 {
116 // Find the index at which the variable is stored
117 const unsigned u_nodal_index = u_index_axi_adv_diff();
118
119 // Number of timsteps (past & present)
120 const unsigned n_time = time_stepper_pt->ntstorage();
121
122 for (unsigned t = 0; t < n_time; t++)
123 {
124 dudt +=
125 time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
126 }
127 }
128 return dudt;
129 }
130
131 /// Disable ALE, i.e. assert the mesh is not moving -- you do this
132 /// at your own risk!
134 {
135 ALE_is_disabled = true;
136 }
137
138
139 /// (Re-)enable ALE, i.e. take possible mesh motion into account
140 /// when evaluating the time-derivative. Note: By default, ALE is
141 /// enabled, at the expense of possibly creating unnecessary work
142 /// in problems where the mesh is, in fact, stationary.
144 {
145 ALE_is_disabled = false;
146 }
147
148
149 /// Number of scalars/fields output by this element. Reimplements
150 /// broken virtual function in base class.
151 unsigned nscalar_paraview() const
152 {
153 return 4;
154 }
155
156 /// Write values of the i-th scalar field at the plot points. Needs
157 /// to be implemented for each new specific element type.
158 void scalar_value_paraview(std::ofstream& file_out,
159 const unsigned& i,
160 const unsigned& nplot) const
161 {
162 // Vector of local coordinates
163 Vector<double> s(2);
164
165 // Loop over plot points
166 unsigned num_plot_points = nplot_points_paraview(nplot);
167 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
168 {
169 // Get local coordinates of plot point
170 get_s_plot(iplot, nplot, s);
171
172 // Get Eulerian coordinate of plot point
173 Vector<double> x(2);
174 interpolated_x(s, x);
175
176 // Winds
177 if (i < 3)
178 {
179 // Get the wind
180 Vector<double> wind(3);
181
182 // Dummy ipt argument needed... ?
183 unsigned ipt = 0;
184 get_wind_axi_adv_diff(ipt, s, x, wind);
185
186 file_out << wind[i] << std::endl;
187 }
188 // Advection Diffusion
189 else if (i == 3)
190 {
191 file_out << this->interpolated_u_axi_adv_diff(s) << std::endl;
192 }
193 // Never get here
194 else
195 {
196 std::stringstream error_stream;
197 error_stream << "Advection Diffusion Elements only store 4 fields "
198 << std::endl;
199 throw OomphLibError(error_stream.str(),
200 OOMPH_CURRENT_FUNCTION,
201 OOMPH_EXCEPTION_LOCATION);
202 }
203 }
204 }
205
206 /// Name of the i-th scalar field. Default implementation
207 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
208 /// overloaded with more meaningful names in specific elements.
209 std::string scalar_name_paraview(const unsigned& i) const
210 {
211 // Winds
212 if (i < 3)
213 {
214 return "Wind " + StringConversion::to_string(i);
215 }
216 // Advection Diffusion field
217 else if (i == 3)
218 {
219 return "Advection Diffusion";
220 }
221 // Never get here
222 else
223 {
224 std::stringstream error_stream;
225 error_stream << "Advection Diffusion Elements only store 4 fields "
226 << std::endl;
227 throw OomphLibError(
228 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
229 // Dummy return
230 return " ";
231 }
232 }
233
234 /// Output with default number of plot points
235 void output(std::ostream& outfile)
236 {
237 unsigned nplot = 5;
238 output(outfile, nplot);
239 }
240
241 /// Output FE representation of soln: r,z,u at
242 /// nplot^2 plot points
243 void output(std::ostream& outfile, const unsigned& nplot);
244
245
246 /// C_style output with default number of plot points
247 void output(FILE* file_pt)
248 {
249 unsigned n_plot = 5;
250 output(file_pt, n_plot);
251 }
252
253 /// C-style output FE representation of soln: r,z,u at
254 /// n_plot^2 plot points
255 void output(FILE* file_pt, const unsigned& n_plot);
256
257
258 /// Output exact soln: r,z,u_exact at nplot^2 plot points
259 void output_fct(std::ostream& outfile,
260 const unsigned& nplot,
262
263 /// Get error against and norm of exact solution
264 void compute_error(std::ostream& outfile,
266 double& error,
267 double& norm);
268
269 /// Access function: Pointer to source function
271 {
272 return Source_fct_pt;
273 }
274
275
276 /// Access function: Pointer to source function. Const version
278 {
279 return Source_fct_pt;
280 }
281
282
283 /// Access function: Pointer to wind function
285 {
286 return Wind_fct_pt;
287 }
288
289
290 /// Access function: Pointer to wind function. Const version
292 {
293 return Wind_fct_pt;
294 }
295
296 // Access functions for the physical constants
297
298 /// Peclet number
299 inline const double& pe() const
300 {
301 return *Pe_pt;
302 }
303
304 /// Pointer to Peclet number
305 inline double*& pe_pt()
306 {
307 return Pe_pt;
308 }
309
310 /// Peclet number multiplied by Strouhal number
311 inline const double& pe_st() const
312 {
313 return *PeSt_pt;
314 }
315
316 /// Pointer to Peclet number multipled by Strouha number
317 inline double*& pe_st_pt()
318 {
319 return PeSt_pt;
320 }
321
322 /// Peclet number multiplied by Strouhal number
323 inline const double& d() const
324 {
325 return *D_pt;
326 }
327
328 /// Pointer to Peclet number multipled by Strouha number
329 inline double*& d_pt()
330 {
331 return D_pt;
332 }
333
334 /// Get source term at (Eulerian) position x. This function is
335 /// virtual to allow overloading in multi-physics problems where
336 /// the strength of the source function might be determined by
337 /// another system of equations
338 inline virtual void get_source_axi_adv_diff(const unsigned& ipt,
339 const Vector<double>& x,
340 double& source) const
341 {
342 // If no source function has been set, return zero
343 if (Source_fct_pt == 0)
344 {
345 source = 0.0;
346 }
347 else
348 {
349 // Get source strength
350 (*Source_fct_pt)(x, source);
351 }
352 }
353
354 /// Get wind at (Eulerian) position x and/or local coordinate s.
355 /// This function is
356 /// virtual to allow overloading in multi-physics problems where
357 /// the wind function might be determined by
358 /// another system of equations
359 inline virtual void get_wind_axi_adv_diff(const unsigned& ipt,
360 const Vector<double>& s,
361 const Vector<double>& x,
362 Vector<double>& wind) const
363 {
364 // If no wind function has been set, return zero
365 if (Wind_fct_pt == 0)
366 {
367 for (unsigned i = 0; i < 3; i++)
368 {
369 wind[i] = 0.0;
370 }
371 }
372 else
373 {
374 // Get wind
375 (*Wind_fct_pt)(x, wind);
376 }
377 }
378
379 /// Get flux: [du/dr,du/dz]
380 void get_flux(const Vector<double>& s, Vector<double>& flux) const
381 {
382 // Find out how many nodes there are in the element
383 const unsigned n_node = nnode();
384
385 // Get the nodal index at which the unknown is stored
386 const unsigned u_nodal_index = u_index_axi_adv_diff();
387
388 // Set up memory for the shape and test functions
389 Shape psi(n_node);
390 DShape dpsidx(n_node, 2);
391
392 // Call the derivatives of the shape and test functions
393 dshape_eulerian(s, psi, dpsidx);
394
395 // Initialise to zero
396 for (unsigned j = 0; j < 2; j++)
397 {
398 flux[j] = 0.0;
399 }
400
401 // Loop over nodes
402 for (unsigned l = 0; l < n_node; l++)
403 {
404 const double u_value = this->nodal_value(l, u_nodal_index);
405 // Add in the derivative directions
406 flux[0] += u_value * dpsidx(l, 0);
407 flux[1] += u_value * dpsidx(l, 1);
408 }
409 }
410
411
412 /// Add the element's contribution to its residual vector (wrapper)
414 {
415 // Call the generic residuals function with flag set to 0 and using
416 // a dummy matrix
418 residuals,
421 0);
422 }
423
424
425 /// Add the element's contribution to its residual vector and
426 /// the element Jacobian matrix (wrapper)
428 DenseMatrix<double>& jacobian)
429 {
430 // Call the generic routine with the flag set to 1
432 residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
433 }
434
435
436 /// Return FE representation of function value u(s) at local coordinate s
437 inline double interpolated_u_axi_adv_diff(const Vector<double>& s) const
438 {
439 // Find number of nodes
440 const unsigned n_node = nnode();
441
442 // Get the nodal index at which the unknown is stored
443 const unsigned u_nodal_index = u_index_axi_adv_diff();
444
445 // Local shape function
446 Shape psi(n_node);
447
448 // Find values of shape function
449 shape(s, psi);
450
451 // Initialise value of u
452 double interpolated_u = 0.0;
453
454 // Loop over the local nodes and sum
455 for (unsigned l = 0; l < n_node; l++)
456 {
457 interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
458 }
459
460 return (interpolated_u);
461 }
462
463
464 /// Return derivative of u at point s with respect to all data
465 /// that can affect its value.
466 /// In addition, return the global equation numbers corresponding to the
467 /// data. This is virtual so that it can be overloaded in the
468 /// refineable version
470 const Vector<double>& s,
471 Vector<double>& du_ddata,
472 Vector<unsigned>& global_eqn_number)
473 {
474 // Find number of nodes
475 const unsigned n_node = nnode();
476
477 // Get the nodal index at which the unknown is stored
478 const unsigned u_nodal_index = u_index_axi_adv_diff();
479
480 // Local shape function
481 Shape psi(n_node);
482
483 // Find values of shape function
484 shape(s, psi);
485
486 // Find the number of dofs associated with interpolated u
487 unsigned n_u_dof = 0;
488 for (unsigned l = 0; l < n_node; l++)
489 {
490 int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
491 // If it's positive add to the count
492 if (global_eqn >= 0)
493 {
494 ++n_u_dof;
495 }
496 }
497
498 // Now resize the storage schemes
499 du_ddata.resize(n_u_dof, 0.0);
500 global_eqn_number.resize(n_u_dof, 0);
501
502 // Loop over the nodes again and set the derivatives
503 unsigned count = 0;
504 for (unsigned l = 0; l < n_node; l++)
505 {
506 // Get the global equation number
507 int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
508 // If it's positive
509 if (global_eqn >= 0)
510 {
511 // Set the global equation number
512 global_eqn_number[count] = global_eqn;
513 // Set the derivative with respect to the unknown
514 du_ddata[count] = psi[l];
515 // Increase the counter
516 ++count;
517 }
518 }
519 }
520
521
522 /// Self-test: Return 0 for OK
523 unsigned self_test();
524
525 protected:
526 /// Shape/test functions and derivs w.r.t. to global coords at
527 /// local coord. s; return Jacobian of mapping
529 const Vector<double>& s,
530 Shape& psi,
531 DShape& dpsidx,
532 Shape& test,
533 DShape& dtestdx) const = 0;
534
535 /// Shape/test functions and derivs w.r.t. to global coords at
536 /// integration point ipt; return Jacobian of mapping
538 const unsigned& ipt,
539 Shape& psi,
540 DShape& dpsidx,
541 Shape& test,
542 DShape& dtestdx) const = 0;
543
544 /// Add the element's contribution to its residual vector only
545 /// (if flag=and/or element Jacobian matrix
547 Vector<double>& residuals,
548 DenseMatrix<double>& jacobian,
549 DenseMatrix<double>& mass_matrix,
550 unsigned flag);
551
552 // Physical constants
553
554 /// Pointer to global Peclet number
555 double* Pe_pt;
556
557 /// Pointer to global Peclet number multiplied by Strouhal number
558 double* PeSt_pt;
559
560 /// Pointer to global Diffusion parameter
561 double* D_pt;
562
563 /// Pointer to source function:
565
566 /// Pointer to wind function:
568
569 /// Boolean flag to indicate whether AlE formulation is disable
571
572 private:
573 /// Static default value for the Peclet number
575
576 /// Static default value for the Peclet number
578 ;
579
580
581 }; // End class AxisymAdvectionDiffusionEquations
582
583
584 /// ////////////////////////////////////////////////////////////////////////
585 /// ////////////////////////////////////////////////////////////////////////
586 /// ////////////////////////////////////////////////////////////////////////
587
588
589 //======================================================================
590 /// QAxisymAdvectionDiffusionElement elements are
591 /// linear/quadrilateral/brick-shaped Axisymmetric Advection Diffusion
592 /// elements with isoparametric interpolation for the function.
593 //======================================================================
594 template<unsigned NNODE_1D>
596 : public virtual QElement<2, NNODE_1D>,
598 {
599 private:
600 /// Static array of ints to hold number of variables at
601 /// nodes: Initial_Nvalue[n]
602 static const unsigned Initial_Nvalue;
603
604 public:
605 /// Constructor: Call constructors for QElement and
606 /// Advection Diffusion equations
609 {
610 }
611
612 /// Broken copy constructor
614 const QAxisymAdvectionDiffusionElement<NNODE_1D>& dummy) = delete;
615
616 /// Broken assignment operator
617 /*void operator=(const QAxisymAdvectionDiffusionElement<NNODE_1D>&) =
618 * delete;*/
619
620 /// Required # of `values' (pinned or dofs)
621 /// at node n
622 inline unsigned required_nvalue(const unsigned& n) const
623 {
624 return Initial_Nvalue;
625 }
626
627 /// Output function:
628 /// r,z,u
629 void output(std::ostream& outfile)
630 {
632 }
633
634 /// Output function:
635 /// r,z,u at n_plot^2 plot points
636 void output(std::ostream& outfile, const unsigned& n_plot)
637 {
639 }
640
641
642 /// C-style output function:
643 /// r,z,u
644 void output(FILE* file_pt)
645 {
647 }
648
649 /// C-style output function:
650 /// r,z,u at n_plot^2 plot points
651 void output(FILE* file_pt, const unsigned& n_plot)
652 {
654 }
655
656 /// Output function for an exact solution:
657 /// r,z,u_exact at n_plot^2 plot points
658 void output_fct(std::ostream& outfile,
659 const unsigned& n_plot,
661 {
663 outfile, n_plot, exact_soln_pt);
664 }
665
666
667 protected:
668 /// Shape, test functions & derivs. w.r.t. to global coords. Return
669 /// Jacobian.
671 const Vector<double>& s,
672 Shape& psi,
673 DShape& dpsidx,
674 Shape& test,
675 DShape& dtestdx) const;
676
677 /// Shape, test functions & derivs. w.r.t. to global coords. at
678 /// integration point ipt. Return Jacobian.
680 const unsigned& ipt,
681 Shape& psi,
682 DShape& dpsidx,
683 Shape& test,
684 DShape& dtestdx) const;
685
686 }; // End class QAxisymAdvectionDiffusionElement
687
688 // Inline functions:
689
690 //======================================================================
691 /// Define the shape functions and test functions and derivatives
692 /// w.r.t. global coordinates and return Jacobian of mapping.
693 ///
694 /// Galerkin: Test functions = shape functions
695 //======================================================================
696
697 template<unsigned NNODE_1D>
698 double QAxisymAdvectionDiffusionElement<
699 NNODE_1D>::dshape_and_dtest_eulerian_axi_adv_diff(const Vector<double>& s,
700 Shape& psi,
701 DShape& dpsidx,
702 Shape& test,
703 DShape& dtestdx) const
704 {
705 // Call the geometrical shape functions and derivatives
706 double J = this->dshape_eulerian(s, psi, dpsidx);
707
708 // Loop over the test functions and derivatives and set them equal to the
709 // shape functions
710 for (unsigned i = 0; i < NNODE_1D; i++)
711 {
712 test[i] = psi[i];
713 for (unsigned j = 0; j < 2; j++)
714 {
715 dtestdx(i, j) = dpsidx(i, j);
716 }
717 }
718
719 // Return the jacobian
720 return J;
721 }
722
723
724 //======================================================================
725 /// Define the shape functions and test functions and derivatives
726 /// w.r.t. global coordinates and return Jacobian of mapping.
727 ///
728 /// Galerkin: Test functions = shape functions
729 //======================================================================
730
731 template<unsigned NNODE_1D>
734 Shape& psi,
735 DShape& dpsidx,
736 Shape& test,
737 DShape& dtestdx) const
738 {
739 // Call the geometrical shape functions and derivatives
740 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
741
742 // Set the test functions equal to the shape functions (pointer copy)
743 test = psi;
744 dtestdx = dpsidx;
745
746 // Return the jacobian
747 return J;
748 }
749
750 /// /////////////////////////////////////////////////////////////////////
751 /// /////////////////////////////////////////////////////////////////////
752 /// /////////////////////////////////////////////////////////////////////
753
754 template<unsigned NNODE_1D>
756 : public virtual QElement<1, NNODE_1D>
757 {
758 public:
759 /// Constructor: Call the constructor for the
760 /// appropriate lower-dimensional QElement
761 FaceGeometry() : QElement<1, NNODE_1D>() {}
762 };
763
764
765 /// /////////////////////////////////////////////////////////////////////
766 /// /////////////////////////////////////////////////////////////////////
767 /// /////////////////////////////////////////////////////////////////////
768
769 //======================================================================
770 /// A class for elements that allow the imposition of an
771 /// applied Robin boundary condition on the boundaries of Steady
772 /// Axisymmnetric Advection Diffusion Flux elements.
773 /// \f[ -\Delta u \cdot \mathbf{n} + \alpha(r,z) u = \beta(r,z) \f]
774 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
775 /// policy class.
776 //======================================================================
777 /// AT THE MOMENT THESE ARE IN SPHERICAL POLARS --- NOT CONSISTENT
778
779 // template <class ELEMENT>
780 // class AxisymAdvectionDiffusionFluxElement :
781 // public virtual FaceGeometry<ELEMENT>,
782 // public virtual FaceElement
783 // {
784
785 // public:
786
787
788 // /// Function pointer to the prescribed-beta function fct(x,beta(x))
789 // --
790 // /// x is a Vector!
791 // typedef void (*AxisymAdvectionDiffusionPrescribedBetaFctPt)(
792 // const Vector<double>& x,
793 // double& beta);
794
795 // /// Function pointer to the prescribed-alpha function
796 // fct(x,alpha(x)) --
797 // /// x is a Vector!
798 // typedef void (*AxisymAdvectionDiffusionPrescribedAlphaFctPt)(
799 // const Vector<double>& x,
800 // double& alpha);
801
802
803 // /// Constructor, takes the pointer to the "bulk" element
804 // /// and the index of the face to be created
805 // AxisymAdvectionDiffusionFluxElement(FiniteElement* const &bulk_el_pt,
806 // const int &face_index);
807
808
809 // /// Broken empty constructor
810 // AxisymAdvectionDiffusionFluxElement()
811 // {
812 // throw OomphLibError(
813 // "Don't call empty constructor for AxisymAdvectionDiffusionFluxElement",
814 // OOMPH_CURRENT_FUNCTION,
815 // OOMPH_EXCEPTION_LOCATION);
816 // }
817
818 // /// Broken copy constructor
819 // AxisymAdvectionDiffusionFluxElement(
820 // const AxisymAdvectionDiffusionFluxElement& dummy) = delete;
821
822 // /// Broken assignment operator
823 // /*void operator=(const AxisymAdvectionDiffusionFluxElement&) = delete;*/
824
825 // /// Access function for the prescribed-beta function pointer
826 // AxisymAdvectionDiffusionPrescribedBetaFctPt& beta_fct_pt()
827 // {
828 // return Beta_fct_pt;
829 // }
830
831 // /// Access function for the prescribed-alpha function pointer
832 // AxisymAdvectionDiffusionPrescribedAlphaFctPt& alpha_fct_pt()
833 // {
834 // return Alpha_fct_pt;
835 // }
836
837
838 // /// Add the element's contribution to its residual vector
839 // inline void fill_in_contribution_to_residuals(Vector<double> &residuals)
840 // {
841 // //Call the generic residuals function with flag set to 0
842 // //using a dummy matrix
843 // fill_in_generic_residual_contribution_axi_adv_diff_flux(
844 // residuals,
845 // GeneralisedElement::Dummy_matrix,0);
846 // }
847
848
849 // /// Add the element's contribution to its residual vector and
850 // /// its Jacobian matrix
851 // inline void fill_in_contribution_to_jacobian(Vector<double> &residuals,
852 // DenseMatrix<double>
853 // &jacobian)
854 // {
855 // //Call the generic routine with the flag set to 1
856 // fill_in_generic_residual_contribution_axi_adv_diff_flux(residuals,
857 // jacobian,
858 // 1);
859 // }
860
861 // /// Specify the value of nodal zeta from the face geometry
862 // /// The "global" intrinsic coordinate of the element when
863 // /// viewed as part of a geometric object should be given by
864 // /// the FaceElement representation, by default (needed to break
865 // /// indeterminacy if bulk element is SolidElement)
866 // double zeta_nodal(const unsigned &n, const unsigned &k,
867 // const unsigned &i) const
868 // {return FaceElement::zeta_nodal(n,k,i);}
869
870
871 // /// Output function -- forward to broken version in FiniteElement
872 // /// until somebody decides what exactly they want to plot here...
873 // void output(std::ostream &outfile)
874 // {
875 // FiniteElement::output(outfile);
876 // }
877
878 // /// Output function -- forward to broken version in FiniteElement
879 // /// until somebody decides what exactly they want to plot here...
880 // void output(std::ostream &outfile, const unsigned &nplot)
881 // {
882 // FiniteElement::output(outfile,nplot);
883 // }
884
885
886 // protected:
887
888 // /// Function to compute the shape and test functions and to return
889 // /// the Jacobian of mapping between local and global (Eulerian)
890 // /// coordinates
891 // inline double shape_and_test(const Vector<double> &s,
892 // Shape &psi,
893 // Shape &test) const
894 // {
895 // //Find number of nodes
896 // unsigned n_node = nnode();
897
898 // //Get the shape functions
899 // shape(s,psi);
900
901 // //Set the test functions to be the same as the shape functions
902 // for(unsigned i=0;i<n_node;i++)
903 // {
904 // test[i] = psi[i];
905 // }
906
907 // //Return the value of the jacobian
908 // return J_eulerian(s);
909 // }
910
911
912 // /// Function to compute the shape and test functions and to return
913 // /// the Jacobian of mapping between local and global (Eulerian)
914 // /// coordinates
915 // inline double shape_and_test_at_knot(const unsigned &ipt,
916 // Shape &psi,
917 // Shape &test) const
918 // {
919 // //Find number of nodes
920 // unsigned n_node = nnode();
921
922 // //Get the shape functions
923 // shape_at_knot(ipt,psi);
924
925 // //Set the test functions to be the same as the shape functions
926 // for(unsigned i=0;i<n_node;i++)
927 // {
928 // test[i] = psi[i];
929 // }
930
931 // //Return the value of the jacobian
932 // return J_eulerian_at_knot(ipt);
933 // }
934
935 // /// Function to calculate the prescribed beta at a given spatial
936 // /// position
937 // void get_beta(const Vector<double>& x, double& beta)
938 // {
939 // //If the function pointer is zero return zero
940 // if(Beta_fct_pt == 0)
941 // {
942 // beta = 0.0;
943 // }
944 // //Otherwise call the function
945 // else
946 // {
947 // (*Beta_fct_pt)(x,beta);
948 // }
949 // }
950
951 // /// Function to calculate the prescribed alpha at a given spatial
952 // /// position
953 // void get_alpha(const Vector<double>& x, double& alpha)
954 // {
955 // //If the function pointer is zero return zero
956 // if(Alpha_fct_pt == 0)
957 // {
958 // alpha = 0.0;
959 // }
960 // //Otherwise call the function
961 // else
962 // {
963 // (*Alpha_fct_pt)(x,alpha);
964 // }
965 // }
966
967 // private:
968
969
970 // /// Add the element's contribution to its residual vector.
971 // /// flag=1(or 0): do (or don't) compute the Jacobian as well.
972 // void fill_in_generic_residual_contribution_axi_adv_diff_flux(
973 // Vector<double> &residuals, DenseMatrix<double> &jacobian,
974 // unsigned flag);
975
976
977 // /// Function pointer to the (global) prescribed-beta function
978 // AxisymAdvectionDiffusionPrescribedBetaFctPt Beta_fct_pt;
979
980 // /// Function pointer to the (global) prescribed-alpha function
981 // AxisymAdvectionDiffusionPrescribedAlphaFctPt Alpha_fct_pt;
982
983 // /// The index at which the unknown is stored at the nodes
984 // unsigned U_index_adv_diff;
985
986
987 // };//End class AxisymAdvectionDiffusionFluxElement
988
989
990 // ///////////////////////////////////////////////////////////////////////
991 // ///////////////////////////////////////////////////////////////////////
992 // ///////////////////////////////////////////////////////////////////////
993
994
995 // //===========================================================================
996 // /// Constructor, takes the pointer to the "bulk" element and the
997 // index
998 // /// of the face to be created
999 // //===========================================================================
1000 // template<class ELEMENT>
1001 // AxisymAdvectionDiffusionFluxElement<ELEMENT>::
1002 // AxisymAdvectionDiffusionFluxElement(FiniteElement* const &bulk_el_pt,
1003 // const int &face_index) :
1004 // FaceGeometry<ELEMENT>(), FaceElement()
1005
1006 // {
1007 // //Let the bulk element build the FaceElement, i.e. setup the pointers
1008 // //to its nodes (by referring to the appropriate nodes in the bulk
1009 // //element), etc.
1010 // bulk_el_pt->build_face_element(face_index,this);
1011
1012 // #ifdef PARANOID
1013 // {
1014 // //Check that the element is not a refineable 3d element
1015 // ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
1016 // //If it's three-d
1017 // if(elem_pt->dim()==3)
1018 // {
1019 // //Is it refineable
1020 // RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
1021 // if(ref_el_pt!=0)
1022 // {
1023 // if (this->has_hanging_nodes())
1024 // {
1025 // throw OomphLibError(
1026 // "This flux element will not work correctly if nodes are
1027 // hanging\n", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1028 // }
1029 // }
1030 // }
1031 // }
1032 // #endif
1033
1034 // //Initialise the prescribed-beta function pointer to zero
1035 // Beta_fct_pt = 0;
1036
1037 // //Set up U_index_adv_diff. Initialise to zero, which probably won't change
1038 // //in most cases, oh well, the price we pay for generality
1039 // U_index_adv_diff = 0;
1040
1041 // //Cast to the appropriate AdvectionDiffusionEquation so that we can
1042 // //find the index at which the variable is stored
1043 // //We assume that the dimension of the full problem is the same
1044 // //as the dimension of the node, if this is not the case you will have
1045 // //to write custom elements, sorry
1046 // AxisymAdvectionDiffusionEquations* eqn_pt =
1047 // dynamic_cast<AxisymAdvectionDiffusionEquations*>(bulk_el_pt);
1048
1049 // //If the cast has failed die
1050 // if(eqn_pt==0)
1051 // {
1052 // std::string error_string =
1053 // "Bulk element must inherit from AxisymAdvectionDiffusionEquations.";
1054 // error_string +=
1055 // "Nodes are two dimensional, but cannot cast the bulk element to\n";
1056 // error_string += "AxisymAdvectionDiffusionEquations<2>\n.";
1057 // error_string +=
1058 // "If you desire this functionality, you must implement it yourself\n";
1059
1060 // throw OomphLibError(
1061 // error_string,
1062 // OOMPH_CURRENT_FUNCTION,
1063 // OOMPH_EXCEPTION_LOCATION);
1064 // }
1065 // else
1066 // {
1067 // //Read the index from the (cast) bulk element.
1068 // U_index_adv_diff = eqn_pt->u_index_axi_adv_diff();
1069 // }
1070
1071
1072 // }
1073
1074
1075 // //===========================================================================
1076 // /// Compute the element's residual vector and the (zero) Jacobian
1077 // /// matrix for the Robin boundary condition:
1078 // /// \f[ \Delta u \cdot \mathbf{n} + \alpha (\mathbf{x}) = \beta (\mathbf{x}) \f]
1079 // //===========================================================================
1080 // template<class ELEMENT>
1081 // void AxisymAdvectionDiffusionFluxElement<ELEMENT>::
1082 // fill_in_generic_residual_contribution_axi_adv_diff_flux(
1083 // Vector<double> &residuals,
1084 // DenseMatrix<double> &jacobian,
1085 // unsigned flag)
1086 // {
1087 // //Find out how many nodes there are
1088 // const unsigned n_node = nnode();
1089
1090 // // Locally cache the index at which the variable is stored
1091 // const unsigned u_index_axi_adv_diff = U_index_adv_diff;
1092
1093 // //Set up memory for the shape and test functions
1094 // Shape psif(n_node), testf(n_node);
1095
1096 // //Set the value of n_intpt
1097 // const unsigned n_intpt = integral_pt()->nweight();
1098
1099 // //Set the Vector to hold local coordinates
1100 // Vector<double> s(1);
1101
1102 // //Integers used to store the local equation number and local unknown
1103 // //indices for the residuals and jacobians
1104 // int local_eqn=0, local_unknown=0;
1105
1106 // //Loop over the integration points
1107 // //--------------------------------
1108 // for(unsigned ipt=0;ipt<n_intpt;ipt++)
1109 // {
1110
1111 // //Assign values of s
1112 // for(unsigned i=0;i<1;i++)
1113 // {
1114 // s[i] = integral_pt()->knot(ipt,i);
1115 // }
1116
1117 // //Get the integral weight
1118 // double w = integral_pt()->weight(ipt);
1119
1120 // //Find the shape and test functions and return the Jacobian
1121 // //of the mapping
1122 // double J = shape_and_test(s,psif,testf);
1123
1124 // //Premultiply the weights and the Jacobian
1125 // double W = w*J;
1126
1127 // //Calculate local values of the solution and its derivatives
1128 // //Allocate
1129 // double interpolated_u=0.0;
1130 // Vector<double> interpolated_x(2,0.0);
1131
1132 // //Calculate position
1133 // for(unsigned l=0;l<n_node;l++)
1134 // {
1135 // //Get the value at the node
1136 // double u_value = raw_nodal_value(l,u_index_axi_adv_diff);
1137 // interpolated_u += u_value*psif(l);
1138 // //Loop over coordinate direction
1139 // for(unsigned i=0;i<2;i++)
1140 // {
1141 // interpolated_x[i] += nodal_position(l,i)*psif(l);
1142 // }
1143 // }
1144
1145 // //Get the imposed beta (beta=flux when alpha=0.0)
1146 // double beta;
1147 // get_beta(interpolated_x,beta);
1148
1149 // //Get the imposed alpha
1150 // double alpha;
1151 // get_alpha(interpolated_x,alpha);
1152
1153 // //calculate the area weighting dS = r^{2} sin theta dr dtheta
1154 // // r = x[0] and theta = x[1]
1155 // double dS = interpolated_x[0]*interpolated_x[0]*sin(interpolated_x[1]);
1156
1157 // //Now add to the appropriate equations
1158
1159 // //Loop over the test functions
1160 // for(unsigned l=0;l<n_node;l++)
1161 // {
1162 // //Set the local equation number
1163 // local_eqn = nodal_local_eqn(l,u_index_axi_adv_diff);
1164 // /*IF it's not a boundary condition*/
1165 // if(local_eqn >= 0)
1166 // {
1167 // //Add the prescribed beta terms
1168 // residuals[local_eqn] -= dS*(beta - alpha*interpolated_u)*testf(l)*W;
1169
1170 // //Calculate the Jacobian
1171 // //----------------------
1172 // if ( (flag) && (alpha!=0.0) )
1173 // {
1174 // //Loop over the velocity shape functions again
1175 // for(unsigned l2=0;l2<n_node;l2++)
1176 // {
1177 // //Set the number of the unknown
1178 // local_unknown = nodal_local_eqn(l2,u_index_axi_adv_diff);
1179
1180 // //If at a non-zero degree of freedom add in the entry
1181 // if(local_unknown >= 0)
1182 // {
1183 // jacobian(local_eqn,local_unknown) +=
1184 // dS*alpha*psif[l2]*testf[l]*W;
1185 // }
1186 // }
1187 // }
1188
1189 // }
1190 // } //end loop over test functions
1191
1192 // } //end loop over integration points
1193
1194 // }//end fill_in_generic_residual_contribution_adv_diff_flux
1195
1196 // } //End AxisymAdvectionDiffusionFluxElement class
1197} // namespace oomph
1198#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 all elements that solve the Advection Diffusion equations in a cylindrical polar coordina...
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
AxisymAdvectionDiffusionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
virtual void fill_in_generic_residual_contribution_axi_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.
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: [du/dr,du/dz].
AxisymAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
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(* AxisymAdvectionDiffusionSourceFctPt)(const Vector< double > &x, double &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
const double & d() const
Peclet number multiplied by Strouhal number.
void(* AxisymAdvectionDiffusionWindFctPt)(const Vector< double > &x, Vector< double > &wind)
Function pointer to wind function fct(x,w(x)) – x is a Vector!
virtual void get_wind_axi_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...
virtual unsigned u_index_axi_adv_diff() const
Broken assignment operator.
double du_dt_axi_adv_diff(const unsigned &n) const
du/dt at local node n. Uses suitably interpolated value for hanging nodes.
const double & pe_st() const
Peclet number multiplied by Strouhal number.
double *& d_pt()
Pointer to Peclet number multipled by Strouha number.
virtual double dshape_and_dtest_eulerian_at_knot_axi_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(FILE *file_pt)
C_style output with default number of plot points.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,z,u_exact at nplot^2 plot points.
AxisymAdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
virtual void dinterpolated_u_axi_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....
AxisymAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
static double Default_peclet_number
Static default value for the Peclet number.
double interpolated_u_axi_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
virtual double dshape_and_dtest_eulerian_axi_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 compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
AxisymAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
AxisymAdvectionDiffusionEquations()
Constructor: Initialise the Source_fct_pt and Wind_fct_pt to null and set (pointer to) Peclet number ...
AxisymAdvectionDiffusionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
static double Default_diffusion_parameter
Static default value for the Peclet number.
bool ALE_is_disabled
Boolean flag to indicate whether AlE formulation is disable.
virtual void get_source_axi_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...
AxisymAdvectionDiffusionEquations(const AxisymAdvectionDiffusionEquations &dummy)=delete
Broken copy constructor.
double * D_pt
Pointer to global Diffusion parameter.
void output(std::ostream &outfile)
Output with default number of plot points.
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative....
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one,...
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
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2862
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 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...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
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
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
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
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....
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
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.
QAxisymAdvectionDiffusionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,z,u at n_plot^2 plot points.
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.
double dshape_and_dtest_eulerian_at_knot_axi_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....
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,z,u at n_plot^2 plot points.
QAxisymAdvectionDiffusionElement(const QAxisymAdvectionDiffusionElement< NNODE_1D > &dummy)=delete
Broken copy constructor.
double dshape_and_dtest_eulerian_axi_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.
void output(FILE *file_pt)
C-style output function: r,z,u.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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.
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types).
//////////////////////////////////////////////////////////////////// ////////////////////////////////...