gen_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
27#ifndef OOMPH_GEN_AXISYM_ADV_DIFF_ELEMENTS_HEADER
28#define OOMPH_GEN_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/oomph_utilities.h"
40
41namespace oomph
42{
43 //=============================================================
44 // DOXYERROR: the maths below is wrong somehow but I have no idea what it's
45 // supposed to be doing!
46 /// A class for all elements that solve the Advection
47 /// Diffusion equations in conservative form using isoparametric elements
48 /// in a cylindrical polar coordinate system.
49 /// \mbox{\boldmath$\nabla\cdot$} \left(
50 /// Pe \mbox{\boldmath$w$}(\mbox{\boldmath$x$}) u
51 /// - D(\mbox{\boldmath$x$)\mbox{\boldmath$\nabla$} u\right)
52 /// = f(\mbox{\boldmath$x$})
53 /// This contains the generic maths. Shape functions, geometric
54 /// mapping etc. must get implemented in derived class.
55 //=============================================================
57 : public virtual FiniteElement
58 {
59 public:
60 /// Function pointer to source function fct(x,f(x)) --
61 /// x is a Vector!
62 typedef void (*GeneralisedAxisymAdvectionDiffusionSourceFctPt)(
63 const Vector<double>& x, double& f);
64
65 /// Function pointer to wind function fct(x,w(x)) --
66 /// x is a Vector!
67 typedef void (*GeneralisedAxisymAdvectionDiffusionWindFctPt)(
69
70
71 /// Function pointer to a diffusivity function
74
75 /// Constructor: Initialise the Source_fct_pt and Wind_fct_pt
76 /// to null and set (pointer to) Peclet number to default
78 : Source_fct_pt(0),
83 {
84 // Set Peclet number to default
86 // Set Peclet Strouhal number to default
88 }
89
90 /// Broken copy constructor
93
94 /// Broken assignment operator
95 // Commented out broken assignment operator because this can lead to a
96 // conflict warning when used in the virtual inheritence hierarchy.
97 // Essentially the compiler doesn't realise that two separate
98 // implementations of the broken function are the same and so, quite
99 // rightly, it shouts.
100 /*void operator=(const GeneralisedAxisymAdvectionDiffusionEquations&) =
101 * delete;*/
102
103 /// Return the index at which the unknown value
104 /// is stored. The default value, 0, is appropriate for single-physics
105 /// problems, when there is only one variable, the value that satisfies
106 /// the advection-diffusion equation.
107 /// In derived multi-physics elements, this function should be overloaded
108 /// to reflect the chosen storage scheme. Note that these equations require
109 /// that the unknown is always stored at the same index at each node.
110 virtual inline unsigned u_index_cons_axisym_adv_diff() const
111 {
112 return 0;
113 }
114
115 /// du/dt at local node n.
116 /// Uses suitably interpolated value for hanging nodes.
117 double du_dt_cons_axisym_adv_diff(const unsigned& n) const
118 {
119 // Get the data's timestepper
121
122 // Initialise dudt
123 double dudt = 0.0;
124 // Loop over the timesteps, if there is a non Steady timestepper
126 {
127 // Find the index at which the variable is stored
128 const unsigned u_nodal_index = u_index_cons_axisym_adv_diff();
129
130 // Number of timsteps (past & present)
131 const unsigned n_time = time_stepper_pt->ntstorage();
132
133 for (unsigned t = 0; t < n_time; t++)
134 {
135 dudt +=
136 time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
137 }
138 }
139 return dudt;
140 }
141
142 /// Disable ALE, i.e. assert the mesh is not moving -- you do this
143 /// at your own risk!
145 {
146 ALE_is_disabled = true;
147 }
148
149
150 /// (Re-)enable ALE, i.e. take possible mesh motion into account
151 /// when evaluating the time-derivative. Note: By default, ALE is
152 /// enabled, at the expense of possibly creating unnecessary work
153 /// in problems where the mesh is, in fact, stationary.
155 {
156 ALE_is_disabled = false;
157 }
158
159
160 /// Output with default number of plot points
161 void output(std::ostream& outfile)
162 {
163 unsigned nplot = 5;
164 output(outfile, nplot);
165 }
166
167 /// Output FE representation of soln: r,z,u at
168 /// nplot^2 plot points
169 void output(std::ostream& outfile, const unsigned& nplot);
170
171 /// C_style output with default number of plot points
172 void output(FILE* file_pt)
173 {
174 unsigned n_plot = 5;
175 output(file_pt, n_plot);
176 }
177
178 /// C-style output FE representation of soln: r,z,u at
179 /// n_plot^2 plot points
180 void output(FILE* file_pt, const unsigned& n_plot);
181
182
183 /// Output exact soln: r,z,u_exact at nplot^2 plot points
184 void output_fct(std::ostream& outfile,
185 const unsigned& nplot,
187
188 /// Output exact soln: r,z,,u_exact at
189 /// nplot^2 plot points (dummy time-dependent version to
190 /// keep intel compiler happy)
191 virtual void output_fct(
192 std::ostream& outfile,
193 const unsigned& nplot,
194 const double& time,
196 {
197 throw OomphLibError("There is no time-dependent output_fct() for "
198 "Advection Diffusion elements",
199 OOMPH_CURRENT_FUNCTION,
200 OOMPH_EXCEPTION_LOCATION);
201 }
202
203
204 /// Get error against and norm of exact solution
205 void compute_error(std::ostream& outfile,
207 double& error,
208 double& norm);
209
210
211 /// Dummy, time dependent error checker
212 void compute_error(std::ostream& outfile,
214 const double& time,
215 double& error,
216 double& norm)
217 {
218 throw OomphLibError(
219 "No time-dependent compute_error() for Advection Diffusion elements",
220 OOMPH_CURRENT_FUNCTION,
221 OOMPH_EXCEPTION_LOCATION);
222 }
223
224 /// Integrate the concentration over the element
225 double integrate_u();
226
227
228 /// Access function: Pointer to source function
229 GeneralisedAxisymAdvectionDiffusionSourceFctPt& source_fct_pt()
230 {
231 return Source_fct_pt;
232 }
233
234
235 /// Access function: Pointer to source function. Const version
236 GeneralisedAxisymAdvectionDiffusionSourceFctPt source_fct_pt() const
237 {
238 return Source_fct_pt;
239 }
240
241
242 /// Access function: Pointer to wind function
243 GeneralisedAxisymAdvectionDiffusionWindFctPt& wind_fct_pt()
244 {
245 return Wind_fct_pt;
246 }
247
248
249 /// Access function: Pointer to wind function. Const version
250 GeneralisedAxisymAdvectionDiffusionWindFctPt wind_fct_pt() const
251 {
252 return Wind_fct_pt;
253 }
254
255
256 /// Access function: Pointer to additional (conservative) wind function
257 GeneralisedAxisymAdvectionDiffusionWindFctPt& conserved_wind_fct_pt()
258 {
260 }
261
262
263 /// Access function: Pointer to additional (conservative)
264 /// wind function.
265 /// Const version
266 GeneralisedAxisymAdvectionDiffusionWindFctPt conserved_wind_fct_pt() const
267 {
269 }
270
271 /// Access function: Pointer to diffusion function
273 {
274 return Diff_fct_pt;
275 }
276
277 /// Access function: Pointer to diffusion function. Const version
279 {
280 return Diff_fct_pt;
281 }
282
283 /// Peclet number
284 const double& pe() const
285 {
286 return *Pe_pt;
287 }
288
289 /// Pointer to Peclet number
290 double*& pe_pt()
291 {
292 return Pe_pt;
293 }
294
295 /// Peclet number multiplied by Strouhal number
296 const double& pe_st() const
297 {
298 return *PeSt_pt;
299 }
300
301 /// Pointer to Peclet number multipled by Strouha number
302 double*& pe_st_pt()
303 {
304 return PeSt_pt;
305 }
306
307 /// Get source term at (Eulerian) position x. This function is
308 /// virtual to allow overloading in multi-physics problems where
309 /// the strength of the source function might be determined by
310 /// another system of equations
311 inline virtual void get_source_cons_axisym_adv_diff(const unsigned& ipt,
312 const Vector<double>& x,
313 double& source) const
314 {
315 // If no source function has been set, return zero
316 if (Source_fct_pt == 0)
317 {
318 source = 0.0;
319 }
320 else
321 {
322 // Get source strength
323 (*Source_fct_pt)(x, source);
324 }
325 }
326
327 /// Get wind at (Eulerian) position x and/or local coordinate s.
328 /// This function is
329 /// virtual to allow overloading in multi-physics problems where
330 /// the wind function might be determined by
331 /// another system of equations
332 inline virtual void get_wind_cons_axisym_adv_diff(
333 const unsigned& ipt,
334 const Vector<double>& s,
335 const Vector<double>& x,
336 Vector<double>& wind) const
337 {
338 // If no wind function has been set, return zero
339 // There are three components of the wind, but only two matter
340 if (Wind_fct_pt == 0)
341 {
342 for (unsigned i = 0; i < 3; i++)
343 {
344 wind[i] = 0.0;
345 }
346 }
347 else
348 {
349 // Get wind
350 (*Wind_fct_pt)(x, wind);
351 }
352 }
353
354
355 /// Get additional (conservative)
356 /// wind at (Eulerian) position x and/or local coordinate s.
357 /// This function is
358 /// virtual to allow overloading in multi-physics problems where
359 /// the wind function might be determined by
360 /// another system of equations
362 const unsigned& ipt,
363 const Vector<double>& s,
364 const Vector<double>& x,
365 Vector<double>& wind) const
366 {
367 // If no wind function has been set, return zero
368 if (Conserved_wind_fct_pt == 0)
369 {
370 for (unsigned i = 0; i < 3; i++)
371 {
372 wind[i] = 0.0;
373 }
374 }
375 else
376 {
377 // Get wind
378 (*Conserved_wind_fct_pt)(x, wind);
379 }
380 }
381
382
383 /// Get diffusivity tensor at (Eulerian) position
384 /// x and/or local coordinate s.
385 /// This function is
386 /// virtual to allow overloading in multi-physics problems where
387 /// the wind function might be determined by
388 /// another system of equations
389 inline virtual void get_diff_cons_axisym_adv_diff(
390 const unsigned& ipt,
391 const Vector<double>& s,
392 const Vector<double>& x,
393 DenseMatrix<double>& D) const
394 {
395 // If no wind function has been set, return identity
396 // Again three components, but not all of them matter
397 // Those in the theta direction are ignored
398 if (Diff_fct_pt == 0)
399 {
400 for (unsigned i = 0; i < 3; i++)
401 {
402 for (unsigned j = 0; j < 3; j++)
403 {
404 if (i == j)
405 {
406 D(i, j) = 1.0;
407 }
408 else
409 {
410 D(i, j) = 0.0;
411 }
412 }
413 }
414 }
415 else
416 {
417 // Get diffusivity tensor
418 (*Diff_fct_pt)(x, D);
419 }
420 }
421
422
423 /// Get flux: \f$\mbox{flux}[i] = \mbox{d}u / \mbox{d}x_i \f$
424 void get_flux(const Vector<double>& s, Vector<double>& flux) const
425 {
426 // Find out how many nodes there are in the element
427 const unsigned n_node = this->nnode();
428
429 // Get the nodal index at which the unknown is stored
430 const unsigned u_nodal_index = this->u_index_cons_axisym_adv_diff();
431
432 // Set up memory for the shape and test functions
433 Shape psi(n_node);
434 DShape dpsidx(n_node, 2);
435
436 // Call the derivatives of the shape and test functions
437 dshape_eulerian(s, psi, dpsidx);
438
439 // Initialise to zero
440 for (unsigned j = 0; j < 2; j++)
441 {
442 flux[j] = 0.0;
443 }
444
445 // Loop over nodes
446 for (unsigned l = 0; l < n_node; l++)
447 {
448 const double u_value = this->nodal_value(l, u_nodal_index);
449 // Add in the derivative directions
450 flux[0] += u_value * dpsidx(l, 0);
451 flux[1] += u_value * dpsidx(l, 1);
452 }
453 }
454
455 /// Get flux: \f$\mbox{flux}[i] = \mbox{d}u / \mbox{d}x_i \f$
457 Vector<double>& total_flux) const
458 {
459 // Find out how many nodes there are in the element
460 const unsigned n_node = nnode();
461
462 // Get the nodal index at which the unknown is stored
463 const unsigned u_nodal_index = u_index_cons_axisym_adv_diff();
464
465 // Set up memory for the shape and test functions
466 Shape psi(n_node);
467 DShape dpsidx(n_node, 2);
468
469 // Call the derivatives of the shape and test functions
470 dshape_eulerian(s, psi, dpsidx);
471
472 // Storage for the Eulerian position
474 // Storage for the concentration
475 double interpolated_u = 0.0;
476 // Storage for the derivatives of the concentration
477 Vector<double> interpolated_dudx(2, 0.0);
478
479 // Loop over nodes
480 for (unsigned l = 0; l < n_node; l++)
481 {
482 // Get the value at the node
483 const double u_value = this->nodal_value(l, u_nodal_index);
484 interpolated_u += u_value * psi(l);
485 // Loop over directions
486 for (unsigned j = 0; j < 2; j++)
487 {
488 interpolated_x[j] += this->nodal_position(l, j) * psi(l);
489 interpolated_dudx[j] += u_value * dpsidx(l, j);
490 }
491 }
492
493 // Dummy integration point
494 unsigned ipt = 0;
495
496 // Get the conserved wind (non-divergence free)
497 Vector<double> conserved_wind(3);
499 ipt, s, interpolated_x, conserved_wind);
500
501 // Get diffusivity tensor
504
505 // Calculate the total flux made up of the diffusive flux
506 // and the conserved wind only bother with the
507 // first two components in each case becuase there can be no
508 // variation in the azimuthal direction
509 for (unsigned i = 0; i < 2; i++)
510 {
511 total_flux[i] = 0.0;
512 for (unsigned j = 0; j < 2; j++)
513 {
514 total_flux[i] += D(i, j) * interpolated_dudx[j];
515 }
516 total_flux[i] -= conserved_wind[i] * interpolated_u;
517 }
518 }
519
520
521 /// Add the element's contribution to its residual vector (wrapper)
523 {
524 // Call the generic residuals function with flag set to 0 and using
525 // a dummy matrix
527 residuals,
530 0);
531 }
532
533
534 /// Add the element's contribution to its residual vector and
535 /// the element Jacobian matrix (wrapper)
537 DenseMatrix<double>& jacobian)
538 {
539 // Call the generic routine with the flag set to 1
541 residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
542 }
543
544
545 /// Add the element's contribution to its residuals vector,
546 /// jacobian matrix and mass matrix
548 Vector<double>& residuals,
549 DenseMatrix<double>& jacobian,
550 DenseMatrix<double>& mass_matrix)
551 {
552 // Call the generic routine with the flag set to 2
554 residuals, jacobian, mass_matrix, 2);
555 }
556
557
558 /// Return FE representation of function value u(s) at local coordinate s
560 const Vector<double>& s) const
561 {
562 // Find number of nodes
563 unsigned n_node = nnode();
564
565 // Get the nodal index at which the unknown is stored
566 unsigned u_nodal_index = u_index_cons_axisym_adv_diff();
567
568 // Local shape function
569 Shape psi(n_node);
570
571 // Find values of shape function
572 shape(s, psi);
573
574 // Initialise value of u
575 double interpolated_u = 0.0;
576
577 // Loop over the local nodes and sum
578 for (unsigned l = 0; l < n_node; l++)
579 {
580 interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
581 }
582
583 return (interpolated_u);
584 }
585
586
587 /// Self-test: Return 0 for OK
588 unsigned self_test();
589
590 protected:
591 /// Shape/test functions and derivs w.r.t. to global coords at
592 /// local coord. s; return Jacobian of mapping
594 const Vector<double>& s,
595 Shape& psi,
596 DShape& dpsidx,
597 Shape& test,
598 DShape& dtestdx) const = 0;
599
600 /// Shape/test functions and derivs w.r.t. to global coords at
601 /// integration point ipt; return Jacobian of mapping
603 const unsigned& ipt,
604 Shape& psi,
605 DShape& dpsidx,
606 Shape& test,
607 DShape& dtestdx) const = 0;
608
609 /// Add the element's contribution to its residual vector only
610 /// (if flag=and/or element Jacobian matrix
612 Vector<double>& residuals,
613 DenseMatrix<double>& jacobian,
614 DenseMatrix<double>& mass_matrix,
615 unsigned flag);
616
617 /// Pointer to global Peclet number
618 double* Pe_pt;
619
620 /// Pointer to global Peclet number multiplied by Strouhal number
621 double* PeSt_pt;
622
623 /// Pointer to source function:
624 GeneralisedAxisymAdvectionDiffusionSourceFctPt Source_fct_pt;
625
626 /// Pointer to wind function:
627 GeneralisedAxisymAdvectionDiffusionWindFctPt Wind_fct_pt;
628
629 /// Pointer to additional (conservative) wind function:
630 GeneralisedAxisymAdvectionDiffusionWindFctPt Conserved_wind_fct_pt;
631
632 /// Pointer to diffusivity funciton
634
635 /// Boolean flag to indicate if ALE formulation is disabled when
636 /// time-derivatives are computed. Only set to false if you're sure
637 /// that the mesh is stationary.
639
640 private:
641 /// Static default value for the Peclet number
643 };
644
645
646 /// ////////////////////////////////////////////////////////////////////////
647 /// ////////////////////////////////////////////////////////////////////////
648 /// ////////////////////////////////////////////////////////////////////////
649
650
651 //======================================================================
652 /// QGeneralisedAxisymAdvectionDiffusionElement elements are
653 /// linear/quadrilateral/brick-shaped Advection Diffusion elements with
654 /// isoparametric interpolation for the function.
655 //======================================================================
656 template<unsigned NNODE_1D>
658 : public virtual QElement<2, NNODE_1D>,
660 {
661 private:
662 /// Static array of ints to hold number of variables at
663 /// nodes: Initial_Nvalue[n]
664 static const unsigned Initial_Nvalue;
665
666 public:
667 /// Constructor: Call constructors for QElement and
668 /// Advection Diffusion equations
671 {
672 }
673
674 /// Broken copy constructor
677 delete;
678
679 /// Broken assignment operator
680 /*void operator=(const
681 QGeneralisedAxisymAdvectionDiffusionElement<NNODE_1D>&) = delete;*/
682
683 /// Required # of `values' (pinned or dofs)
684 /// at node n
685 inline unsigned required_nvalue(const unsigned& n) const
686 {
687 return Initial_Nvalue;
688 }
689
690 /// Output function:
691 /// r,z,u
692 void output(std::ostream& outfile)
693 {
695 }
696
697 /// Output function:
698 /// r,z,u at n_plot^2 plot points
699 void output(std::ostream& outfile, const unsigned& n_plot)
700 {
702 }
703
704
705 /// C-style output function:
706 /// r,z,u
707 void output(FILE* file_pt)
708 {
710 }
711
712 /// C-style output function:
713 /// r,z,u at n_plot^2 plot points
714 void output(FILE* file_pt, const unsigned& n_plot)
715 {
717 }
718
719 /// Output function for an exact solution:
720 /// r,z,u_exact at n_plot^2 plot points
721 void output_fct(std::ostream& outfile,
722 const unsigned& n_plot,
724 {
726 outfile, n_plot, exact_soln_pt);
727 }
728
729
730 /// Output function for a time-dependent exact solution.
731 /// r,z,u_exact at n_plot^2 plot points
732 /// (Calls the steady version)
733 void output_fct(std::ostream& outfile,
734 const unsigned& n_plot,
735 const double& time,
737 {
739 outfile, n_plot, time, exact_soln_pt);
740 }
741
742
743 protected:
744 /// Shape, test functions & derivs. w.r.t. to global coords. Return
745 /// Jacobian.
747 const Vector<double>& s,
748 Shape& psi,
749 DShape& dpsidx,
750 Shape& test,
751 DShape& dtestdx) const;
752
753 /// Shape, test functions & derivs. w.r.t. to global coords. at
754 /// integration point ipt. Return Jacobian.
756 const unsigned& ipt,
757 Shape& psi,
758 DShape& dpsidx,
759 Shape& test,
760 DShape& dtestdx) const;
761 };
762
763 // Inline functions:
764
765
766 //======================================================================
767 /// Define the shape functions and test functions and derivatives
768 /// w.r.t. global coordinates and return Jacobian of mapping.
769 ///
770 /// Galerkin: Test functions = shape functions
771 //======================================================================
772 template<unsigned NNODE_1D>
775 Shape& psi,
776 DShape& dpsidx,
777 Shape& test,
778 DShape& dtestdx) const
779 {
780 // Call the geometrical shape functions and derivatives
781 double J = this->dshape_eulerian(s, psi, dpsidx);
782
783 // Loop over the test functions and derivatives and set them equal to the
784 // shape functions
785 for (unsigned i = 0; i < NNODE_1D; i++)
786 {
787 test[i] = psi[i];
788 for (unsigned j = 0; j < 2; j++)
789 {
790 dtestdx(i, j) = dpsidx(i, j);
791 }
792 }
793
794 // Return the jacobian
795 return J;
796 }
797
798
799 //======================================================================
800 /// Define the shape functions and test functions and derivatives
801 /// w.r.t. global coordinates and return Jacobian of mapping.
802 ///
803 /// Galerkin: Test functions = shape functions
804 //======================================================================
805 template<unsigned NNODE_1D>
808 const unsigned& ipt,
809 Shape& psi,
810 DShape& dpsidx,
811 Shape& test,
812 DShape& dtestdx) const
813 {
814 // Call the geometrical shape functions and derivatives
815 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
816
817 // Set the test functions equal to the shape functions (pointer copy)
818 test = psi;
819 dtestdx = dpsidx;
820
821 // Return the jacobian
822 return J;
823 }
824
825
826 /// /////////////////////////////////////////////////////////////////////
827 /// /////////////////////////////////////////////////////////////////////
828 /// /////////////////////////////////////////////////////////////////////
829
830
831 //=======================================================================
832 /// Face geometry for the QGeneralisedAxisymAdvectionDiffusionElement
833 /// elements: The spatial dimension of the face elements is one lower than
834 /// that of the bulk element but they have the same number of points along
835 /// their 1D edges.
836 //=======================================================================
837 template<unsigned NNODE_1D>
839 : public virtual QElement<1, NNODE_1D>
840 {
841 public:
842 /// Constructor: Call the constructor for the
843 /// appropriate lower-dimensional QElement
844 FaceGeometry() : QElement<1, NNODE_1D>() {}
845 };
846
847} // namespace oomph
848
849#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
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 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
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
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
const double & pe_st() const
Peclet number multiplied by Strouhal number.
GeneralisedAxisymAdvectionDiffusionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
GeneralisedAxisymAdvectionDiffusionDiffFctPt Diff_fct_pt
Pointer to diffusivity funciton.
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
double interpolated_u_cons_axisym_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
void output(FILE *file_pt)
C_style output with default number of plot points.
double du_dt_cons_axisym_adv_diff(const unsigned &n) const
du/dt at local node n. Uses suitably interpolated value for hanging nodes.
virtual void get_conserved_wind_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Get additional (conservative) wind at (Eulerian) position x and/or local coordinate s....
GeneralisedAxisymAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
virtual double dshape_and_dtest_eulerian_at_knot_cons_axisym_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 ...
GeneralisedAxisymAdvectionDiffusionDiffFctPt diff_fct_pt() const
Access function: Pointer to diffusion function. Const version.
virtual void get_diff_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &D) const
Get diffusivity tensor at (Eulerian) position x and/or local coordinate s. This function is virtual t...
virtual void output_fct(std::ostream &outfile, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,z,,u_exact at nplot^2 plot points (dummy time-dependent version to keep intel co...
void get_total_flux(const Vector< double > &s, Vector< double > &total_flux) const
Get flux: .
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
GeneralisedAxisymAdvectionDiffusionDiffFctPt & diff_fct_pt()
Access function: Pointer to diffusion function.
static double Default_peclet_number
Static default value for the Peclet number.
GeneralisedAxisymAdvectionDiffusionWindFctPt conserved_wind_fct_pt() const
Access function: Pointer to additional (conservative) wind function. Const version.
virtual void get_wind_cons_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 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.
GeneralisedAxisymAdvectionDiffusionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
GeneralisedAxisymAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
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 residuals vector, jacobian matrix and mass matrix.
GeneralisedAxisymAdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
double integrate_u()
Integrate the concentration over the element.
virtual void fill_in_generic_residual_contribution_cons_axisym_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.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
GeneralisedAxisymAdvectionDiffusionWindFctPt & conserved_wind_fct_pt()
Access function: Pointer to additional (conservative) wind function.
Function pointer to wind function Vector< double > & wind
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
virtual double dshape_and_dtest_eulerian_cons_axisym_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.
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 enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative....
virtual void get_source_cons_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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
GeneralisedAxisymAdvectionDiffusionWindFctPt Conserved_wind_fct_pt
Pointer to additional (conservative) wind function:
Function pointer to a diffusivity function typedef void(* GeneralisedAxisymAdvectionDiffusionDiffFctPt)(const Vector< double > &x, DenseMatrix< double > &D)
Broken copy constructor GeneralisedAxisymAdvectionDiffusionEquations(const GeneralisedAxisymAdvectionDiffusionEquations &dummy)=delete
GeneralisedAxisymAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
void output(std::ostream &outfile)
Output with default number of plot points.
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h: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_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: r,z,u_exact at n_plot^2 plot points.
void output(std::ostream &outfile)
Output 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.
double dshape_and_dtest_eulerian_at_knot_cons_axisym_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_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. r,z,u_exact at n_plot^2 plot points (Calls the s...
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
double dshape_and_dtest_eulerian_cons_axisym_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.
QGeneralisedAxisymAdvectionDiffusionElement(const QGeneralisedAxisymAdvectionDiffusionElement< NNODE_1D > &dummy)=delete
Broken copy constructor.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,z,u at n_plot^2 plot points.
QGeneralisedAxisymAdvectionDiffusionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...