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_ADV_DIFF_ELEMENTS_HEADER
28#define OOMPH_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 using isoparametric elements.
46 /// \f[ \frac{\partial^2 u}{\partial x_i^2} = Pe w_i(x_k) \frac{\partial u}{\partial x_i} + 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 double& f);
58
59
60 /// Function pointer to wind function fct(x,w(x)) --
61 /// x is a Vector!
63 Vector<double>& wind);
64
65
66 /// Constructor: Initialise the Source_fct_pt and Wind_fct_pt
67 /// to null and set (pointer to) Peclet number to default
70 {
71 // Set Peclet number to default
73 // Set Peclet Strouhal number to default
75 }
76
77 /// Broken copy constructor
79 delete;
80
81 /// Broken assignment operator
83
84 /// Return the index at which the unknown value
85 /// is stored. The default value, 0, is appropriate for single-physics
86 /// problems, when there is only one variable, the value that satisfies
87 /// the advection-diffusion equation.
88 /// In derived multi-physics elements, this function should be overloaded
89 /// to reflect the chosen storage scheme. Note that these equations require
90 /// that the unknown is always stored at the same index at each node.
91 virtual inline unsigned u_index_adv_diff() const
92 {
93 return 0;
94 }
95
96 /// du/dt at local node n.
97 /// Uses suitably interpolated value for hanging nodes.
98 double du_dt_adv_diff(const unsigned& n) const
99 {
100 // Get the data's timestepper
102
103 // Initialise dudt
104 double dudt = 0.0;
105 // Loop over the timesteps, if there is a non Steady timestepper
107 {
108 // Find the index at which the variable is stored
109 const unsigned u_nodal_index = u_index_adv_diff();
110
111 // Number of timsteps (past & present)
112 const unsigned n_time = time_stepper_pt->ntstorage();
113
114 for (unsigned t = 0; t < n_time; t++)
115 {
116 dudt +=
117 time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
118 }
119 }
120 return dudt;
121 }
122
123 /// Disable ALE, i.e. assert the mesh is not moving -- you do this
124 /// at your own risk!
126 {
127 ALE_is_disabled = true;
128 }
129
130
131 /// (Re-)enable ALE, i.e. take possible mesh motion into account
132 /// when evaluating the time-derivative. Note: By default, ALE is
133 /// enabled, at the expense of possibly creating unnecessary work
134 /// in problems where the mesh is, in fact, stationary.
136 {
137 ALE_is_disabled = false;
138 }
139
140 /// Number of scalars/fields output by this element. Reimplements
141 /// broken virtual function in base class.
142 unsigned nscalar_paraview() const
143 {
144 return DIM + 1;
145 }
146
147 /// Write values of the i-th scalar field at the plot points. Needs
148 /// to be implemented for each new specific element type.
149 void scalar_value_paraview(std::ofstream& file_out,
150 const unsigned& i,
151 const unsigned& nplot) const
152 {
153 // Vector of local coordinates
154 Vector<double> s(DIM);
155
156 // Loop over plot points
157 unsigned num_plot_points = nplot_points_paraview(nplot);
158 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
159 {
160 // Get local coordinates of plot point
161 get_s_plot(iplot, nplot, s);
162
163 // Get Eulerian coordinate of plot point
164 Vector<double> x(DIM);
165 interpolated_x(s, x);
166
167 if (i < DIM)
168 {
169 // Get the wind
170 Vector<double> wind(DIM);
171
172 // Dummy ipt argument needed... ?
173 unsigned ipt = 0;
174 get_wind_adv_diff(ipt, s, x, wind);
175
176 file_out << wind[i] << std::endl;
177 }
178 // Advection Diffusion
179 else if (i == DIM)
180 {
181 file_out << interpolated_u_adv_diff(s) << std::endl;
182 }
183 // Never get here
184 else
185 {
186 std::stringstream error_stream;
187 error_stream << "Advection Diffusion Elements only store " << DIM + 1
188 << " fields " << std::endl;
189 throw OomphLibError(error_stream.str(),
190 OOMPH_CURRENT_FUNCTION,
191 OOMPH_EXCEPTION_LOCATION);
192 }
193 }
194 }
195
196 /// Name of the i-th scalar field. Default implementation
197 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
198 /// overloaded with more meaningful names in specific elements.
199 std::string scalar_name_paraview(const unsigned& i) const
200 {
201 // Winds
202 if (i < DIM)
203 {
204 return "Wind " + StringConversion::to_string(i);
205 }
206 // Advection Diffusion field
207 else if (i == DIM)
208 {
209 return "Advection Diffusion";
210 }
211 // Never get here
212 else
213 {
214 std::stringstream error_stream;
215 error_stream << "Advection Diffusion Elements only store " << DIM + 1
216 << " fields" << std::endl;
217 throw OomphLibError(
218 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
219 // Dummy return
220 return " ";
221 }
222 }
223
224 /// Output with default number of plot points
225 void output(std::ostream& outfile)
226 {
227 unsigned nplot = 5;
228 output(outfile, nplot);
229 }
230
231 /// Output FE representation of soln: x,y,u or x,y,z,u at
232 /// nplot^DIM plot points
233 void output(std::ostream& outfile, const unsigned& nplot);
234
235
236 /// C_style output with default number of plot points
237 void output(FILE* file_pt)
238 {
239 unsigned n_plot = 5;
240 output(file_pt, n_plot);
241 }
242
243 /// C-style output FE representation of soln: x,y,u or x,y,z,u at
244 /// n_plot^DIM plot points
245 void output(FILE* file_pt, const unsigned& n_plot);
246
247
248 /// Output exact soln: x,y,u_exact or x,y,z,u_exact at nplot^DIM plot points
249 void output_fct(std::ostream& outfile,
250 const unsigned& nplot,
252
253 /// Output exact soln: x,y,u_exact or x,y,z,u_exact at
254 /// nplot^DIM plot points (dummy time-dependent version to
255 /// keep intel compiler happy)
256 virtual void output_fct(
257 std::ostream& outfile,
258 const unsigned& nplot,
259 const double& time,
261 {
262 throw OomphLibError("There is no time-dependent output_fct() for "
263 "Advection Diffusion elements",
264 OOMPH_CURRENT_FUNCTION,
265 OOMPH_EXCEPTION_LOCATION);
266 }
267
268
269 /// Get error against and norm of exact solution
270 void compute_error(std::ostream& outfile,
272 double& error,
273 double& norm);
274
275
276 /// Dummy, time dependent error checker
277 void compute_error(std::ostream& outfile,
279 const double& time,
280 double& error,
281 double& norm)
282 {
283 throw OomphLibError(
284 "No time-dependent compute_error() for Advection Diffusion elements",
285 OOMPH_CURRENT_FUNCTION,
286 OOMPH_EXCEPTION_LOCATION);
287 }
288
289
290 /// Access function: Pointer to source function
292 {
293 return Source_fct_pt;
294 }
295
296
297 /// Access function: Pointer to source function. Const version
299 {
300 return Source_fct_pt;
301 }
302
303
304 /// Access function: Pointer to wind function
306 {
307 return Wind_fct_pt;
308 }
309
310
311 /// Access function: Pointer to wind function. Const version
313 {
314 return Wind_fct_pt;
315 }
316
317 /// Peclet number
318 const double& pe() const
319 {
320 return *Pe_pt;
321 }
322
323 /// Pointer to Peclet number
324 double*& pe_pt()
325 {
326 return Pe_pt;
327 }
328
329 /// Peclet number multiplied by Strouhal number
330 const double& pe_st() const
331 {
332 return *PeSt_pt;
333 }
334
335 /// Pointer to Peclet number multipled by Strouha number
336 double*& pe_st_pt()
337 {
338 return PeSt_pt;
339 }
340
341 /// Get source term at (Eulerian) position x. This function is
342 /// virtual to allow overloading in multi-physics problems where
343 /// the strength of the source function might be determined by
344 /// another system of equations
345 inline virtual void get_source_adv_diff(const unsigned& ipt,
346 const Vector<double>& x,
347 double& source) const
348 {
349 // If no source function has been set, return zero
350 if (Source_fct_pt == 0)
351 {
352 source = 0.0;
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(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 wind
382 (*Wind_fct_pt)(x, wind);
383 }
384 }
385
386 /// Get flux: \f$\mbox{flux}[i] = \mbox{d}u / \mbox{d}x_i \f$
387 void get_flux(const Vector<double>& s, Vector<double>& flux) const
388 {
389 // Find out how many nodes there are in the element
390 unsigned n_node = nnode();
391
392 // Get the nodal index at which the unknown is stored
393 unsigned u_nodal_index = u_index_adv_diff();
394
395 // Set up memory for the shape and test functions
396 Shape psi(n_node);
397 DShape dpsidx(n_node, DIM);
398
399 // Call the derivatives of the shape and test functions
400 dshape_eulerian(s, psi, dpsidx);
401
402 // Initialise to zero
403 for (unsigned j = 0; j < DIM; j++)
404 {
405 flux[j] = 0.0;
406 }
407
408 // Loop over nodes
409 for (unsigned l = 0; l < n_node; l++)
410 {
411 // Loop over derivative directions
412 for (unsigned j = 0; j < DIM; j++)
413 {
414 flux[j] += nodal_value(l, u_nodal_index) * dpsidx(l, j);
415 }
416 }
417 }
418
419
420 /// Add the element's contribution to its residual vector (wrapper)
422 {
423 // Call the generic residuals function with flag set to 0 and using
424 // a dummy matrix
426 residuals,
429 0);
430 }
431
432
433 /// Add the element's contribution to its residual vector and
434 /// the element Jacobian matrix (wrapper)
436 DenseMatrix<double>& jacobian)
437 {
438 // Call the generic routine with the flag set to 1
440 residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
441 }
442
443
444 /// Add the element's contribution to its residuals vector,
445 /// jacobian matrix and mass matrix
447 Vector<double>& residuals,
448 DenseMatrix<double>& jacobian,
449 DenseMatrix<double>& mass_matrix)
450 {
451 // Call the generic routine with the flag set to 2
453 residuals, jacobian, mass_matrix, 2);
454 }
455
456
457 /// Return FE representation of function value u(s) at local coordinate s
458 inline double interpolated_u_adv_diff(const Vector<double>& s) const
459 {
460 // Find number of nodes
461 unsigned n_node = nnode();
462
463 // Get the nodal index at which the unknown is stored
464 unsigned u_nodal_index = u_index_adv_diff();
465
466 // Local shape function
467 Shape psi(n_node);
468
469 // Find values of shape function
470 shape(s, psi);
471
472 // Initialise value of u
473 double interpolated_u = 0.0;
474
475 // Loop over the local nodes and sum
476 for (unsigned l = 0; l < n_node; l++)
477 {
478 interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
479 }
480
481 return (interpolated_u);
482 }
483
484
485 /// Return derivative of u at point s with respect to all data
486 /// that can affect its value.
487 /// In addition, return the global equation numbers corresponding to the
488 /// data. This is virtual so that it can be overloaded in the
489 /// refineable version
491 const Vector<double>& s,
492 Vector<double>& du_ddata,
493 Vector<unsigned>& global_eqn_number)
494 {
495 // Find number of nodes
496 const unsigned n_node = nnode();
497
498 // Get the nodal index at which the unknown is stored
499 const unsigned u_nodal_index = u_index_adv_diff();
500
501 // Local shape function
502 Shape psi(n_node);
503
504 // Find values of shape function
505 shape(s, psi);
506
507 // Find the number of dofs associated with interpolated u
508 unsigned n_u_dof = 0;
509 for (unsigned l = 0; l < n_node; l++)
510 {
511 int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
512 // If it's positive add to the count
513 if (global_eqn >= 0)
514 {
515 ++n_u_dof;
516 }
517 }
518
519 // Now resize the storage schemes
520 du_ddata.resize(n_u_dof, 0.0);
521 global_eqn_number.resize(n_u_dof, 0);
522
523 // Loop over the nodes again and set the derivatives
524 unsigned count = 0;
525 for (unsigned l = 0; l < n_node; l++)
526 {
527 // Get the global equation number
528 int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
529 // If it's positive
530 if (global_eqn >= 0)
531 {
532 // Set the global equation number
533 global_eqn_number[count] = global_eqn;
534 // Set the derivative with respect to the unknown
535 du_ddata[count] = psi[l];
536 // Increase the counter
537 ++count;
538 }
539 }
540 }
541
542
543 /// Self-test: Return 0 for OK
544 unsigned self_test();
545
546 protected:
547 /// Shape/test functions and derivs w.r.t. to global coords at
548 /// local coord. s; return Jacobian of mapping
550 const Vector<double>& s,
551 Shape& psi,
552 DShape& dpsidx,
553 Shape& test,
554 DShape& dtestdx) const = 0;
555
556 /// Shape/test functions and derivs w.r.t. to global coords at
557 /// integration point ipt; return Jacobian of mapping
559 const unsigned& ipt,
560 Shape& psi,
561 DShape& dpsidx,
562 Shape& test,
563 DShape& dtestdx) const = 0;
564
565 /// Add the element's contribution to its residual vector only
566 /// (if flag=and/or element Jacobian matrix
568 Vector<double>& residuals,
569 DenseMatrix<double>& jacobian,
570 DenseMatrix<double>& mass_matrix,
571 unsigned flag);
572
573 /// Pointer to global Peclet number
574 double* Pe_pt;
575
576 /// Pointer to global Peclet number multiplied by Strouhal number
577 double* PeSt_pt;
578
579 /// Pointer to source function:
581
582 /// Pointer to wind function:
584
585 /// Boolean flag to indicate if ALE formulation is disabled when
586 /// time-derivatives are computed. Only set to true if you're sure
587 /// that the mesh is stationary.
589
590 private:
591 /// Static default value for the Peclet number
593 };
594
595
596 /// ////////////////////////////////////////////////////////////////////////
597 /// ////////////////////////////////////////////////////////////////////////
598 /// ////////////////////////////////////////////////////////////////////////
599
600
601 //======================================================================
602 /// QAdvectionDiffusionElement elements are
603 /// linear/quadrilateral/brick-shaped Advection Diffusion elements with
604 /// isoparametric interpolation for the function.
605 //======================================================================
606 template<unsigned DIM, unsigned NNODE_1D>
608 : public virtual QElement<DIM, NNODE_1D>,
609 public virtual AdvectionDiffusionEquations<DIM>
610 {
611 private:
612 /// Static array of ints to hold number of variables at
613 /// nodes: Initial_Nvalue[n]
614 static const unsigned Initial_Nvalue;
615
616 public:
617 /// Constructor: Call constructors for QElement and
618 /// Advection Diffusion equations
620 : QElement<DIM, NNODE_1D>(), AdvectionDiffusionEquations<DIM>()
621 {
622 }
623
624 /// Broken copy constructor
626 const QAdvectionDiffusionElement<DIM, NNODE_1D>& dummy) = delete;
627
628 /// Broken assignment operator
630
631 /// Required # of `values' (pinned or dofs)
632 /// at node n
633 inline unsigned required_nvalue(const unsigned& n) const
634 {
635 return Initial_Nvalue;
636 }
637
638 /// Output function:
639 /// x,y,u or x,y,z,u
640 void output(std::ostream& outfile)
641 {
643 }
644
645 /// Output function:
646 /// x,y,u or x,y,z,u at n_plot^DIM plot points
647 void output(std::ostream& outfile, const unsigned& n_plot)
648 {
650 }
651
652
653 /// C-style output function:
654 /// x,y,u or x,y,z,u
655 void output(FILE* file_pt)
656 {
658 }
659
660 /// C-style output function:
661 /// x,y,u or x,y,z,u at n_plot^DIM plot points
662 void output(FILE* file_pt, const unsigned& n_plot)
663 {
665 }
666
667 /// Output function for an exact solution:
668 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
669 void output_fct(std::ostream& outfile,
670 const unsigned& n_plot,
672 {
674 outfile, n_plot, exact_soln_pt);
675 }
676
677
678 /// Output function for a time-dependent exact solution.
679 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
680 /// (Calls the steady version)
681 void output_fct(std::ostream& outfile,
682 const unsigned& n_plot,
683 const double& time,
685 {
687 outfile, n_plot, time, exact_soln_pt);
688 }
689
690
691 protected:
692 /// Shape, test functions & derivs. w.r.t. to global coords. Return
693 /// Jacobian.
695 Shape& psi,
696 DShape& dpsidx,
697 Shape& test,
698 DShape& dtestdx) const;
699
700 /// Shape, test functions & derivs. w.r.t. to global coords. at
701 /// integration point ipt. Return Jacobian.
703 const unsigned& ipt,
704 Shape& psi,
705 DShape& dpsidx,
706 Shape& test,
707 DShape& dtestdx) const;
708 };
709
710 // Inline functions:
711
712
713 //======================================================================
714 /// Define the shape functions and test functions and derivatives
715 /// w.r.t. global coordinates and return Jacobian of mapping.
716 ///
717 /// Galerkin: Test functions = shape functions
718 //======================================================================
719 template<unsigned DIM, unsigned NNODE_1D>
722 Shape& psi,
723 DShape& dpsidx,
724 Shape& test,
725 DShape& dtestdx) const
726 {
727 // Call the geometrical shape functions and derivatives
728 double J = this->dshape_eulerian(s, psi, dpsidx);
729
730 // Loop over the test functions and derivatives and set them equal to the
731 // shape functions
732 for (unsigned i = 0; i < NNODE_1D; i++)
733 {
734 test[i] = psi[i];
735 for (unsigned j = 0; j < DIM; j++)
736 {
737 dtestdx(i, j) = dpsidx(i, j);
738 }
739 }
740
741 // Return the jacobian
742 return J;
743 }
744
745
746 //======================================================================
747 /// Define the shape functions and test functions and derivatives
748 /// w.r.t. global coordinates and return Jacobian of mapping.
749 ///
750 /// Galerkin: Test functions = shape functions
751 //======================================================================
752 template<unsigned DIM, unsigned NNODE_1D>
755 Shape& psi,
756 DShape& dpsidx,
757 Shape& test,
758 DShape& dtestdx) const
759 {
760 // Call the geometrical shape functions and derivatives
761 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
762
763 // Set the test functions equal to the shape functions (pointer copy)
764 test = psi;
765 dtestdx = dpsidx;
766
767 // Return the jacobian
768 return J;
769 }
770
771
772 /// /////////////////////////////////////////////////////////////////////
773 /// /////////////////////////////////////////////////////////////////////
774 /// /////////////////////////////////////////////////////////////////////
775
776
777 //=======================================================================
778 /// Face geometry for the QAdvectionDiffusionElement elements:
779 /// The spatial dimension of the face elements is one lower than that
780 /// of the bulk element but they have the same number of points along
781 /// their 1D edges.
782 //=======================================================================
783 template<unsigned DIM, unsigned NNODE_1D>
785 : public virtual QElement<DIM - 1, NNODE_1D>
786 {
787 public:
788 /// Constructor: Call the constructor for the
789 /// appropriate lower-dimensional QElement
790 FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
791 };
792
793
794 /// /////////////////////////////////////////////////////////////////////
795 /// /////////////////////////////////////////////////////////////////////
796 /// /////////////////////////////////////////////////////////////////////
797
798
799 //=======================================================================
800 /// Face geometry for the 1D QAdvectionDiffusion elements: Point elements
801 //=======================================================================
802 template<unsigned NNODE_1D>
804 : public virtual PointElement
805 {
806 public:
807 /// Constructor: Call the constructor for the
808 /// appropriate lower-dimensional QElement
810 };
811
812} // namespace oomph
813
814#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A class for all elements that solve the Advection Diffusion equations using isoparametric elements.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
void output(FILE *file_pt)
C_style output with default number of plot points.
const double & pe_st() const
Peclet number multiplied by Strouhal number.
AdvectionDiffusionEquations()
Constructor: Initialise the Source_fct_pt and Wind_fct_pt to null and set (pointer to) Peclet number ...
double interpolated_u_adv_diff(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
double du_dt_adv_diff(const unsigned &n) const
du/dt at local node n. Uses suitably interpolated value for hanging nodes.
void output(std::ostream &outfile)
Output with default number of plot points.
virtual void get_wind_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(* AdvectionDiffusionWindFctPt)(const Vector< double > &x, Vector< double > &wind)
Function pointer to wind function fct(x,w(x)) – x is a Vector!
virtual double dshape_and_dtest_eulerian_at_knot_adv_diff(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one,...
virtual void fill_in_generic_residual_contribution_adv_diff(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add the element's contribution to its residual vector only (if flag=and/or element Jacobian matrix.
double * Pe_pt
Pointer to global Peclet number.
virtual void get_source_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...
double *& pe_pt()
Pointer to Peclet number.
AdvectionDiffusionEquations(const AdvectionDiffusionEquations &dummy)=delete
Broken copy constructor.
AdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and the element Jacobian matrix (wrapper)
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class.
AdvectionDiffusionSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
static double Default_peclet_number
Static default value for the Peclet number.
virtual void dinterpolated_u_adv_diff_ddata(const Vector< double > &s, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
Return derivative of u at point s with respect to all data that can affect its value....
AdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
void operator=(const AdvectionDiffusionEquations &)=delete
Broken assignment operator.
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative....
const double & pe() const
Peclet number.
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
AdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
void(* AdvectionDiffusionSourceFctPt)(const Vector< double > &x, double &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
virtual unsigned u_index_adv_diff() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
AdvectionDiffusionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
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.
virtual double dshape_and_dtest_eulerian_adv_diff(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
unsigned self_test()
Self-test: Return 0 for OK.
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_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
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.
AdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
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...
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2862
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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 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
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
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 ...
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.
QAdvectionDiffusionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
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.
unsigned required_nvalue(const unsigned &n) const
Required # of ‘values’ (pinned or dofs) at node n.
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
QAdvectionDiffusionElement(const QAdvectionDiffusionElement< DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
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.
static const unsigned Initial_Nvalue
Static array of ints to hold number of variables at nodes: Initial_Nvalue[n].
double dshape_and_dtest_eulerian_adv_diff(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
void operator=(const QAdvectionDiffusionElement< DIM, NNODE_1D > &)=delete
Broken assignment operator.
double dshape_and_dtest_eulerian_at_knot_adv_diff(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:389
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types).
//////////////////////////////////////////////////////////////////// ////////////////////////////////...