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-2023 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 
40 namespace 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
76  : FluxTransportEquations<DIM>(),
79  {
80  // Set the default value of gamma
82  }
83 
84  /// Destructor
85  virtual ~EulerEquations()
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;
95  Average_gradient = 0;
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
167  Vector<double> interpolated_x(DIM, 0.0);
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>
427  class FaceGeometry<QSpectralEulerElement<DIM, 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
488  void numerical_flux(const Vector<double>& 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_gradient(const unsigned &i, const unsigned &j)
Return access to the average gradient.
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.
unsigned nflux() const
DIM momentum-components, a density and an energy are transported.
double * Average_gradient
Storage for the approximated gradients.
unsigned required_nvalue(const unsigned &n) const
The number of unknowns at each node is the number of flux components.
double & average_prim_value(const unsigned &i)
Return the average values.
static double Default_Gamma_Value
double *const & gamma_pt() const
Access function for the pointer to gamma (const version)
EulerEquations()
Return the flux derivatives as a function of the unknowns.
double *& gamma_pt()
Access function for the pointer to gamma.
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
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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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
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...