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