euler_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_EULER_ELEMENTS_HEADER
29#define OOMPH_EULER_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/Qspectral_elements.h"
38#include "../generic/dg_elements.h"
39
40namespace oomph
41{
42 //==============================================================
43 /// Base class for Euler equations
44 //=============================================================
45 template<unsigned DIM>
47 {
48 // Default value for gamma (initialised to 1.4)
49 static double Default_Gamma_Value;
50
51 // Pointer to the specific gas constant
52 double* Gamma_pt;
53
54 /// Storage for the average primitive values
56
57 /// Storage for the approximated gradients
59
60 protected:
61 /// DIM momentum-components, a density and an energy are transported
62 inline unsigned nflux() const
63 {
64 return DIM + 2;
65 }
66
67 /// Return the flux as a function of the unknown
68 void flux(const Vector<double>& u, DenseMatrix<double>& f);
69
70 /// Return the flux derivatives as a function of the unknowns
71 // void dflux_du(const Vector<double> &u, RankThreeTensor<double> &df_du);
72
73 public:
74 /// Constructor
79 {
80 // Set the default value of gamma
82 }
83
84 /// Destructor
86 {
87 if (this->Average_prim_value != 0)
88 {
89 delete[] Average_prim_value;
91 }
92 if (this->Average_gradient != 0)
93 {
94 delete[] Average_gradient;
96 }
97 }
98
99 /// Calculate the pressure from the unknowns
100 double pressure(const Vector<double>& u) const;
101
102 /// Return the value of gamma
103 double gamma() const
104 {
105 return *Gamma_pt;
106 }
107
108 /// Access function for the pointer to gamma
109 double*& gamma_pt()
110 {
111 return Gamma_pt;
112 }
113
114 /// Access function for the pointer to gamma (const version)
115 double* const& gamma_pt() const
116 {
117 return Gamma_pt;
118 }
119
120 /// The number of unknowns at each node is the number of flux
121 /// components
122 inline unsigned required_nvalue(const unsigned& n) const
123 {
124 return DIM + 2;
125 }
126
127 /// Compute the error and norm of solution integrated over the element
128 /// Does not plot the error in the outfile
130 std::ostream& outfile,
132 const double& t,
133 Vector<double>& error,
134 Vector<double>& norm)
135 {
136 // Find the number of fluxes
137 const unsigned n_flux = this->nflux();
138 // Find the number of nodes
139 const unsigned n_node = this->nnode();
140 // Storage for the shape function and derivatives of shape function
141 Shape psi(n_node);
142 DShape dpsidx(n_node, DIM);
143
144 // Find the number of integration points
145 unsigned n_intpt = this->integral_pt()->nweight();
146
147 error.resize(n_flux);
148 norm.resize(n_flux);
149 for (unsigned i = 0; i < n_flux; i++)
150 {
151 error[i] = 0.0;
152 norm[i] = 0.0;
153 }
154
155 Vector<double> s(DIM);
156
157 // Loop over the integration points
158 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
159 {
160 // Get the shape functions at the knot
161 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
162
163 // Get the integral weight
164 double W = this->integral_pt()->weight(ipt) * J;
165
166 // Storage for the local functions
168 Vector<double> interpolated_u(n_flux, 0.0);
169
170 // Loop over the shape functions
171 for (unsigned l = 0; l < n_node; l++)
172 {
173 // Locally cache the shape function
174 const double psi_ = psi(l);
175 for (unsigned i = 0; i < DIM; i++)
176 {
177 interpolated_x[i] += this->nodal_position(l, i) * psi_;
178 }
179
180 for (unsigned i = 0; i < n_flux; i++)
181 {
182 // Calculate the velocity and tangent vector
183 interpolated_u[i] += this->nodal_value(l, i) * psi_;
184 }
185 }
186
187 // Now get the exact solution at this value of x and t
188 Vector<double> exact_u(n_flux, 0.0);
189 (*exact_solution_pt)(t, interpolated_x, exact_u);
190
191 // Loop over the unknowns
192 for (unsigned i = 0; i < n_flux; i++)
193 {
194 error[i] += pow((interpolated_u[i] - exact_u[i]), 2.0) * W;
195 norm[i] += interpolated_u[i] * interpolated_u[i] * W;
196 }
197 }
198 }
199
200
201 /// Output function:
202 /// x,y,u or x,y,z,u
203 void output(std::ostream& outfile)
204 {
205 unsigned nplot = 5;
206 output(outfile, nplot);
207 }
208
209 /// Output function:
210 /// x,y,u or x,y,z,u at n_plot^DIM plot points
211 void output(std::ostream& outfile, const unsigned& n_plot);
212
214 {
215 // Find the number of fluxes
216 const unsigned n_flux = this->nflux();
217 // Resize the memory if necessary
218 if (this->Average_prim_value == 0)
219 {
220 this->Average_prim_value = new double[n_flux];
221 }
222 // Will move this one day
223 if (this->Average_gradient == 0)
224 {
225 this->Average_gradient = new double[n_flux * DIM];
226 }
227 }
228
229 /// Return access to the average gradient
230 double& average_gradient(const unsigned& i, const unsigned& j)
231 {
232 if (Average_gradient == 0)
233 {
234 throw OomphLibError("Averages not calculated yet",
235 OOMPH_CURRENT_FUNCTION,
236 OOMPH_EXCEPTION_LOCATION);
237 }
238 return Average_gradient[DIM * i + j];
239 }
240
241 /// Return the average values
242 double& average_prim_value(const unsigned& i)
243 {
244 if (Average_prim_value == 0)
245 {
246 throw OomphLibError("Averages not calculated yet",
247 OOMPH_CURRENT_FUNCTION,
248 OOMPH_EXCEPTION_LOCATION);
249 }
250 return Average_prim_value[i];
251 }
252 };
253
254
255 template<unsigned DIM, unsigned NNODE_1D>
256 class QSpectralEulerElement : public virtual QSpectralElement<DIM, NNODE_1D>,
257 public virtual EulerEquations<DIM>
258 {
259 public:
260 /// Constructor: Call constructors for QElement and
261 /// Advection Diffusion equations
263 : QSpectralElement<DIM, NNODE_1D>(), EulerEquations<DIM>()
264 {
265 }
266
267 /// Broken copy constructor
269 delete;
270
271 /// Broken assignment operator
272 // Commented out broken assignment operator because this can lead to a
273 // conflict warning when used in the virtual inheritence hierarchy.
274 // Essentially the compiler doesn't realise that two separate
275 // implementations of the broken function are the same and so, quite
276 // rightly, it shouts.
277 /*void operator=(
278 const QSpectralEulerElement<DIM,NNODE_1D>&) = delete;*/
279
280
281 /// Output function:
282 /// x,y,u or x,y,z,u
283 void output(std::ostream& outfile)
284 {
286 }
287
288 /// Output function:
289 /// x,y,u or x,y,z,u at n_plot^DIM plot points
290 void output(std::ostream& outfile, const unsigned& n_plot)
291 {
292 EulerEquations<DIM>::output(outfile, n_plot);
293 }
294
295
296 /*/// C-style output function:
297 /// x,y,u or x,y,z,u
298 void output(FILE* file_pt)
299 {
300 EulerEquations<NFLUX,DIM>::output(file_pt);
301 }
302
303 /// C-style output function:
304 /// x,y,u or x,y,z,u at n_plot^DIM plot points
305 void output(FILE* file_pt, const unsigned &n_plot)
306 {
307 EulerEquations<NFLUX,DIM>::output(file_pt,n_plot);
308 }
309
310 /// Output function for an exact solution:
311 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
312 void output_fct(std::ostream &outfile, const unsigned &n_plot,
313 FiniteElement::SteadyExactSolutionFctPt
314 exact_soln_pt)
315 {
316 EulerEquations<NFLUX,DIM>::
317 output_fct(outfile,n_plot,exact_soln_pt);}
318
319
320 /// Output function for a time-dependent exact solution.
321 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
322 /// (Calls the steady version)
323 void output_fct(std::ostream &outfile, const unsigned &n_plot,
324 const double& time,
325 FiniteElement::UnsteadyExactSolutionFctPt
326 exact_soln_pt)
327 {
328 EulerEquations<NFLUX,DIM>::
329 output_fct(outfile,n_plot,time,exact_soln_pt);
330 }*/
331
332
333 protected:
334 /// Shape, test functions & derivs. w.r.t. to global coords. Return
335 /// Jacobian.
337 const Vector<double>& s,
338 Shape& psi,
339 DShape& dpsidx,
340 Shape& test,
341 DShape& dtestdx) const;
342
343 /// Shape, test functions & derivs. w.r.t. to global coords. at
344 /// integration point ipt. Return Jacobian.
346 const unsigned& ipt,
347 Shape& psi,
348 DShape& dpsidx,
349 Shape& test,
350 DShape& dtestdx) const;
351 };
352
353 // Inline functions:
354
355
356 //======================================================================
357 /// Define the shape functions and test functions and derivatives
358 /// w.r.t. global coordinates and return Jacobian of mapping.
359 ///
360 /// Galerkin: Test functions = shape functions
361 //======================================================================
362 template<unsigned DIM, unsigned NNODE_1D>
365 Shape& psi,
366 DShape& dpsidx,
367 Shape& test,
368 DShape& dtestdx) const
369 {
370 // Call the geometrical shape functions and derivatives
371 double J = this->dshape_eulerian(s, psi, dpsidx);
372
373 // Loop over the test functions and derivatives and set them equal to the
374 // shape functions
375 for (unsigned i = 0; i < NNODE_1D; i++)
376 {
377 test[i] = psi[i];
378 for (unsigned j = 0; j < DIM; j++)
379 {
380 dtestdx(i, j) = dpsidx(i, j);
381 }
382 }
383
384 // Return the jacobian
385 return J;
386 }
387
388
389 //======================================================================
390 /// Define the shape functions and test functions and derivatives
391 /// w.r.t. global coordinates and return Jacobian of mapping.
392 ///
393 /// Galerkin: Test functions = shape functions
394 //======================================================================
395 template<unsigned DIM, unsigned NNODE_1D>
398 Shape& psi,
399 DShape& dpsidx,
400 Shape& test,
401 DShape& dtestdx) const
402 {
403 // Call the geometrical shape functions and derivatives
404 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
405
406 // Set the test functions equal to the shape functions (pointer copy)
407 test = psi;
408 dtestdx = dpsidx;
409
410 // Return the jacobian
411 return J;
412 }
413
414
415 /// /////////////////////////////////////////////////////////////////////
416 /// /////////////////////////////////////////////////////////////////////
417 /// /////////////////////////////////////////////////////////////////////
418
419
420 //=======================================================================
421 /// Face geometry for the QEulerElement elements:
422 /// The spatial dimension of the face elements is one lower than that
423 /// of the bulk element but they have the same number of points along
424 /// their 1D edges.
425 //=======================================================================
426 template<unsigned DIM, unsigned NNODE_1D>
428 : public virtual QSpectralElement<DIM - 1, NNODE_1D>
429 {
430 public:
431 /// Constructor: Call the constructor for the
432 /// appropriate lower-dimensional QElement
433 FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
434 };
435
436
437 //======================================================================
438 /// FaceElement for Discontinuous Galerkin Problems
439 //======================================================================
440 template<class ELEMENT>
441 class DGEulerFaceElement : public virtual FaceGeometry<ELEMENT>,
442 public virtual DGFaceElement
443 {
444 unsigned Nflux;
445
447
448 public:
449 /// Constructor
450 DGEulerFaceElement(FiniteElement* const& element_pt, const int& face_index)
451 : FaceGeometry<ELEMENT>(), DGFaceElement()
452 {
453 // Attach geometric information to the element
454 // N.B. This function also assigns nbulk_value from required_nvalue
455 // of the bulk element.
456 element_pt->build_face_element(face_index, this);
457 // Set the value of the flux
458 Nflux = 2 + element_pt->dim();
459 // Use the Face integration scheme of the element for 2D elements
460 if (Nflux > 3)
461 {
463 dynamic_cast<ELEMENT*>(element_pt)->face_integration_pt());
464 }
465 }
466
467 /// Specify the value of nodal zeta from the face geometry
468 /// The "global" intrinsic coordinate of the element when
469 /// viewed as part of a geometric object should be given by
470 /// the FaceElement representation, by default (needed to break
471 /// indeterminacy if bulk element is SolidElement)
472 double zeta_nodal(const unsigned& n,
473 const unsigned& k,
474 const unsigned& i) const
475 {
476 return FaceElement::zeta_nodal(n, k, i);
477 }
478
479
480 // There is a single required n_flux
481 unsigned required_nflux()
482 {
483 return Nflux;
484 }
485
486 /// Calculate the normal flux, so this is the dot product of the
487 /// numerical flux with n_out
489 const Vector<double>& u_int,
490 const Vector<double>& u_ext,
491 Vector<double>& flux)
492 {
493 // Let's follow the yellow book and use local Lax-Friedrichs
494 // This is almost certainly not the best flux to use ---
495 // further investigation is required here.
496 // Cache the bulk element
497 ELEMENT* cast_bulk_element_pt =
498 dynamic_cast<ELEMENT*>(this->bulk_element_pt());
499
500 // Find the dimension of the problem
501 const unsigned dim = cast_bulk_element_pt->dim();
502 // Find the number of variables
503 const unsigned n_flux = this->Nflux;
504
505 const double gamma = cast_bulk_element_pt->gamma();
506
507 // Now let's find the fluxes
508 DenseMatrix<double> flux_int(n_flux, dim);
509 DenseMatrix<double> flux_ext(n_flux, dim);
510
511 cast_bulk_element_pt->flux(u_int, flux_int);
512 cast_bulk_element_pt->flux(u_ext, flux_ext);
513
514 DenseMatrix<double> flux_av(n_flux, dim);
515 for (unsigned i = 0; i < n_flux; i++)
516 {
517 for (unsigned j = 0; j < dim; j++)
518 {
519 flux_av(i, j) = 0.5 * (flux_int(i, j) + flux_ext(i, j));
520 }
521 }
522
523 flux.initialise(0.0);
524 // Now set the first part of the numerical flux
525 for (unsigned i = 0; i < n_flux; i++)
526 {
527 for (unsigned j = 0; j < dim; j++)
528 {
529 flux[i] += flux_av(i, j) * n_out[j];
530 }
531 }
532
533 // Now let's find the normal jumps in the fluxes
534 Vector<double> jump(n_flux);
535 // Now we take the normal jump in the quantities
536 // The first two are scalars
537 for (unsigned i = 0; i < 2; i++)
538 {
539 jump[i] = u_int[i] - u_ext[i];
540 }
541
542 // The next are vectors so treat them accordingly ??
543 // By taking the component dotted with the normal component
544 double velocity_jump = 0.0;
545 for (unsigned j = 0; j < dim; j++)
546 {
547 velocity_jump += (u_int[2 + j] - u_ext[2 + j]) * n_out[j];
548 }
549 // Now we multiply by the outer normal to get the vector flux
550 for (unsigned j = 0; j < dim; j++)
551 {
552 jump[2 + j] = velocity_jump * n_out[j];
553 }
554
555
556 // Let's find the Roe average
557 /* Vector<double> u_average(n_flux);
558 double sum = sqrt(u_int[0]) + sqrt(u_ext[0]);
559 for(unsigned i=0;i<n_flux;i++)
560 {
561 u_average[i] = (sqrt(u_int[0])*u_int[i] + sqrt(u_ext[0])*u_ext[i])/sum;
562 }
563
564 //Find the average pressure
565 double p = cast_bulk_element_pt->pressure(u_average);
566
567 //Now do the normal velocity
568 double vel = 0.0;
569 for(unsigned j=0;j<dim;j++)
570 {
571 //vel += u_average[2+j]*u_average[2+j];
572 vel += u_average[2+j]*n_out[j];
573 }
574 //vel = sqrt(vel)/u_average[0];
575 vel /= u_average[0];
576
577 double arg = std::fabs(gamma*p/u_average[0]);*/
578
579 // Let's calculate the internal and external pressures and then enthalpies
580 double p_int = cast_bulk_element_pt->pressure(u_int);
581 double p_ext = cast_bulk_element_pt->pressure(u_ext);
582
583 // Limit the pressures to zero if necessary, but keep the energy the
584 // same
585 if (p_int < 0)
586 {
587 oomph_info << "Negative int pressure" << p_int << "\n";
588 p_int = 0.0;
589 }
590 if (p_ext < 0)
591 {
592 oomph_info << "Negative ext pressure " << p_ext << "\n";
593 p_ext = 0.0;
594 }
595
596 double H_int = (u_int[1] + p_int) / u_int[0];
597 double H_ext = (u_ext[1] + p_ext) / u_ext[0];
598
599 // Also calculate the velocities
600 Vector<double> vel_int(2), vel_ext(2);
601 for (unsigned j = 0; j < dim; j++)
602 {
603 vel_int[j] = u_int[2 + j] / u_int[0];
604 vel_ext[j] = u_ext[2 + j] / u_ext[0];
605 }
606
607 // Now we calculate the Roe averaged values
608 Vector<double> vel_average(dim);
609 double s_int = sqrt(u_int[0]);
610 double s_ext = sqrt(u_ext[0]);
611 double sum = s_int + s_ext;
612
613 // Velocities
614 for (unsigned j = 0; j < dim; j++)
615 {
616 vel_average[j] = (s_int * vel_int[j] + s_ext * vel_ext[j]) / sum;
617 }
618
619 // The enthalpy
620 double H_average = (s_int * H_int + s_ext * H_ext) / sum;
621
622 // Now calculate the square of the sound speed
623 double arg = H_average;
624 for (unsigned j = 0; j < dim; j++)
625 {
626 arg -= 0.5 * vel_average[j];
627 }
628 arg *= (gamma - 1.0);
629 // Get the local sound speed
630 if (arg < 0.0)
631 {
632 oomph_info << "Square of sound speed is negative!\n";
633 arg = 0.0;
634 }
635 double a = sqrt(arg);
636
637 // Calculate the normal average velocity
638 // Now do the normal velocity
639 double vel = 0.0;
640 for (unsigned j = 0; j < dim; j++)
641 {
642 vel += vel_average[j] * n_out[j];
643 }
644
645 Vector<double> eigA(3, 0.0);
646
647 eigA[0] = vel - a;
648 eigA[1] = vel;
649 eigA[2] = vel + a;
650
651 double lambda = std::fabs(eigA[0]);
652 for (unsigned i = 1; i < 3; i++)
653 {
654 if (std::fabs(eigA[i]) > lambda)
655 {
656 lambda = std::fabs(eigA[i]);
657 }
658 }
659
660
661 // Find the magnitude of the external velocity
662 /*double vel_mag = 0.0;
663 for(unsigned j=0;j<dim;j++) {vel_mag += u_ext[2+j]*u_ext[2+j];}
664 vel_mag = sqrt(vel_mag)/u_ext[0];
665
666 //Get the pressure
667 double p = cast_bulk_element_pt->pressure(u_ext);
668 double lambda_ext = vel_mag + sqrt(std::fabs(gamma*p/u_ext[0]));
669
670 //Let's do the same for the internal one
671 vel_mag = 0.0;
672 for(unsigned j=0;j<dim;j++) {vel_mag += u_int[2+j]*u_int[2+j];}
673 vel_mag = sqrt(vel_mag)/u_int[0];
674
675 //Get the pressure
676 p = cast_bulk_element_pt->pressure(u_int);
677 double lambda_int = vel_mag + sqrt(std::fabs(gamma*p/u_int[0]));
678
679 //Now take the largest
680 double lambda = lambda_int;
681 if(lambda_ext > lambda) {lambda = lambda_ext;}*/
682
683 for (unsigned i = 0; i < n_flux; i++)
684 {
685 flux[i] += 0.5 * lambda * jump[i];
686 }
687 }
688 };
689
690
691 //======================================================================
692 /// FaceElement for Discontinuous Galerkin Problems with reflection
693 /// boundary conditions
694 //======================================================================
695 template<class ELEMENT>
696 class DGEulerFaceReflectionElement : public virtual FaceGeometry<ELEMENT>,
697 public virtual DGFaceElement
698 {
699 unsigned Nflux;
700
701 public:
702 /// Constructor
704 const int& face_index)
705 : FaceGeometry<ELEMENT>(), DGFaceElement()
706 {
707 // Attach geometric information to the element
708 // N.B. This function also assigns nbulk_value from required_nvalue
709 // of the bulk element.
710 element_pt->build_face_element(face_index, this);
711 // Set the value of the flux
712 Nflux = 2 + element_pt->dim();
713 }
714
715
716 /// Specify the value of nodal zeta from the face geometry
717 /// The "global" intrinsic coordinate of the element when
718 /// viewed as part of a geometric object should be given by
719 /// the FaceElement representation, by default (needed to break
720 /// indeterminacy if bulk element is SolidElement)
721 double zeta_nodal(const unsigned& n,
722 const unsigned& k,
723 const unsigned& i) const
724 {
725 return FaceElement::zeta_nodal(n, k, i);
726 }
727
728 // There is a single required n_flux
729 unsigned required_nflux()
730 {
731 return Nflux;
732 }
733
734 /// We overload interpolated_u to reflect
736 {
737 // Get the standard interpolated_u
739
740 // Now do the reflection condition for the velocities
741 // Find dot product of normal and velocities
742 const unsigned nodal_dim = this->nodal_dimension();
743 Vector<double> n(nodal_dim);
744 this->outer_unit_normal(s, n);
745
746 double dot = 0.0;
747 for (unsigned j = 0; j < nodal_dim; j++)
748 {
749 dot += n[j] * u[2 + j];
750 }
751
752 // Now subtract
753 for (unsigned j = 0; j < nodal_dim; j++)
754 {
755 u[2 + j] -= 2.0 * dot * n[j];
756 }
757 }
758 };
759
760
761 //=================================================================
762 /// General DGEulerClass. Establish the template parameters
763 //===================================================================
764 template<unsigned DIM, unsigned NNODE_1D>
766 {
767 };
768
769
770 //==================================================================
771 // Specialization for 1D DG Advection element
772 //==================================================================
773 template<unsigned NNODE_1D>
774 class DGSpectralEulerElement<1, NNODE_1D>
775 : public QSpectralEulerElement<1, NNODE_1D>, public DGElement
776 {
777 friend class DGEulerFaceElement<DGSpectralEulerElement<1, NNODE_1D>>;
778
780
781 public:
782 /// Overload the required number of fluxes for the DGElement
783 unsigned required_nflux()
784 {
785 return this->nflux();
786 }
787
788 // Calculate averages
789 void calculate_element_averages(double*& average_value)
790 {
792 }
793
794
795 // Constructor
797 {
798 // Need to up the order of integration for the accurate resolution
799 // of the quadratic non-linearities
800 this->set_integration_scheme(new Gauss<1, NNODE_1D>);
801 }
802
803 // Dummy
805 {
806 return 0;
807 }
808
809
811
813 {
814 // Make the two faces
815 Face_element_pt.resize(2);
816 // Make the face on the left
817 Face_element_pt[0] =
819 // Make the face on the right
820 Face_element_pt[1] =
822 }
823
824
825 /// Compute the residuals for the Navier--Stokes equations;
826 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
828 Vector<double>& residuals,
829 DenseMatrix<double>& jacobian,
830 DenseMatrix<double>& mass_matrix,
831 unsigned flag)
832 {
835 residuals, jacobian, mass_matrix, flag);
836
837 this->add_flux_contributions_to_residuals(residuals, jacobian, flag);
838 }
839
840
841 //============================================================================
842 /// Function that returns the current value of the residuals
843 /// multiplied by the inverse mass matrix (virtual so that it can be
844 /// overloaded specific elements in which time and memory saving tricks can
845 /// be applied)
846 //============================================================================
847 /* void get_inverse_mass_matrix_times_residuals(Vector<double> &minv_res)
848 {
849 //Now let's assemble stuff
850 const unsigned n_dof = this->ndof();
851
852 //Resize and initialise the vector that will holds the residuals
853 minv_res.resize(n_dof);
854 for(unsigned n=0;n<n_dof;n++) {minv_res[n] = 0.0;}
855
856 //If we are recycling the mass matrix
857 if(Mass_matrix_reuse_is_enabled && Mass_matrix_has_been_computed)
858 {
859 //Get the residuals
860 this->fill_in_contribution_to_residuals(minv_res);
861 }
862 //Otherwise
863 else
864 {
865 //Temporary mass matrix
866 DenseDoubleMatrix M(n_dof,n_dof,0.0);
867
868 //Get the local mass matrix and residuals
869 this->fill_in_contribution_to_mass_matrix(minv_res,M);
870
871 //Store the diagonal entries
872 Inverse_mass_diagonal.clear();
873 for(unsigned n=0;n<n_dof;n++)
874 {Inverse_mass_diagonal.push_back(1.0/M(n,n));}
875
876 //The mass matrix has been computed
877 Mass_matrix_has_been_computed=true;
878 }
879
880 for(unsigned n=0;n<n_dof;n++) {minv_res[n] *= Inverse_mass_diagonal[n];}
881
882 }*/
883 };
884
885
886 //=======================================================================
887 /// Face geometry of the 1D DG elements
888 //=======================================================================
889 template<unsigned NNODE_1D>
891 : public virtual PointElement
892 {
893 public:
895 };
896
897
898 //==================================================================
899 /// Specialisation for 2D DG Elements
900 //==================================================================
901 template<unsigned NNODE_1D>
902 class DGSpectralEulerElement<2, NNODE_1D>
903 : public QSpectralEulerElement<2, NNODE_1D>, public DGElement
904 {
905 friend class DGEulerFaceElement<DGSpectralEulerElement<2, NNODE_1D>>;
906
908
909 public:
910 /// Overload the required number of fluxes for the DGElement
911 unsigned required_nflux()
912 {
913 return this->nflux();
914 }
915
916 // Calculate averages
917 void calculate_element_averages(double*& average_value)
918 {
920 }
921
922 // Constructor
924 {
925 // Need to up the order of integration for the accurate resolution
926 // of the quadratic non-linearities
927 this->set_integration_scheme(
929 }
930
932
934 {
935 return &Default_face_integration_scheme;
936 }
937
939 {
940 Face_element_pt.resize(4);
941 Face_element_pt[0] =
943 Face_element_pt[1] =
945 Face_element_pt[2] =
947 Face_element_pt[3] =
949 }
950
951
952 /// Compute the residuals for the Navier--Stokes equations;
953 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
955 Vector<double>& residuals,
956 DenseMatrix<double>& jacobian,
957 DenseMatrix<double>& mass_matrix,
958 unsigned flag)
959 {
962 residuals, jacobian, mass_matrix, flag);
963
964 this->add_flux_contributions_to_residuals(residuals, jacobian, flag);
965 }
966 };
967
968
969 //=======================================================================
970 /// Face geometry of the DG elements
971 //=======================================================================
972 template<unsigned NNODE_1D>
974 : public virtual QSpectralElement<1, NNODE_1D>
975 {
976 public:
977 FaceGeometry() : QSpectralElement<1, NNODE_1D>() {}
978 };
979
980
981} // namespace oomph
982
983#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
FaceElement for Discontinuous Galerkin Problems.
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_out.
unsigned required_nflux()
Set the number of flux components.
DGEulerFaceElement(FiniteElement *const &element_pt, const int &face_index)
Constructor.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
FaceElement for Discontinuous Galerkin Problems with reflection boundary conditions.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
void interpolated_u(const Vector< double > &s, Vector< double > &u)
We overload interpolated_u to reflect.
unsigned required_nflux()
Set the number of flux components.
DGEulerFaceReflectionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor.
Base class for Discontinuous Galerkin Faces. These are responsible for calculating the normal fluxes ...
Definition: dg_elements.h:50
virtual void interpolated_u(const Vector< double > &s, Vector< double > &f)
Return the interpolated values of the unknown fluxes.
Definition: dg_elements.cc:196
void calculate_element_averages(double *&average_value)
Calculate the averages in the element.
unsigned required_nflux()
Overload the required number of fluxes for the DGElement.
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()
Overload the required number of fluxes for the DGElement.
void calculate_element_averages(double *&average_value)
Calculate the averages in the element.
static Gauss< 1, NNODE_1D > Default_face_integration_scheme
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 DGEulerClass. 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
Base class for Euler equations.
virtual ~EulerEquations()
Destructor.
double & average_prim_value(const unsigned &i)
Return the average values.
double pressure(const Vector< double > &u) const
Calculate the pressure from the unknowns.
void flux(const Vector< double > &u, DenseMatrix< double > &f)
Return the flux as a function of the unknown.
double * Average_prim_value
Storage for the average primitive values.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_solution_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 output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
double gamma() const
Return the value of gamma.
double & average_gradient(const unsigned &i, const unsigned &j)
Return access to the average gradient.
unsigned nflux() const
DIM momentum-components, a density and an energy are transported.
double * Average_gradient
Storage for the approximated gradients.
double *& gamma_pt()
Access function for the pointer to gamma.
unsigned required_nvalue(const unsigned &n) const
The number of unknowns at each node is the number of flux components.
double *const & gamma_pt() const
Access function for the pointer to gamma (const version)
static double Default_Gamma_Value
EulerEquations()
Return the flux derivatives as a function of the unknowns.
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:6006
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary,...
Definition: elements.h:4497
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
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3210
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.
Class for multidimensional Gauss Lobatto Legendre integration rules empty - just establishes template...
Definition: integral.h:1281
Class for multidimensional Gaussian integration rules.
Definition: integral.h:145
Generic class for numerical integration schemes:
Definition: integral.h:49
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
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.
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....
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.
QSpectralEulerElement()
Constructor: Call constructors for QElement and Advection Diffusion equations.
QSpectralEulerElement(const QSpectralEulerElement< DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
void output(std::ostream &outfile)
Broken assignment operator.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
void initialise(const _Tp &__value)
Iterate over all values and set to the desired value.
Definition: Vector.h:167
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...