scalar_advection_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 scalar advection elements
27
28#ifndef OOMPH_SCALAR_ADVECTION_ELEMENTS_HEADER
29#define OOMPH_SCALAR_ADVECTION_ELEMENTS_HEADER
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
37#include "../generic/Qelements.h"
38#include "../generic/Qspectral_elements.h"
39#include "../generic/dg_elements.h"
40
41namespace oomph
42{
43 //==============================================================
44 /// Base class for advection equations
45 //=============================================================
46 template<unsigned DIM>
48 {
49 /// Typedef for a wind function as a possible function of position
50 typedef void (*ScalarAdvectionWindFctPt)(const Vector<double>& x,
51 Vector<double>& wind);
52
53 /// Function pointer to the wind function
55
56 protected:
57 /// A single flux is interpolated
58 inline unsigned nflux() const
59 {
60 return 1;
61 }
62
63 /// Return the flux as a function of the unknown
64 void flux(const Vector<double>& u, DenseMatrix<double>& f);
65
66 /// Return the flux derivatives as a function of the unknowns
67 void dflux_du(const Vector<double>& u, RankThreeTensor<double>& df_du);
68
69 /// Return the wind at a given position
70 inline virtual void get_wind_scalar_adv(const unsigned& ipt,
71 const Vector<double>& s,
72 const Vector<double>& x,
73 Vector<double>& wind) const
74 {
75 // If no wind function has been set, return zero
76 if (Wind_fct_pt == 0)
77 {
78 for (unsigned i = 0; i < DIM; i++)
79 {
80 wind[i] = 0.0;
81 }
82 }
83 else
84 {
85 // Get wind
86 (*Wind_fct_pt)(x, wind);
87 }
88 }
89
90 public:
91 /// Constructor
93 {
94 }
95
96 /// Access function: Pointer to wind function
98 {
99 return Wind_fct_pt;
100 }
101
102 /// Access function: Pointer to wind function. Const version
104 {
105 return Wind_fct_pt;
106 }
107
108 /// The number of unknowns at each node is the number of values
109 unsigned required_nvalue(const unsigned& n) const
110 {
111 return 1;
112 }
113
114 /// Compute the error and norm of solution integrated over the element
115 /// Does not plot the error in the outfile
117 std::ostream& outfile,
119 const double& t,
120 Vector<double>& error,
121 Vector<double>& norm)
122 {
123 // Find the number of fluxes
124 const unsigned n_flux = this->nflux();
125 // Find the number of nodes
126 const unsigned n_node = this->nnode();
127 // Storage for the shape function and derivatives of shape function
128 Shape psi(n_node);
129 DShape dpsidx(n_node, DIM);
130
131 // Find the number of integration points
132 unsigned n_intpt = this->integral_pt()->nweight();
133
134 error.resize(n_flux);
135 norm.resize(n_flux);
136 for (unsigned i = 0; i < n_flux; i++)
137 {
138 error[i] = 0.0;
139 norm[i] = 0.0;
140 }
141
142 Vector<double> s(DIM);
143
144 // Loop over the integration points
145 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
146 {
147 // Get the shape functions at the knot
148 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
149
150 // Get the integral weight
151 double W = this->integral_pt()->weight(ipt) * J;
152
153 // Storage for the local functions
155 Vector<double> interpolated_u(n_flux, 0.0);
156
157 // Loop over the shape functions
158 for (unsigned l = 0; l < n_node; l++)
159 {
160 // Locally cache the shape function
161 const double psi_ = psi(l);
162 for (unsigned i = 0; i < DIM; i++)
163 {
164 interpolated_x[i] += this->nodal_position(l, i) * psi_;
165 }
166
167 for (unsigned i = 0; i < n_flux; i++)
168 {
169 // Calculate the velocity and tangent vector
170 interpolated_u[i] += this->nodal_value(l, i) * psi_;
171 }
172 }
173
174 // Get the global wind
175 Vector<double> wind(DIM);
176 this->get_wind_scalar_adv(ipt, s, interpolated_x, wind);
177
178 // Rescale the position
179 for (unsigned i = 0; i < DIM; i++)
180 {
181 interpolated_x[i] -= wind[i] * t;
182 }
183
184 // Now get the initial condition at this value of x
185 Vector<double> exact_u(n_flux, 0.0);
186 (*initial_condition_pt)(0.0, interpolated_x, exact_u);
187
188 // Loop over the unknowns
189 for (unsigned i = 0; i < n_flux; i++)
190 {
191 error[i] += pow((interpolated_u[i] - exact_u[i]), 2.0) * W;
192 norm[i] += interpolated_u[i] * interpolated_u[i] * W;
193 }
194 }
195 }
196 };
197
198
199 template<unsigned DIM, unsigned NNODE_1D>
201 : public virtual QSpectralElement<DIM, NNODE_1D>,
202 public virtual ScalarAdvectionEquations<DIM>
203 {
204 public:
205 /// Constructor: Call constructors for QElement and
206 /// Advection Diffusion equations
208 : QSpectralElement<DIM, NNODE_1D>(), ScalarAdvectionEquations<DIM>()
209 {
210 }
211
212 /// Broken copy constructor
215
216 /// Broken assignment operator
217 // Commented out broken assignment operator because this can lead to a
218 // conflict warning when used in the virtual inheritence hierarchy.
219 // Essentially the compiler doesn't realise that two separate
220 // implementations of the broken function are the same and so, quite
221 // rightly, it shouts.
222 /*void operator=(
223 const QSpectralScalarAdvectionElement<DIM,NNODE_1D>&) = delete;*/
224
225 /// Required # of `values' (pinned or dofs)
226 /// at node n
227 inline unsigned required_nvalue(const unsigned& n) const
228 {
229 return 1;
230 }
231
232 /// Output function:
233 /// x,y,u or x,y,z,u
234 void output(std::ostream& outfile)
235 {
237 }
238
239 /// Output function:
240 /// x,y,u or x,y,z,u at n_plot^DIM plot points
241 void output(std::ostream& outfile, const unsigned& n_plot)
242 {
244 }
245
246
247 /*/// C-style output function:
248 /// x,y,u or x,y,z,u
249 void output(FILE* file_pt)
250 {
251 ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt);
252 }
253
254 /// C-style output function:
255 /// x,y,u or x,y,z,u at n_plot^DIM plot points
256 void output(FILE* file_pt, const unsigned &n_plot)
257 {
258 ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt,n_plot);
259 }
260
261 /// Output function for an exact solution:
262 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
263 void output_fct(std::ostream &outfile, const unsigned &n_plot,
264 FiniteElement::SteadyExactSolutionFctPt
265 exact_soln_pt)
266 {
267 ScalarAdvectionEquations<NFLUX,DIM>::
268 output_fct(outfile,n_plot,exact_soln_pt);}
269
270
271 /// Output function for a time-dependent exact solution.
272 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
273 /// (Calls the steady version)
274 void output_fct(std::ostream &outfile, const unsigned &n_plot,
275 const double& time,
276 FiniteElement::UnsteadyExactSolutionFctPt
277 exact_soln_pt)
278 {
279 ScalarAdvectionEquations<NFLUX,DIM>::
280 output_fct(outfile,n_plot,time,exact_soln_pt);
281 }*/
282
283
284 protected:
285 /// Shape, test functions & derivs. w.r.t. to global coords. Return
286 /// Jacobian.
288 const Vector<double>& s,
289 Shape& psi,
290 DShape& dpsidx,
291 Shape& test,
292 DShape& dtestdx) const;
293
294 /// Shape, test functions & derivs. w.r.t. to global coords. at
295 /// integration point ipt. Return Jacobian.
297 const unsigned& ipt,
298 Shape& psi,
299 DShape& dpsidx,
300 Shape& test,
301 DShape& dtestdx) const;
302 };
303
304 // Inline functions:
305
306
307 //======================================================================
308 /// Define the shape functions and test functions and derivatives
309 /// w.r.t. global coordinates and return Jacobian of mapping.
310 ///
311 /// Galerkin: Test functions = shape functions
312 //======================================================================
313 template<unsigned DIM, unsigned NNODE_1D>
316 Shape& psi,
317 DShape& dpsidx,
318 Shape& test,
319 DShape& dtestdx) const
320 {
321 // Call the geometrical shape functions and derivatives
322 double J = this->dshape_eulerian(s, psi, dpsidx);
323
324 // Loop over the test functions and derivatives and set them equal to the
325 // shape functions
326 for (unsigned i = 0; i < NNODE_1D; i++)
327 {
328 test[i] = psi[i];
329 for (unsigned j = 0; j < DIM; j++)
330 {
331 dtestdx(i, j) = dpsidx(i, j);
332 }
333 }
334
335 // Return the jacobian
336 return J;
337 }
338
339
340 //======================================================================
341 /// Define the shape functions and test functions and derivatives
342 /// w.r.t. global coordinates and return Jacobian of mapping.
343 ///
344 /// Galerkin: Test functions = shape functions
345 //======================================================================
346 template<unsigned DIM, unsigned NNODE_1D>
349 Shape& psi,
350 DShape& dpsidx,
351 Shape& test,
352 DShape& dtestdx) const
353 {
354 // Call the geometrical shape functions and derivatives
355 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
356
357 // Set the test functions equal to the shape functions (pointer copy)
358 test = psi;
359 dtestdx = dpsidx;
360
361 // Return the jacobian
362 return J;
363 }
364
365
366 /// /////////////////////////////////////////////////////////////////////
367 /// /////////////////////////////////////////////////////////////////////
368 /// /////////////////////////////////////////////////////////////////////
369
370
371 //=======================================================================
372 /// Face geometry for the QScalarAdvectionElement elements:
373 /// The spatial dimension of the face elements is one lower than that
374 /// of the bulk element but they have the same number of points along
375 /// their 1D edges.
376 //=======================================================================
377 template<unsigned DIM, unsigned NNODE_1D>
379 : public virtual QSpectralElement<DIM - 1, NNODE_1D>
380 {
381 public:
382 /// Constructor: Call the constructor for the
383 /// appropriate lower-dimensional QElement
384 FaceGeometry() : QSpectralElement<DIM - 1, NNODE_1D>() {}
385 };
386
387
388 //======================================================================
389 /// FaceElement for Discontinuous Galerkin Problems
390 //======================================================================
391 template<class ELEMENT>
392 class DGScalarAdvectionFaceElement : public virtual FaceGeometry<ELEMENT>,
393 public virtual DGFaceElement
394 {
395 public:
396 /// Constructor
398 const int& face_index)
399 : FaceGeometry<ELEMENT>(), DGFaceElement()
400 {
401 // Attach geometric information to the element
402 // N.B. This function also assigns nbulk_value from required_nvalue
403 // of the bulk element.
404 element_pt->build_face_element(face_index, this);
405 }
406
407
408 // There is a single required n_flux
409 unsigned required_nflux()
410 {
411 return 1;
412 }
413
414 /// Calculate the normal flux, so this is the dot product of the
415 /// numerical flux with n_in
417 const Vector<double>& u_int,
418 const Vector<double>& u_ext,
419 Vector<double>& flux)
420 {
421 const unsigned dim = this->nodal_dimension();
422 Vector<double> Wind(dim);
423 Vector<double> s, x;
424 // Dummy integration point
425 unsigned ipt = 0;
426 dynamic_cast<ELEMENT*>(this->bulk_element_pt())
427 ->get_wind_scalar_adv(ipt, s, x, Wind);
428
429 // Now we can work this out for standard upwind advection
430 double dot = 0.0;
431 for (unsigned i = 0; i < dim; i++)
432 {
433 dot += Wind[i] * n_out[i];
434 }
435
436 const unsigned n_value = this->required_nflux();
437 if (dot >= 0.0)
438 {
439 for (unsigned n = 0; n < n_value; n++)
440 {
441 flux[n] = dot * u_int[n];
442 }
443 }
444 else
445 {
446 for (unsigned n = 0; n < n_value; n++)
447 {
448 flux[n] = dot * u_ext[n];
449 }
450 }
451
452 // flux[0] = 0.5*(u_int[0]+u_ext[0])*n_out[0];
453 }
454 };
455
456 //=================================================================
457 /// General DGScalarAdvectionClass. Establish the template parameters
458 //===================================================================
459 template<unsigned DIM, unsigned NNODE_1D>
461 {
462 };
463
464
465 //==================================================================
466 // Specialization for 1D DG Advection element
467 //==================================================================
468 template<unsigned NNODE_1D>
470 : public QSpectralScalarAdvectionElement<1, NNODE_1D>, public DGElement
471 {
472 friend class DGScalarAdvectionFaceElement<
474
476
477 public:
478 // There is a single required n_flux
479 unsigned required_nflux()
480 {
481 return 1;
482 }
483
484 // Calculate averages
485 void calculate_element_averages(double*& average_value)
486 {
488 }
489
490 // Constructor
493 {
494 }
495
497
499 {
500 // Make the two faces
501 Face_element_pt.resize(2);
502 // Make the face on the left
503 Face_element_pt[0] = new DGScalarAdvectionFaceElement<
505 // Make the face on the right
506 Face_element_pt[1] = new DGScalarAdvectionFaceElement<
508 }
509
510
511 /// Compute the residuals for the Navier--Stokes equations;
512 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
514 Vector<double>& residuals,
515 DenseMatrix<double>& jacobian,
516 DenseMatrix<double>& mass_matrix,
517 unsigned flag)
518 {
521 residuals, jacobian, mass_matrix, flag);
522
523 this->add_flux_contributions_to_residuals(residuals, jacobian, flag);
524 }
525
526
527 //============================================================================
528 /// Function that returns the current value of the residuals
529 /// multiplied by the inverse mass matrix (virtual so that it can be
530 /// overloaded specific elements in which time and memory saving tricks can
531 /// be applied)
532 //============================================================================
534 {
535 // If there are external data this is not going to work
536 if (nexternal_data() > 0)
537 {
538 std::ostringstream error_stream;
539 error_stream
540 << "Cannot use a discontinuous formulation for the mass matrix when\n"
541 << "there are external data.\n "
542 << "Do not call Problem::enable_discontinuous_formulation()\n";
543
544 throw OomphLibError(
545 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
546 }
547
548
549 // Now let's assemble stuff
550 const unsigned n_dof = this->ndof();
551
552 // Resize and initialise the vector that will holds the residuals
553 minv_res.resize(n_dof);
554 for (unsigned n = 0; n < n_dof; n++)
555 {
556 minv_res[n] = 0.0;
557 }
558
559 // If we are recycling the mass matrix
560 if (Mass_matrix_reuse_is_enabled && Mass_matrix_has_been_computed)
561 {
562 // Get the residuals
563 this->fill_in_contribution_to_residuals(minv_res);
564 }
565 // Otherwise
566 else
567 {
568 // Temporary mass matrix
569 DenseDoubleMatrix M(n_dof, n_dof, 0.0);
570
571 // Get the local mass matrix and residuals
572 this->fill_in_contribution_to_mass_matrix(minv_res, M);
573
574 // Store the diagonal entries
575 Inverse_mass_diagonal.clear();
576 for (unsigned n = 0; n < n_dof; n++)
577 {
578 Inverse_mass_diagonal.push_back(1.0 / M(n, n));
579 }
580
581 // The mass matrix has been computed
582 Mass_matrix_has_been_computed = true;
583 }
584
585 for (unsigned n = 0; n < n_dof; n++)
586 {
587 minv_res[n] *= Inverse_mass_diagonal[n];
588 }
589 }
590 };
591
592
593 //=======================================================================
594 /// Face geometry of the 1D DG elements
595 //=======================================================================
596 template<unsigned NNODE_1D>
598 : public virtual PointElement
599 {
600 public:
602 };
603
604
605 //==================================================================
606 /// Specialisation for 2D DG Elements
607 //==================================================================
608 template<unsigned NNODE_1D>
610 : public QSpectralScalarAdvectionElement<2, NNODE_1D>, public DGElement
611 {
612 friend class DGScalarAdvectionFaceElement<
614
615 public:
616 // Calculate averages
617 void calculate_element_averages(double*& average_value)
618 {
620 }
621
622 // There is a single required n_flux
623 unsigned required_nflux()
624 {
625 return 1;
626 }
627
628 // Constructor
631 {
632 }
633
635
637 {
638 Face_element_pt.resize(4);
639 Face_element_pt[0] = new DGScalarAdvectionFaceElement<
641 Face_element_pt[1] = new DGScalarAdvectionFaceElement<
643 Face_element_pt[2] = new DGScalarAdvectionFaceElement<
645 Face_element_pt[3] = new DGScalarAdvectionFaceElement<
647 }
648
649
650 /// Compute the residuals for the Navier--Stokes equations;
651 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
653 Vector<double>& residuals,
654 DenseMatrix<double>& jacobian,
655 DenseMatrix<double>& mass_matrix,
656 unsigned flag)
657 {
660 residuals, jacobian, mass_matrix, flag);
661
662 this->add_flux_contributions_to_residuals(residuals, jacobian, flag);
663 }
664 };
665
666
667 //=======================================================================
668 /// Face geometry of the DG elements
669 //=======================================================================
670 template<unsigned NNODE_1D>
672 : public virtual QSpectralElement<1, NNODE_1D>
673 {
674 public:
675 FaceGeometry() : QSpectralElement<1, NNODE_1D>() {}
676 };
677
678
679 //=============================================================
680 /// Non-spectral version of the classes
681 //============================================================
682 template<unsigned DIM, unsigned NNODE_1D>
683 class QScalarAdvectionElement : public virtual QElement<DIM, NNODE_1D>,
684 public virtual ScalarAdvectionEquations<DIM>
685 {
686 public:
687 /// Constructor: Call constructors for QElement and
688 /// Advection Diffusion equations
690 : QElement<DIM, NNODE_1D>(), ScalarAdvectionEquations<DIM>()
691 {
692 }
693
694 /// Broken copy constructor
696 const QScalarAdvectionElement<DIM, NNODE_1D>& dummy) = delete;
697
698 /// Broken assignment operator
699 /*void operator=(
700 const QScalarAdvectionElement<DIM,NNODE_1D>&) = delete;*/
701
702 /// Required # of `values' (pinned or dofs)
703 /// at node n
704 inline unsigned required_nvalue(const unsigned& n) const
705 {
706 return 1;
707 }
708
709 /// Output function:
710 /// x,y,u or x,y,z,u
711 void output(std::ostream& outfile)
712 {
714 }
715
716 /// Output function:
717 /// x,y,u or x,y,z,u at n_plot^DIM plot points
718 void output(std::ostream& outfile, const unsigned& n_plot)
719 {
721 }
722
723
724 /*/// C-style output function:
725 /// x,y,u or x,y,z,u
726 void output(FILE* file_pt)
727 {
728 ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt);
729 }
730
731 /// C-style output function:
732 /// x,y,u or x,y,z,u at n_plot^DIM plot points
733 void output(FILE* file_pt, const unsigned &n_plot)
734 {
735 ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt,n_plot);
736 }
737
738 /// Output function for an exact solution:
739 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
740 void output_fct(std::ostream &outfile, const unsigned &n_plot,
741 FiniteElement::SteadyExactSolutionFctPt
742 exact_soln_pt)
743 {
744 ScalarAdvectionEquations<NFLUX,DIM>::
745 output_fct(outfile,n_plot,exact_soln_pt);}
746
747
748 /// Output function for a time-dependent exact solution.
749 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
750 /// (Calls the steady version)
751 void output_fct(std::ostream &outfile, const unsigned &n_plot,
752 const double& time,
753 FiniteElement::UnsteadyExactSolutionFctPt
754 exact_soln_pt)
755 {
756 ScalarAdvectionEquations<NFLUX,DIM>::
757 output_fct(outfile,n_plot,time,exact_soln_pt);
758 }*/
759
760
761 protected:
762 /// Shape, test functions & derivs. w.r.t. to global coords. Return
763 /// Jacobian.
765 const Vector<double>& s,
766 Shape& psi,
767 DShape& dpsidx,
768 Shape& test,
769 DShape& dtestdx) const;
770
771 /// Shape, test functions & derivs. w.r.t. to global coords. at
772 /// integration point ipt. Return Jacobian.
774 const unsigned& ipt,
775 Shape& psi,
776 DShape& dpsidx,
777 Shape& test,
778 DShape& dtestdx) const;
779 };
780
781 // Inline functions:
782
783
784 //======================================================================
785 /// Define the shape functions and test functions and derivatives
786 /// w.r.t. global coordinates and return Jacobian of mapping.
787 ///
788 /// Galerkin: Test functions = shape functions
789 //======================================================================
790 template<unsigned DIM, unsigned NNODE_1D>
793 Shape& psi,
794 DShape& dpsidx,
795 Shape& test,
796 DShape& dtestdx) const
797 {
798 // Call the geometrical shape functions and derivatives
799 double J = this->dshape_eulerian(s, psi, dpsidx);
800
801 // Loop over the test functions and derivatives and set them equal to the
802 // shape functions
803 for (unsigned i = 0; i < NNODE_1D; i++)
804 {
805 test[i] = psi[i];
806 for (unsigned j = 0; j < DIM; j++)
807 {
808 dtestdx(i, j) = dpsidx(i, j);
809 }
810 }
811
812 // Return the jacobian
813 return J;
814 }
815
816
817 //======================================================================
818 /// Define the shape functions and test functions and derivatives
819 /// w.r.t. global coordinates and return Jacobian of mapping.
820 ///
821 /// Galerkin: Test functions = shape functions
822 //======================================================================
823 template<unsigned DIM, unsigned NNODE_1D>
826 Shape& psi,
827 DShape& dpsidx,
828 Shape& test,
829 DShape& dtestdx) const
830 {
831 // Call the geometrical shape functions and derivatives
832 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
833
834 // Set the test functions equal to the shape functions (pointer copy)
835 test = psi;
836 dtestdx = dpsidx;
837
838 // Return the jacobian
839 return J;
840 }
841
842
843 /// /////////////////////////////////////////////////////////////////////
844 /// /////////////////////////////////////////////////////////////////////
845 /// /////////////////////////////////////////////////////////////////////
846
847
848 //=======================================================================
849 /// Face geometry for the QScalarAdvectionElement elements:
850 /// The spatial dimension of the face elements is one lower than that
851 /// of the bulk element but they have the same number of points along
852 /// their 1D edges.
853 //=======================================================================
854 template<unsigned DIM, unsigned NNODE_1D>
856 : public virtual QElement<DIM - 1, NNODE_1D>
857 {
858 public:
859 /// Constructor: Call the constructor for the
860 /// appropriate lower-dimensional QElement
861 FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
862 };
863
864
865 //=================================================================
866 /// General DGScalarAdvectionClass. Establish the template parameters
867 //===================================================================
868 template<unsigned DIM, unsigned NNODE_1D>
870 {
871 };
872
873
874 //==================================================================
875 // Specialization for 1D DG Advection element
876 //==================================================================
877 template<unsigned NNODE_1D>
878 class DGScalarAdvectionElement<1, NNODE_1D>
879 : public QScalarAdvectionElement<1, NNODE_1D>, public DGElement
880 {
881 friend class DGScalarAdvectionFaceElement<
882 DGScalarAdvectionElement<1, NNODE_1D>>;
883
884 public:
885 // There is a single required n_flux
886 unsigned required_nflux()
887 {
888 return 1;
889 }
890
891 // Calculate averages
892 void calculate_element_averages(double*& average_value)
893 {
895 }
896
897 // Constructor
899 : QScalarAdvectionElement<1, NNODE_1D>(), DGElement()
900 {
901 }
902
904
906 {
907 // Make the two faces
908 Face_element_pt.resize(2);
909 // Make the face on the left
910 Face_element_pt[0] =
912 this, -1);
913 // Make the face on the right
914 Face_element_pt[1] =
916 this, +1);
917 }
918
919
920 /// Compute the residuals for the Navier--Stokes equations;
921 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
923 Vector<double>& residuals,
924 DenseMatrix<double>& jacobian,
925 DenseMatrix<double>& mass_matrix,
926 unsigned flag)
927 {
930 residuals, jacobian, mass_matrix, flag);
931
932 this->add_flux_contributions_to_residuals(residuals, jacobian, flag);
933 }
934 };
935
936
937 //=======================================================================
938 /// Face geometry of the 1D DG elements
939 //=======================================================================
940 template<unsigned NNODE_1D>
942 : public virtual PointElement
943 {
944 public:
946 };
947
948
949 //==================================================================
950 /// Specialisation for 2D DG Elements
951 //==================================================================
952 template<unsigned NNODE_1D>
953 class DGScalarAdvectionElement<2, NNODE_1D>
954 : public QScalarAdvectionElement<2, NNODE_1D>, public DGElement
955 {
956 friend class DGScalarAdvectionFaceElement<
957 DGScalarAdvectionElement<2, NNODE_1D>>;
958
959 public:
960 // There is a single required n_flux
961 unsigned required_nflux()
962 {
963 return 1;
964 }
965
966 // Calculate averages
967 void calculate_element_averages(double*& average_value)
968 {
970 }
971
972 // Constructor
974 : QScalarAdvectionElement<2, NNODE_1D>(), DGElement()
975 {
976 }
977
979
981 {
982 Face_element_pt.resize(4);
983 Face_element_pt[0] =
985 this, 2);
986 Face_element_pt[1] =
988 this, 1);
989 Face_element_pt[2] =
991 this, -2);
992 Face_element_pt[3] =
994 this, -1);
995 }
996
997
998 /// Compute the residuals for the Navier--Stokes equations;
999 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
1001 Vector<double>& residuals,
1002 DenseMatrix<double>& jacobian,
1003 DenseMatrix<double>& mass_matrix,
1004 unsigned flag)
1005 {
1008 residuals, jacobian, mass_matrix, flag);
1009
1010 this->add_flux_contributions_to_residuals(residuals, jacobian, flag);
1011 }
1012 };
1013
1014
1015 //=======================================================================
1016 /// Face geometry of the DG elements
1017 //=======================================================================
1018 template<unsigned NNODE_1D>
1020 : public virtual QElement<1, NNODE_1D>
1021 {
1022 public:
1023 FaceGeometry() : QElement<1, NNODE_1D>() {}
1024 };
1025
1026
1027} // namespace oomph
1028
1029#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 Base class for DGElements.
Definition: dg_elements.h:161
Base class for Discontinuous Galerkin Faces. These are responsible for calculating the normal fluxes ...
Definition: dg_elements.h:50
unsigned required_nflux()
Set the number of flux components.
void calculate_element_averages(double *&average_value)
Calculate the averages in the element.
void fill_in_generic_residual_contribution_flux_transport(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
void fill_in_generic_residual_contribution_flux_transport(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
unsigned required_nflux()
Set the number of flux components.
void calculate_element_averages(double *&average_value)
Calculate the averages in the element.
General DGScalarAdvectionClass. Establish the template parameters.
FaceElement for Discontinuous Galerkin Problems.
DGScalarAdvectionFaceElement(FiniteElement *const &element_pt, const int &face_index)
Constructor.
unsigned required_nflux()
Set the number of flux components.
void numerical_flux(const Vector< double > &n_out, const Vector< double > &u_int, const Vector< double > &u_ext, Vector< double > &flux)
Calculate the normal flux, so this is the dot product of the numerical flux with n_in.
unsigned required_nflux()
Set the number of flux components.
void calculate_element_averages(double *&average_value)
Calculate the averages in the element.
void fill_in_generic_residual_contribution_flux_transport(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
void get_inverse_mass_matrix_times_residuals(Vector< double > &minv_res)
Function that returns the current value of the residuals multiplied by the inverse mass matrix (virtu...
unsigned required_nflux()
Set the number of flux components.
void calculate_element_averages(double *&average_value)
Calculate the averages in the element.
void fill_in_generic_residual_contribution_flux_transport(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
General DGScalarAdvectionClass. Establish the template parameters.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1271
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4735
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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 double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3325
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 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
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
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement....
Definition: elements.cc:5132
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
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
Base class for the flux transport equations templated by the dimension DIM. The equations that are so...
virtual void fill_in_generic_residual_contribution_flux_transport(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
void calculate_element_averages(double *&average_values)
Compute the average values of the fluxes.
void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: elements.h:3439
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
Non-spectral version of the classes.
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.
QScalarAdvectionElement(const QScalarAdvectionElement< DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
double dshape_and_dtest_eulerian_flux_transport(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.
double dshape_and_dtest_eulerian_at_knot_flux_transport(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....
QScalarAdvectionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
General QLegendreElement class.
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.
QSpectralScalarAdvectionElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
double dshape_and_dtest_eulerian_flux_transport(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.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
QSpectralScalarAdvectionElement(const QSpectralScalarAdvectionElement< DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
double dshape_and_dtest_eulerian_at_knot_flux_transport(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: matrices.h:1370
Base class for advection equations.
ScalarAdvectionWindFctPt Wind_fct_pt
Function pointer to the wind function.
unsigned required_nvalue(const unsigned &n) const
The number of unknowns at each node is the number of values.
ScalarAdvectionWindFctPt wind_fct_pt() const
Access function: Pointer to wind function. Const version.
ScalarAdvectionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
unsigned nflux() const
A single flux is interpolated.
virtual void get_wind_scalar_adv(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Return the wind at a given position.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt initial_condition_pt, const double &t, Vector< double > &error, Vector< double > &norm)
Compute the error and norm of solution integrated over the element Does not plot the error in the out...
void dflux_du(const Vector< double > &u, RankThreeTensor< double > &df_du)
Return the flux derivatives as a function of the unknowns.
void flux(const Vector< double > &u, DenseMatrix< double > &f)
Return the flux as a function of the unknown.
void(* ScalarAdvectionWindFctPt)(const Vector< double > &x, Vector< double > &wind)
Typedef for a wind function as a possible function of position.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
double dot(const Vector< double > &a, const Vector< double > &b)
Probably not always best/fastest because not optimised for dimension but useful...
Definition: Vector.h:291
//////////////////////////////////////////////////////////////////// ////////////////////////////////...