advection_diffusion_reaction_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_ADV_DIFF_REACT_ELEMENTS_HEADER
28#define OOMPH_ADV_DIFF_REACT_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/projection.h"
38#include "../generic/nodes.h"
39#include "../generic/Qelements.h"
40#include "../generic/oomph_utilities.h"
41
42namespace oomph
43{
44 //=============================================================
45 /// A class for all elements that solve the Advection
46 /// Diffusion Reaction equations using isoparametric elements.
47 /// \f[ \tau_{i} \frac{\partial C_{i}}{\partial t} + w_{j} \frac{\partial C_{i}}{\partial x_{j}} = D_{i}\frac{\partial^2 C_{i}}{\partial x_j^2} = - R_{i}(C_{i}) + f_{i} \f]
48 /// This contains the generic maths. Shape functions, geometric
49 /// mapping etc. must get implemented in derived class.
50 //=============================================================
51 template<unsigned NREAGENT, unsigned DIM>
53 {
54 public:
55 /// Function pointer to source function fct(x,f(x)) --
56 /// x is a Vector!
58 const Vector<double>& x, Vector<double>& f);
59
60 /// Function pointer to reaction terms
63
64 /// Function pointer to derivative of reaction terms
66 const Vector<double>& c, DenseMatrix<double>& dRdC);
67
68
69 /// Function pointer to wind function fct(x,w(x)) --
70 /// x is a Vector!
71 typedef void (*AdvectionDiffusionReactionWindFctPt)(const double& time,
72 const Vector<double>& x,
73 Vector<double>& wind);
74
75
76 /// Constructor: Initialise the Source_fct_pt, Wind_fct_pt,
77 /// Reaction_fct_pt to null and initialise the dimensionless
78 /// timescale and diffusion ratios
80 : Source_fct_pt(0),
81 Wind_fct_pt(0),
84 ALE_is_disabled(false)
85 {
86 // Set diffusion coefficients to default
88 // Set timescales to default
90 }
91
92 /// Broken copy constructor
94 const AdvectionDiffusionReactionEquations& dummy) = delete;
95
96 /// Broken assignment operator
98
99 /// Return the index at which the unknown i-th reagent
100 /// is stored. The default value, i, is appropriate for single-physics
101 /// problems, when there are only i variables, the values that satisfies
102 /// the advection-diffusion-reaction equation.
103 /// In derived multi-physics elements, this function should be overloaded
104 /// to reflect the chosen storage scheme. Note that these equations require
105 /// that the unknown is always stored at the same index at each node.
106 virtual inline unsigned c_index_adv_diff_react(const unsigned& i) const
107 {
108 return i;
109 }
110
111 /// dc_r/dt at local node n.
112 /// Uses suitably interpolated value for hanging nodes.
113 double dc_dt_adv_diff_react(const unsigned& n, const unsigned& r) const
114 {
115 // Get the data's timestepper
117
118 // Initialise dudt
119 double dudt = 0.0;
120 // Loop over the timesteps, if there is a non Steady timestepper
122 {
123 // Find the index at which the variable is stored
124 const unsigned c_nodal_index = c_index_adv_diff_react(r);
125
126 // Number of timsteps (past & present)
127 const unsigned n_time = time_stepper_pt->ntstorage();
128
129 for (unsigned t = 0; t < n_time; t++)
130 {
131 dudt +=
132 time_stepper_pt->weight(1, t) * nodal_value(t, n, c_nodal_index);
133 }
134 }
135 return dudt;
136 }
137
138
139 /// dc/dt at local node n.
140 /// Uses suitably interpolated value for hanging nodes.
141 void dc_dt_adv_diff_react(const unsigned& n, Vector<double>& dc_dt) const
142 {
143 // Get the data's timestepper
145
146 // Initialise to zero
147 for (unsigned r = 0; r < NREAGENT; r++)
148 {
149 dc_dt[r] = 0.0;
150 }
151
152 // Loop over the timesteps, if there is a non Steady timestepper
154 {
155 // Number of timsteps (past & present)
156 const unsigned n_time = time_stepper_pt->ntstorage();
157 // Local storage (cache) for the weights
158 double weight[n_time];
159 for (unsigned t = 0; t < n_time; t++)
160 {
161 weight[t] = time_stepper_pt->weight(1, t);
162 }
163
164 // Loop over the reagents
165 for (unsigned r = 0; r < NREAGENT; r++)
166 {
167 // Find the index at which the variable is stored
168 const unsigned c_nodal_index = c_index_adv_diff_react(r);
169
170 for (unsigned t = 0; t < n_time; t++)
171 {
172 dc_dt[r] += weight[t] * nodal_value(t, n, c_nodal_index);
173 }
174 }
175 }
176 }
177
178 /// Disable ALE, i.e. assert the mesh is not moving -- you do this
179 /// at your own risk!
181 {
182 ALE_is_disabled = true;
183 }
184
185 /// (Re-)enable ALE, i.e. take possible mesh motion into account
186 /// when evaluating the time-derivative. Note: By default, ALE is
187 /// enabled, at the expense of possibly creating unnecessary work
188 /// in problems where the mesh is, in fact, stationary.
190 {
191 ALE_is_disabled = false;
192 }
193
194
195 /// Output with default number of plot points
196 void output(std::ostream& outfile)
197 {
198 unsigned nplot = 5;
199 output(outfile, nplot);
200 }
201
202 /// Output FE representation of soln: x,y,u or x,y,z,u at
203 /// nplot^DIM plot points
204 void output(std::ostream& outfile, const unsigned& nplot);
205
206
207 /// C_style output with default number of plot points
208 void output(FILE* file_pt)
209 {
210 unsigned n_plot = 5;
211 output(file_pt, n_plot);
212 }
213
214 /// C-style output FE representation of soln: x,y,u or x,y,z,u at
215 /// n_plot^DIM plot points
216 void output(FILE* file_pt, const unsigned& n_plot);
217
218
219 /// Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points
220 void output_fct(std::ostream& outfile,
221 const unsigned& nplot,
223
224 /// Output exact soln: x,y,u_exact or x,y,z,u_exact at
225 /// nplot^DIM plot points (dummy time-dependent version to
226 /// keep intel compiler happy)
227 virtual void output_fct(
228 std::ostream& outfile,
229 const unsigned& nplot,
230 const double& time,
232 {
233 throw OomphLibError("There is no time-dependent output_fct() for "
234 "Advection Diffusion elements",
235 OOMPH_CURRENT_FUNCTION,
236 OOMPH_EXCEPTION_LOCATION);
237 }
238
239 /// Compute norm of the solution:
240 /// sum of squares of the L2 norm for each reagent
241 void compute_norm(double& norm);
242
243 /// Get error against and norm of exact solution
244 void compute_error(std::ostream& outfile,
246 double& error,
247 double& norm);
248
249
250 /// Dummy, time dependent error checker
251 void compute_error(std::ostream& outfile,
253 const double& time,
254 double& error,
255 double& norm)
256 {
257 throw OomphLibError(
258 "No time-dependent compute_error() for Advection Diffusion elements",
259 OOMPH_CURRENT_FUNCTION,
260 OOMPH_EXCEPTION_LOCATION);
261 }
262
263
264 /// Access function: Pointer to source function
266 {
267 return Source_fct_pt;
268 }
269
270
271 /// Access function: Pointer to source function. Const version
273 {
274 return Source_fct_pt;
275 }
276
277
278 /// Access function: Pointer to wind function
280 {
281 return Wind_fct_pt;
282 }
283
284 /// Access function: Pointer to reaction function. Const version
286 {
287 return Wind_fct_pt;
288 }
289
290 /// Access function: Pointer to reaction function
292 {
293 return Reaction_fct_pt;
294 }
295
296 /// Access function: Pointer to reaction function. Const version
298 {
299 return Reaction_fct_pt;
300 }
301
302 /// Access function: Pointer to reaction derivatives function
304 {
306 }
307
308 /// Access function: Pointer to reaction function. Const version
310 {
312 }
313
314 /// Vector of diffusion coefficients
315 const Vector<double>& diff() const
316 {
317 return *Diff_pt;
318 }
319
320 /// Pointer to vector of diffusion coefficients
322 {
323 return Diff_pt;
324 }
325
326 /// Vector of dimensionless timescales
327 const Vector<double>& tau() const
328 {
329 return *Tau_pt;
330 }
331
332 /// Pointer to vector of dimensionless timescales
334 {
335 return Tau_pt;
336 }
337
338 /// Get source term at (Eulerian) position x. This function is
339 /// virtual to allow overloading in multi-physics problems where
340 /// the strength of the source function might be determined by
341 /// another system of equations
342 inline virtual void get_source_adv_diff_react(const unsigned& ipt,
343 const Vector<double>& x,
344 Vector<double>& source) const
345 {
346 // If no source function has been set, return zero
347 if (Source_fct_pt == 0)
348 {
349 for (unsigned r = 0; r < NREAGENT; r++)
350 {
351 source[r] = 0.0;
352 }
353 }
354 else
355 {
356 // Get source strength
357 (*Source_fct_pt)(x, source);
358 }
359 }
360
361 /// Get wind at (Eulerian) position x and/or local coordinate s.
362 /// This function is
363 /// virtual to allow overloading in multi-physics problems where
364 /// the wind function might be determined by
365 /// another system of equations
366 inline virtual void get_wind_adv_diff_react(const unsigned& ipt,
367 const Vector<double>& s,
368 const Vector<double>& x,
369 Vector<double>& wind) const
370 {
371 // If no wind function has been set, return zero
372 if (Wind_fct_pt == 0)
373 {
374 for (unsigned i = 0; i < DIM; i++)
375 {
376 wind[i] = 0.0;
377 }
378 }
379 else
380 {
381 // Get continuous time from timestepper of first node
382 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
383
384 // Get wind
385 (*Wind_fct_pt)(time, x, wind);
386 }
387 }
388
389
390 /// Get reaction as a function of the given reagent concentrations
391 /// This function is
392 /// virtual to allow overloading in multi-physics problems where
393 /// the reaction function might be determined by
394 /// another system of equations
395 inline virtual void get_reaction_adv_diff_react(const unsigned& ipt,
396 const Vector<double>& C,
397 Vector<double>& R) const
398 {
399 // If no wind function has been set, return zero
400 if (Reaction_fct_pt == 0)
401 {
402 for (unsigned r = 0; r < NREAGENT; r++)
403 {
404 R[r] = 0.0;
405 }
406 }
407 else
408 {
409 // Get reaction terms
410 (*Reaction_fct_pt)(C, R);
411 }
412 }
413
414 /// Get the derivatives of the reaction terms with respect to the
415 /// concentration variables. If no explicit function pointer is set,
416 /// these will be calculated by finite differences
418 const unsigned& ipt,
419 const Vector<double>& C,
420 DenseMatrix<double>& dRdC) const
421 {
422 // If no reaction pointer set, return zero
423 if (Reaction_fct_pt == 0)
424 {
425 for (unsigned r = 0; r < NREAGENT; r++)
426 {
427 for (unsigned p = 0; p < NREAGENT; p++)
428 {
429 dRdC(r, p) = 0.0;
430 }
431 }
432 }
433 else
434 {
435 // If no function pointer get finite differences
436 if (Reaction_deriv_fct_pt == 0)
437 {
438 // Local copy of the unknowns
439 Vector<double> C_local = C;
440 // Finite differences
441 Vector<double> R(NREAGENT), R_plus(NREAGENT), R_minus(NREAGENT);
442 // Get the initial reaction terms
443 //(*Reaction_fct_pt)(C,R);
445 // Now loop over all the reagents
446 for (unsigned p = 0; p < NREAGENT; p++)
447 {
448 // Store the old value
449 double old_var = C_local[p];
450 // Increment the value
451 C_local[p] += fd_step;
452 // Get the new values
453 (*Reaction_fct_pt)(C_local, R_plus);
454
455 // Reset the value
456 C_local[p] = old_var;
457 // Decrement the value
458 C_local[p] -= fd_step;
459 // Get the new values
460 (*Reaction_fct_pt)(C_local, R_minus);
461
462 // Assemble the column of the jacobian
463 for (unsigned r = 0; r < NREAGENT; r++)
464 {
465 dRdC(r, p) = (R_plus[r] - R_minus[r]) / (2.0 * fd_step);
466 }
467
468 // Reset the value
469 C_local[p] = old_var;
470 }
471 }
472 // Otherwise get the terms from the function
473 else
474 {
475 (*Reaction_deriv_fct_pt)(C, dRdC);
476 }
477 }
478 }
479
480 /// Get flux: \f$\mbox{flux}[DIM r + i] = \mbox{d}C_{r} / \mbox{d}x_i \f$
481 void get_flux(const Vector<double>& s, Vector<double>& flux) const
482 {
483 // Find out how many nodes there are in the element
484 const unsigned n_node = nnode();
485
486 // Set up memory for the shape and test functions
487 Shape psi(n_node);
488 DShape dpsidx(n_node, DIM);
489
490 // Call the derivatives of the shape and test functions
491 dshape_eulerian(s, psi, dpsidx);
492
493 // Initialise to zero
494 for (unsigned j = 0; j < DIM * NREAGENT; j++)
495 {
496 flux[j] = 0.0;
497 }
498
499 // Loop over the reagent terms
500 for (unsigned r = 0; r < NREAGENT; r++)
501 {
502 unsigned c_nodal_index = c_index_adv_diff_react(r);
503
504 // Loop over derivative directions
505 for (unsigned j = 0; j < DIM; j++)
506 {
507 unsigned index = r * DIM + j;
508 // Loop over the nodes
509 for (unsigned l = 0; l < n_node; l++)
510 {
511 flux[index] += nodal_value(l, c_nodal_index) * dpsidx(l, j);
512 }
513 }
514 }
515 }
516
517
518 /// Add the element's contribution to its residual vector (wrapper)
520 {
521 // Call the generic residuals function with flag set to 0 and using
522 // a dummy matrix
524 residuals,
527 0);
528 }
529
530
531 /// Add the element's contribution to its residual vector and
532 /// the element Jacobian matrix (wrapper)
534 DenseMatrix<double>& jacobian)
535 {
536 // Call the generic routine with the flag set to 1
538 residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
539 }
540
541
542 /// Add the element's contribution to its residuals vector,
543 /// jacobian matrix and mass matrix
545 Vector<double>& residuals,
546 DenseMatrix<double>& jacobian,
547 DenseMatrix<double>& mass_matrix)
548 {
549 // Call the generic routine with the flag set to 2
551 residuals, jacobian, mass_matrix, 2);
552 }
553
554
555 /// Return FE representation of function value c_i(s) at local coordinate s
557 const unsigned& i) const
558 {
559 // Find number of nodes
560 unsigned n_node = nnode();
561
562 // Get the nodal index at which the unknown is stored
563 unsigned c_nodal_index = c_index_adv_diff_react(i);
564
565 // Local shape function
566 Shape psi(n_node);
567
568 // Find values of shape function
569 shape(s, psi);
570
571 // Initialise value of u
572 double interpolated_c = 0.0;
573
574 // Loop over the local nodes and sum
575 for (unsigned l = 0; l < n_node; l++)
576 {
577 interpolated_c += nodal_value(l, c_nodal_index) * psi[l];
578 }
579
580 return (interpolated_c);
581 }
582
583
584 /// Self-test: Return 0 for OK
585 unsigned self_test();
586
587 /// Return the integrated reagent concentrations
588 void integrate_reagents(Vector<double>& C) const;
589
590 protected:
591 // Static variable that holds the number of reagents
592 static const unsigned N_reagent;
593
594
595 /// Shape/test functions and derivs w.r.t. to global coords at
596 /// local coord. s; return Jacobian of mapping
598 const Vector<double>& s,
599 Shape& psi,
600 DShape& dpsidx,
601 Shape& test,
602 DShape& dtestdx) const = 0;
603
604 /// Shape/test functions and derivs w.r.t. to global coords at
605 /// integration point ipt; return Jacobian of mapping
607 const unsigned& ipt,
608 Shape& psi,
609 DShape& dpsidx,
610 Shape& test,
611 DShape& dtestdx) const = 0;
612
613 /// Add the element's contribution to its residual vector only
614 /// (if flag=and/or element Jacobian matrix
616 Vector<double>& residuals,
617 DenseMatrix<double>& jacobian,
618 DenseMatrix<double>& mass_matrix,
619 unsigned flag);
620
621 /// Pointer to global diffusion coefficients
623
624 /// Pointer to global timescales
626
627 /// Pointer to source function:
629
630 /// Pointer to wind function:
632
633 /// Pointer to reaction function
635
636 /// Pointer to reaction derivatives
638
639 /// Boolean flag to indicate if ALE formulation is disabled when
640 /// time-derivatives are computed. Only set to true if you're sure
641 /// that the mesh is stationary.
643
644 private:
645 /// Static default value for the dimensionless numbers
647 };
648
649
650 /// ////////////////////////////////////////////////////////////////////////
651 /// ////////////////////////////////////////////////////////////////////////
652 /// ////////////////////////////////////////////////////////////////////////
653
654
655 //======================================================================
656 /// QAdvectionDiffusionReactionElement elements are
657 /// linear/quadrilateral/brick-shaped Advection Diffusion elements with
658 /// isoparametric interpolation for the function.
659 //======================================================================
660 template<unsigned NREAGENT, unsigned DIM, unsigned NNODE_1D>
662 : public virtual QElement<DIM, NNODE_1D>,
663 public virtual AdvectionDiffusionReactionEquations<NREAGENT, DIM>
664 {
665 public:
666 /// Constructor: Call constructors for QElement and
667 /// Advection Diffusion equations
669 : QElement<DIM, NNODE_1D>(),
671 {
672 }
673
674 /// Broken copy constructor
677 dummy) = delete;
678
679 /// Broken assignment operator
682 delete;
683
684 /// Required # of `values' (pinned or dofs)
685 /// at node n
686 inline unsigned required_nvalue(const unsigned& n) const
687 {
688 return NREAGENT;
689 }
690
691 /// Output function:
692 /// x,y,u or x,y,z,u
693 void output(std::ostream& outfile)
694 {
696 }
697
698 /// Output function:
699 /// x,y,u or x,y,z,u at n_plot^DIM plot points
700 void output(std::ostream& outfile, const unsigned& n_plot)
701 {
703 n_plot);
704 }
705
706
707 /// C-style output function:
708 /// x,y,u or x,y,z,u
709 void output(FILE* file_pt)
710 {
712 }
713
714 /// C-style output function:
715 /// x,y,u or x,y,z,u at n_plot^DIM plot points
716 void output(FILE* file_pt, const unsigned& n_plot)
717 {
719 n_plot);
720 }
721
722 /// Output function for an exact solution:
723 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
724 void output_fct(std::ostream& outfile,
725 const unsigned& n_plot,
727 {
729 outfile, n_plot, exact_soln_pt);
730 }
731
732
733 /// Output function for a time-dependent exact solution.
734 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
735 /// (Calls the steady version)
736 void output_fct(std::ostream& outfile,
737 const unsigned& n_plot,
738 const double& time,
740 {
742 outfile, n_plot, time, exact_soln_pt);
743 }
744
745
746 protected:
747 /// Shape, test functions & derivs. w.r.t. to global coords. Return
748 /// Jacobian.
750 const Vector<double>& s,
751 Shape& psi,
752 DShape& dpsidx,
753 Shape& test,
754 DShape& dtestdx) const;
755
756 /// Shape, test functions & derivs. w.r.t. to global coords. at
757 /// integration point ipt. Return Jacobian.
759 const unsigned& ipt,
760 Shape& psi,
761 DShape& dpsidx,
762 Shape& test,
763 DShape& dtestdx) const;
764 };
765
766 // Inline functions:
767
768
769 //======================================================================
770 /// Define the shape functions and test functions and derivatives
771 /// w.r.t. global coordinates and return Jacobian of mapping.
772 ///
773 /// Galerkin: Test functions = shape functions
774 //======================================================================
775 template<unsigned NREAGENT, unsigned DIM, unsigned NNODE_1D>
778 Shape& psi,
779 DShape& dpsidx,
780 Shape& test,
781 DShape& dtestdx) const
782 {
783 // Call the geometrical shape functions and derivatives
784 double J = this->dshape_eulerian(s, psi, dpsidx);
785
786 // Loop over the test functions and derivatives and set them equal to the
787 // shape functions
788 for (unsigned i = 0; i < NNODE_1D; i++)
789 {
790 test[i] = psi[i];
791 for (unsigned j = 0; j < DIM; j++)
792 {
793 dtestdx(i, j) = dpsidx(i, j);
794 }
795 }
796
797 // Return the jacobian
798 return J;
799 }
800
801
802 //======================================================================
803 /// Define the shape functions and test functions and derivatives
804 /// w.r.t. global coordinates and return Jacobian of mapping.
805 ///
806 /// Galerkin: Test functions = shape functions
807 //======================================================================
808 template<unsigned NREAGENT, unsigned DIM, unsigned NNODE_1D>
811 Shape& psi,
812 DShape& dpsidx,
813 Shape& test,
814 DShape& dtestdx) const
815 {
816 // Call the geometrical shape functions and derivatives
817 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
818
819 // Set the test functions equal to the shape functions (pointer copy)
820 test = psi;
821 dtestdx = dpsidx;
822
823 // Return the jacobian
824 return J;
825 }
826
827
828 /// /////////////////////////////////////////////////////////////////////
829 /// /////////////////////////////////////////////////////////////////////
830 /// /////////////////////////////////////////////////////////////////////
831
832
833 //=======================================================================
834 /// Face geometry for the QAdvectionDiffusionReactionElement elements:
835 /// The spatial dimension of the face elements is one lower than that
836 /// of the bulk element but they have the same number of points along
837 /// their 1D edges.
838 //=======================================================================
839 template<unsigned NREAGENT, unsigned DIM, unsigned NNODE_1D>
841 QAdvectionDiffusionReactionElement<NREAGENT, DIM, NNODE_1D>>
842 : public virtual QElement<DIM - 1, NNODE_1D>
843 {
844 public:
845 /// Constructor: Call the constructor for the
846 /// appropriate lower-dimensional QElement
847 FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
848 };
849
850
851 /// /////////////////////////////////////////////////////////////////////
852 /// /////////////////////////////////////////////////////////////////////
853 /// /////////////////////////////////////////////////////////////////////
854
855
856 //=======================================================================
857 /// Face geometry for the 1D QAdvectionDiffusionReaction elements: Point
858 /// elements
859 //=======================================================================
860 template<unsigned NREAGENT, unsigned NNODE_1D>
862 : public virtual PointElement
863 {
864 public:
865 /// Constructor: Call the constructor for the
866 /// appropriate lower-dimensional QElement
868 };
869
870
871 //==========================================================
872 /// AdvectionDiffusionReaction upgraded to become projectable
873 //==========================================================
874 template<class ADR_ELEMENT>
876 : public virtual ProjectableElement<ADR_ELEMENT>
877 {
878 public:
879 /// Specify the values associated with field fld.
880 /// The information is returned in a vector of pairs which comprise
881 /// the Data object and the value within it, that correspond to field fld.
883 {
884 // Create the vector
885 unsigned nnod = this->nnode();
886 Vector<std::pair<Data*, unsigned>> data_values(nnod);
887
888 // Loop over all nodes
889 for (unsigned j = 0; j < nnod; j++)
890 {
891 // Add the data value associated field: The node itself
892 data_values[j] = std::make_pair(this->node_pt(j), fld);
893 }
894
895 // Return the vector
896 return data_values;
897 }
898
899 /// Number of fields to be projected: Just one
901 {
902 return this->N_reagent;
903 }
904
905 /// Number of history values to be stored for fld-th field
906 /// (includes current value!)
907 unsigned nhistory_values_for_projection(const unsigned& fld)
908 {
909 return this->node_pt(0)->ntstorage();
910 }
911
912 /// Number of positional history values
913 /// (Note: count includes current value!)
915 {
916 return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
917 }
918
919 /// Return Jacobian of mapping and shape functions of field fld
920 /// at local coordinate s
921 double jacobian_and_shape_of_field(const unsigned& fld,
922 const Vector<double>& s,
923 Shape& psi)
924 {
925 unsigned n_dim = this->dim();
926 unsigned n_node = this->nnode();
927 Shape test(n_node);
928 DShape dpsidx(n_node, n_dim), dtestdx(n_node, n_dim);
929 double J = this->dshape_and_dtest_eulerian_adv_diff_react(
930 s, psi, dpsidx, test, dtestdx);
931 return J;
932 }
933
934
935 /// Return interpolated field fld at local coordinate s, at time
936 /// level t (t=0: present; t>0: history values)
937 double get_field(const unsigned& t,
938 const unsigned& fld,
939 const Vector<double>& s)
940 {
941 // Find the index at which the variable is stored
942 unsigned c_nodal_index = this->c_index_adv_diff_react(fld);
943
944 // Local shape function
945 unsigned n_node = this->nnode();
946 Shape psi(n_node);
947
948 // Find values of shape function
949 this->shape(s, psi);
950
951 // Initialise value of c
952 double interpolated_c = 0.0;
953
954 // Sum over the local nodes
955 for (unsigned l = 0; l < n_node; l++)
956 {
957 interpolated_c += this->nodal_value(t, l, c_nodal_index) * psi[l];
958 }
959 return interpolated_c;
960 }
961
962
963 /// Return number of values in field fld: One per node
964 unsigned nvalue_of_field(const unsigned& fld)
965 {
966 return this->nnode();
967 }
968
969
970 /// Return local equation number of value j in field fld.
971 int local_equation(const unsigned& fld, const unsigned& j)
972 {
973 const unsigned c_nodal_index = this->c_index_adv_diff_react(fld);
974 return this->nodal_local_eqn(j, c_nodal_index);
975 }
976 };
977
978
979 //=======================================================================
980 /// Face geometry for element is the same as that for the underlying
981 /// wrapped element
982 //=======================================================================
983 template<class ELEMENT>
985 : public virtual FaceGeometry<ELEMENT>
986 {
987 public:
988 FaceGeometry() : FaceGeometry<ELEMENT>() {}
989 };
990
991
992 //=======================================================================
993 /// Face geometry of the Face Geometry for element is the same as
994 /// that for the underlying wrapped element
995 //=======================================================================
996 template<class ELEMENT>
999 : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
1000 {
1001 public:
1003 };
1004
1005
1006} // namespace oomph
1007
1008#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 Reaction equations using isoparametric el...
Vector< double > *& diff_pt()
Pointer to vector of diffusion coefficients.
AdvectionDiffusionReactionEquations()
Constructor: Initialise the Source_fct_pt, Wind_fct_pt, Reaction_fct_pt to null and initialise the di...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
virtual void get_reaction_deriv_adv_diff_react(const unsigned &ipt, const Vector< double > &C, DenseMatrix< double > &dRdC) const
Get the derivatives of the reaction terms with respect to the concentration variables....
AdvectionDiffusionReactionWindFctPt wind_fct_pt() const
Access function: Pointer to reaction function. Const version.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
AdvectionDiffusionReactionSourceFctPt Source_fct_pt
Pointer to source function:
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
AdvectionDiffusionReactionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
virtual double dshape_and_dtest_eulerian_adv_diff_react(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...
AdvectionDiffusionReactionReactionDerivFctPt Reaction_deriv_fct_pt
Pointer to reaction derivatives.
virtual void get_source_adv_diff_react(const unsigned &ipt, const Vector< double > &x, Vector< double > &source) const
Get source term at (Eulerian) position x. This function is virtual to allow overloading in multi-phys...
const Vector< double > & diff() const
Vector of diffusion coefficients.
virtual void get_wind_adv_diff_react(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(* AdvectionDiffusionReactionWindFctPt)(const double &time, const Vector< double > &x, Vector< double > &wind)
Function pointer to wind function fct(x,w(x)) – x is a Vector!
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and the element Jacobian matrix (wrapper)
void(* AdvectionDiffusionReactionSourceFctPt)(const Vector< double > &x, Vector< double > &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
void dc_dt_adv_diff_react(const unsigned &n, Vector< double > &dc_dt) const
dc/dt at local node n. Uses suitably interpolated value for hanging nodes.
AdvectionDiffusionReactionReactionFctPt reaction_fct_pt() const
Access function: Pointer to reaction function. Const version.
AdvectionDiffusionReactionReactionFctPt & reaction_fct_pt()
Access function: Pointer to reaction function.
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.
Vector< double > * Tau_pt
Pointer to global timescales.
virtual void fill_in_generic_residual_contribution_adv_diff_react(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.
AdvectionDiffusionReactionReactionDerivFctPt & reaction_deriv_fct_pt()
Access function: Pointer to reaction derivatives function.
AdvectionDiffusionReactionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
void output(std::ostream &outfile)
Output with default number of plot points.
void(* AdvectionDiffusionReactionReactionFctPt)(const Vector< double > &c, Vector< double > &R)
Function pointer to reaction terms.
void integrate_reagents(Vector< double > &C) const
Return the integrated reagent concentrations.
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative....
double interpolated_c_adv_diff_react(const Vector< double > &s, const unsigned &i) const
Return FE representation of function value c_i(s) at local coordinate s.
AdvectionDiffusionReactionWindFctPt Wind_fct_pt
Pointer to wind function:
AdvectionDiffusionReactionReactionDerivFctPt reaction_deriv_fct_pt() const
Access function: Pointer to reaction function. Const version.
Vector< double > *& tau_pt()
Pointer to vector of dimensionless timescales.
void operator=(const AdvectionDiffusionReactionEquations &)=delete
Broken assignment operator.
static Vector< double > Default_dimensionless_number
Static default value for the dimensionless numbers.
void output(FILE *file_pt)
C_style output with default number of plot points.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
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...
virtual double dshape_and_dtest_eulerian_at_knot_adv_diff_react(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 compute_norm(double &norm)
Compute norm of the solution: sum of squares of the L2 norm for each reagent.
virtual unsigned c_index_adv_diff_react(const unsigned &i) const
Return the index at which the unknown i-th reagent is stored. The default value, i,...
AdvectionDiffusionReactionEquations(const AdvectionDiffusionReactionEquations &dummy)=delete
Broken copy constructor.
AdvectionDiffusionReactionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
double dc_dt_adv_diff_react(const unsigned &n, const unsigned &r) const
dc_r/dt at local node n. Uses suitably interpolated value for hanging nodes.
AdvectionDiffusionReactionReactionFctPt Reaction_fct_pt
Pointer to reaction 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.
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(* AdvectionDiffusionReactionReactionDerivFctPt)(const Vector< double > &c, DenseMatrix< double > &dRdC)
Function pointer to derivative of reaction terms.
Vector< double > * Diff_pt
Pointer to global diffusion coefficients.
const Vector< double > & tau() const
Vector of dimensionless timescales.
virtual void get_reaction_adv_diff_react(const unsigned &ipt, const Vector< double > &C, Vector< double > &R) const
Get reaction as a function of the given reagent concentrations This function is virtual to allow over...
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
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:879
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...
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1198
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
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: elements.h:3439
AdvectionDiffusionReaction upgraded to become projectable.
unsigned nfields_for_projection()
Number of fields to be projected: Just one.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field (includes current value!)
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
Definition: projection.h:183
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
double dshape_and_dtest_eulerian_at_knot_adv_diff_react(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, 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(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
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.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
unsigned required_nvalue(const unsigned &n) const
Required # of ‘values’ (pinned or dofs) at node n.
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.
double dshape_and_dtest_eulerian_adv_diff_react(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 operator=(const QAdvectionDiffusionReactionElement< NREAGENT, DIM, NNODE_1D > &)=delete
Broken assignment operator.
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 ...
QAdvectionDiffusionReactionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
QAdvectionDiffusionReactionElement(const QAdvectionDiffusionReactionElement< NREAGENT, DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
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
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
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
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
//////////////////////////////////////////////////////////////////// ////////////////////////////////...