scalar_advection_elements.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2025 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Header file for scalar advection elements
 
#ifndef OOMPH_SCALAR_ADVECTION_ELEMENTS_HEADER
#define OOMPH_SCALAR_ADVECTION_ELEMENTS_HEADER
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "flux_transport_elements.h"
#include "generic/Qelements.h"
#include "generic/Qspectral_elements.h"
#include "generic/dg_elements.h"
 
namespace oomph
{
  //==============================================================
  /// Base class for advection equations
  //=============================================================
  template<unsigned DIM>
  class ScalarAdvectionEquations : public FluxTransportEquations<DIM>
  {
    /// Typedef for a wind function as a possible function of position
    typedef void (*ScalarAdvectionWindFctPt)(const Vector<double>& x,
                                             Vector<double>& wind);
 
    /// Function pointer to the wind function
    ScalarAdvectionWindFctPt Wind_fct_pt;
 
  protected:
    /// A single flux is interpolated
    inline unsigned nflux() const
    {
      return 1;
    }
 
    /// Return the flux as a function of the unknown
    void flux(const Vector<double>& u, DenseMatrix<double>& f);
 
    /// Return the flux derivatives as a function of the unknowns
    void dflux_du(const Vector<double>& u, RankThreeTensor<double>& df_du);
 
    /// Return the wind at a given position
    inline virtual void get_wind_scalar_adv(const unsigned& ipt,
                                            const Vector<double>& s,
                                            const Vector<double>& x,
                                            Vector<double>& wind) const
    {
      // If no wind function has been set, return zero
      if (Wind_fct_pt == 0)
      {
        for (unsigned i = 0; i < DIM; i++)
        {
          wind[i] = 0.0;
        }
      }
      else
      {
        // Get wind
        (*Wind_fct_pt)(x, wind);
      }
    }
 
  public:
    /// Constructor
    ScalarAdvectionEquations() : FluxTransportEquations<DIM>(), Wind_fct_pt(0)
    {
    }
 
    /// Access function: Pointer to wind function
    ScalarAdvectionWindFctPt& wind_fct_pt()
    {
      return Wind_fct_pt;
    }
 
    /// Access function: Pointer to wind function. Const version
    ScalarAdvectionWindFctPt wind_fct_pt() const
    {
      return Wind_fct_pt;
    }
 
    /// The number of unknowns at each node is the number of values
    unsigned required_nvalue(const unsigned& n) const
    {
      return 1;
    }
 
    /// Compute the error and norm of solution integrated over the element
    /// Does not plot the error in the outfile
    void compute_error(
      std::ostream& outfile,
      FiniteElement::UnsteadyExactSolutionFctPt initial_condition_pt,
      const double& t,
      Vector<double>& error,
      Vector<double>& norm)
    {
      // Find the number of fluxes
      const unsigned n_flux = this->nflux();
      // Find the number of nodes
      const unsigned n_node = this->nnode();
      // Storage for the shape function and derivatives of shape function
      Shape psi(n_node);
      DShape dpsidx(n_node, DIM);
 
      // Find the number of integration points
      unsigned n_intpt = this->integral_pt()->nweight();
 
      error.resize(n_flux);
      norm.resize(n_flux);
      for (unsigned i = 0; i < n_flux; i++)
      {
        error[i] = 0.0;
        norm[i] = 0.0;
      }
 
      Vector<double> s(DIM);
 
      // Loop over the integration points
      for (unsigned ipt = 0; ipt < n_intpt; ipt++)
      {
        // Get the shape functions at the knot
        double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
 
        // Get the integral weight
        double W = this->integral_pt()->weight(ipt) * J;
 
        // Storage for the local functions
        Vector<double> interpolated_x(DIM, 0.0);
        Vector<double> interpolated_u(n_flux, 0.0);
 
        // Loop over the shape functions
        for (unsigned l = 0; l < n_node; l++)
        {
          // Locally cache the shape function
          const double psi_ = psi(l);
          for (unsigned i = 0; i < DIM; i++)
          {
            interpolated_x[i] += this->nodal_position(l, i) * psi_;
          }
 
          for (unsigned i = 0; i < n_flux; i++)
          {
            // Calculate the velocity and tangent vector
            interpolated_u[i] += this->nodal_value(l, i) * psi_;
          }
        }
 
        // Get the global wind
        Vector<double> wind(DIM);
        this->get_wind_scalar_adv(ipt, s, interpolated_x, wind);
 
        // Rescale the position
        for (unsigned i = 0; i < DIM; i++)
        {
          interpolated_x[i] -= wind[i] * t;
        }
 
        // Now get the initial condition at this value of x
        Vector<double> exact_u(n_flux, 0.0);
        (*initial_condition_pt)(0.0, interpolated_x, exact_u);
 
        // Loop over the unknowns
        for (unsigned i = 0; i < n_flux; i++)
        {
          error[i] += pow((interpolated_u[i] - exact_u[i]), 2.0) * W;
          norm[i] += interpolated_u[i] * interpolated_u[i] * W;
        }
      }
    }
  };
 
 
  template<unsigned DIM, unsigned NNODE_1D>
  class QSpectralScalarAdvectionElement
    : public virtual QSpectralElement<DIM, NNODE_1D>,
      public virtual ScalarAdvectionEquations<DIM>
  {
  public:
    /// Constructor: Call constructors for QElement and
    /// Advection Diffusion equations
    QSpectralScalarAdvectionElement()
      : QSpectralElement<DIM, NNODE_1D>(), ScalarAdvectionEquations<DIM>()
    {
    }
 
    /// Broken copy constructor
    QSpectralScalarAdvectionElement(
      const QSpectralScalarAdvectionElement<DIM, NNODE_1D>& dummy) = delete;
 
    /// Broken assignment operator
    // Commented out broken assignment operator because this can lead to a
    // conflict warning when used in the virtual inheritence hierarchy.
    // Essentially the compiler doesn't realise that two separate
    // implementations of the broken function are the same and so, quite
    // rightly, it shouts.
    /*void operator=(
     const QSpectralScalarAdvectionElement<DIM,NNODE_1D>&) = delete;*/
 
    ///  Required  # of `values' (pinned or dofs)
    /// at node n
    inline unsigned required_nvalue(const unsigned& n) const
    {
      return 1;
    }
 
    /// Output function:
    ///  x,y,u   or    x,y,z,u
    void output(std::ostream& outfile)
    {
      ScalarAdvectionEquations<DIM>::output(outfile);
    }
 
    /// Output function:
    ///  x,y,u   or    x,y,z,u at n_plot^DIM plot points
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      ScalarAdvectionEquations<DIM>::output(outfile, n_plot);
    }
 
 
    /*/// C-style output function:
   ///  x,y,u   or    x,y,z,u
   void output(FILE* file_pt)
    {
     ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt);
    }
 
   ///  C-style output function:
   ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
   void output(FILE* file_pt, const unsigned &n_plot)
    {
     ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt,n_plot);
    }
 
   /// Output function for an exact solution:
   ///  x,y,u_exact   or    x,y,z,u_exact at n_plot^DIM plot points
   void output_fct(std::ostream &outfile, const unsigned &n_plot,
                   FiniteElement::SteadyExactSolutionFctPt
                   exact_soln_pt)
    {
     ScalarAdvectionEquations<NFLUX,DIM>::
      output_fct(outfile,n_plot,exact_soln_pt);}
 
 
   /// Output function for a time-dependent exact solution.
   ///  x,y,u_exact   or    x,y,z,u_exact at n_plot^DIM plot points
   /// (Calls the steady version)
   void output_fct(std::ostream &outfile, const unsigned &n_plot,
                   const double& time,
                   FiniteElement::UnsteadyExactSolutionFctPt
                   exact_soln_pt)
    {
     ScalarAdvectionEquations<NFLUX,DIM>::
      output_fct(outfile,n_plot,time,exact_soln_pt);
      }*/
 
 
  protected:
    /// Shape, test functions & derivs. w.r.t. to global coords. Return
    /// Jacobian.
    inline 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. at
    /// integration point ipt. Return Jacobian.
    inline double dshape_and_dtest_eulerian_at_knot_flux_transport(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsidx,
      Shape& test,
      DShape& dtestdx) const;
  };
 
  // Inline functions:
 
 
  //======================================================================
  /// Define the shape functions and test functions and derivatives
  /// w.r.t. global coordinates and return Jacobian of mapping.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double QSpectralScalarAdvectionElement<DIM, NNODE_1D>::
    dshape_and_dtest_eulerian_flux_transport(const Vector<double>& s,
                                             Shape& psi,
                                             DShape& dpsidx,
                                             Shape& test,
                                             DShape& dtestdx) const
  {
    // Call the geometrical shape functions and derivatives
    double J = this->dshape_eulerian(s, psi, dpsidx);
 
    // Loop over the test functions and derivatives and set them equal to the
    // shape functions
    for (unsigned i = 0; i < NNODE_1D; i++)
    {
      test[i] = psi[i];
      for (unsigned j = 0; j < DIM; j++)
      {
        dtestdx(i, j) = dpsidx(i, j);
      }
    }
 
    // Return the jacobian
    return J;
  }
 
 
  //======================================================================
  /// Define the shape functions and test functions and derivatives
  /// w.r.t. global coordinates and return Jacobian of mapping.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double QSpectralScalarAdvectionElement<DIM, NNODE_1D>::
    dshape_and_dtest_eulerian_at_knot_flux_transport(const unsigned& ipt,
                                                     Shape& psi,
                                                     DShape& dpsidx,
                                                     Shape& test,
                                                     DShape& dtestdx) const
  {
    // Call the geometrical shape functions and derivatives
    double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
 
    // Set the test functions equal to the shape functions (pointer copy)
    test = psi;
    dtestdx = dpsidx;
 
    // Return the jacobian
    return J;
  }
 
 
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// Face geometry for the QScalarAdvectionElement elements:
  /// The spatial dimension of the face elements is one lower than that
  /// of the bulk element but they have the same number of points along
  /// their 1D edges.
  //=======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class FaceGeometry<QSpectralScalarAdvectionElement<DIM, NNODE_1D>>
    : public virtual QSpectralElement<DIM - 1, NNODE_1D>
  {
  public:
    /// Constructor: Call the constructor for the
    /// appropriate lower-dimensional QElement
    FaceGeometry() : QSpectralElement<DIM - 1, NNODE_1D>() {}
  };
 
 
  //======================================================================
  /// FaceElement for Discontinuous Galerkin Problems
  //======================================================================
  template<class ELEMENT>
  class DGScalarAdvectionFaceElement : public virtual FaceGeometry<ELEMENT>,
                                       public virtual DGFaceElement
  {
  public:
    /// Constructor
    DGScalarAdvectionFaceElement(FiniteElement* const& element_pt,
                                 const int& face_index)
      : FaceGeometry<ELEMENT>(), DGFaceElement()
    {
      // Attach geometric information to the element
      // N.B. This function also assigns nbulk_value from required_nvalue
      // of the bulk element.
      element_pt->build_face_element(face_index, this);
    }
 
 
    // There is a single required n_flux
    unsigned required_nflux()
    {
      return 1;
    }
 
    /// Calculate the normal flux, so this is the dot product of the
    /// numerical flux with n_in
    void numerical_flux(const Vector<double>& n_out,
                        const Vector<double>& u_int,
                        const Vector<double>& u_ext,
                        Vector<double>& flux)
    {
      const unsigned dim = this->nodal_dimension();
      Vector<double> Wind(dim);
      Vector<double> s, x;
      // Dummy integration point
      unsigned ipt = 0;
      dynamic_cast<ELEMENT*>(this->bulk_element_pt())
        ->get_wind_scalar_adv(ipt, s, x, Wind);
 
      // Now we can work this out for standard upwind advection
      double dot = 0.0;
      for (unsigned i = 0; i < dim; i++)
      {
        dot += Wind[i] * n_out[i];
      }
 
      const unsigned n_value = this->required_nflux();
      if (dot >= 0.0)
      {
        for (unsigned n = 0; n < n_value; n++)
        {
          flux[n] = dot * u_int[n];
        }
      }
      else
      {
        for (unsigned n = 0; n < n_value; n++)
        {
          flux[n] = dot * u_ext[n];
        }
      }
 
      // flux[0] = 0.5*(u_int[0]+u_ext[0])*n_out[0];
    }
  };
 
  //=================================================================
  /// General DGScalarAdvectionClass. Establish the template parameters
  //===================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class DGSpectralScalarAdvectionElement
  {
  };
 
 
  //==================================================================
  // Specialization for 1D DG Advection element
  //==================================================================
  template<unsigned NNODE_1D>
  class DGSpectralScalarAdvectionElement<1, NNODE_1D>
    : public QSpectralScalarAdvectionElement<1, NNODE_1D>, public DGElement
  {
    friend class DGScalarAdvectionFaceElement<
      DGSpectralScalarAdvectionElement<1, NNODE_1D>>;
 
    Vector<double> Inverse_mass_diagonal;
 
  public:
    // There is a single required n_flux
    unsigned required_nflux()
    {
      return 1;
    }
 
    // Calculate averages
    void calculate_element_averages(double*& average_value)
    {
      FluxTransportEquations<1>::calculate_element_averages(average_value);
    }
 
    // Constructor
    DGSpectralScalarAdvectionElement()
      : QSpectralScalarAdvectionElement<1, NNODE_1D>(), DGElement()
    {
    }
 
    ~DGSpectralScalarAdvectionElement() {}
 
    void build_all_faces()
    {
      // Make the two faces
      Face_element_pt.resize(2);
      // Make the face on the left
      Face_element_pt[0] = new DGScalarAdvectionFaceElement<
        DGSpectralScalarAdvectionElement<1, NNODE_1D>>(this, -1);
      // Make the face on the right
      Face_element_pt[1] = new DGScalarAdvectionFaceElement<
        DGSpectralScalarAdvectionElement<1, NNODE_1D>>(this, +1);
    }
 
 
    /// Compute the residuals for the Navier--Stokes equations;
    /// flag=1(or 0): do (or don't) compute the Jacobian as well.
    void fill_in_generic_residual_contribution_flux_transport(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix,
      unsigned flag)
    {
      QSpectralScalarAdvectionElement<1, NNODE_1D>::
        fill_in_generic_residual_contribution_flux_transport(
          residuals, jacobian, mass_matrix, flag);
 
      this->add_flux_contributions_to_residuals(residuals, jacobian, flag);
    }
 
 
    //============================================================================
    /// Function that returns the current value of the residuals
    /// multiplied by the inverse mass matrix (virtual so that it can be
    /// overloaded specific elements in which time and memory saving tricks can
    /// be applied)
    //============================================================================
    void get_inverse_mass_matrix_times_residuals(Vector<double>& minv_res)
    {
      // If there are external data this is not going to work
      if (nexternal_data() > 0)
      {
        std::ostringstream error_stream;
        error_stream
          << "Cannot use a discontinuous formulation for the mass matrix when\n"
          << "there are external data.\n "
          << "Do not call Problem::enable_discontinuous_formulation()\n";
 
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
 
      // Now let's assemble stuff
      const unsigned n_dof = this->ndof();
 
      // Resize and initialise the vector that will holds the residuals
      minv_res.resize(n_dof);
      for (unsigned n = 0; n < n_dof; n++)
      {
        minv_res[n] = 0.0;
      }
 
      // If we are recycling the mass matrix
      if (Mass_matrix_reuse_is_enabled && Mass_matrix_has_been_computed)
      {
        // Get the residuals
        this->fill_in_contribution_to_residuals(minv_res);
      }
      // Otherwise
      else
      {
        // Temporary mass matrix
        DenseDoubleMatrix M(n_dof, n_dof, 0.0);
 
        // Get the local mass matrix and residuals
        this->fill_in_contribution_to_mass_matrix(minv_res, M);
 
        // Store the diagonal entries
        Inverse_mass_diagonal.clear();
        for (unsigned n = 0; n < n_dof; n++)
        {
          Inverse_mass_diagonal.push_back(1.0 / M(n, n));
        }
 
        // The mass matrix has been computed
        Mass_matrix_has_been_computed = true;
      }
 
      for (unsigned n = 0; n < n_dof; n++)
      {
        minv_res[n] *= Inverse_mass_diagonal[n];
      }
    }
  };
 
 
  //=======================================================================
  /// Face geometry of the 1D  DG elements
  //=======================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<DGSpectralScalarAdvectionElement<1, NNODE_1D>>
    : public virtual PointElement
  {
  public:
    FaceGeometry() : PointElement() {}
  };
 
 
  //==================================================================
  /// Specialisation for 2D DG Elements
  //==================================================================
  template<unsigned NNODE_1D>
  class DGSpectralScalarAdvectionElement<2, NNODE_1D>
    : public QSpectralScalarAdvectionElement<2, NNODE_1D>, public DGElement
  {
    friend class DGScalarAdvectionFaceElement<
      DGSpectralScalarAdvectionElement<2, NNODE_1D>>;
 
  public:
    // Calculate averages
    void calculate_element_averages(double*& average_value)
    {
      FluxTransportEquations<2>::calculate_element_averages(average_value);
    }
 
    // There is a single required n_flux
    unsigned required_nflux()
    {
      return 1;
    }
 
    // Constructor
    DGSpectralScalarAdvectionElement()
      : QSpectralScalarAdvectionElement<2, NNODE_1D>(), DGElement()
    {
    }
 
    ~DGSpectralScalarAdvectionElement() {}
 
    void build_all_faces()
    {
      Face_element_pt.resize(4);
      Face_element_pt[0] = new DGScalarAdvectionFaceElement<
        DGSpectralScalarAdvectionElement<2, NNODE_1D>>(this, 2);
      Face_element_pt[1] = new DGScalarAdvectionFaceElement<
        DGSpectralScalarAdvectionElement<2, NNODE_1D>>(this, 1);
      Face_element_pt[2] = new DGScalarAdvectionFaceElement<
        DGSpectralScalarAdvectionElement<2, NNODE_1D>>(this, -2);
      Face_element_pt[3] = new DGScalarAdvectionFaceElement<
        DGSpectralScalarAdvectionElement<2, NNODE_1D>>(this, -1);
    }
 
 
    /// Compute the residuals for the Navier--Stokes equations;
    /// flag=1(or 0): do (or don't) compute the Jacobian as well.
    void fill_in_generic_residual_contribution_flux_transport(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix,
      unsigned flag)
    {
      QSpectralScalarAdvectionElement<2, NNODE_1D>::
        fill_in_generic_residual_contribution_flux_transport(
          residuals, jacobian, mass_matrix, flag);
 
      this->add_flux_contributions_to_residuals(residuals, jacobian, flag);
    }
  };
 
 
  //=======================================================================
  /// Face geometry of the DG elements
  //=======================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<DGSpectralScalarAdvectionElement<2, NNODE_1D>>
    : public virtual QSpectralElement<1, NNODE_1D>
  {
  public:
    FaceGeometry() : QSpectralElement<1, NNODE_1D>() {}
  };
 
 
  //=============================================================
  /// Non-spectral version of the classes
  //============================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class QScalarAdvectionElement : public virtual QElement<DIM, NNODE_1D>,
                                  public virtual ScalarAdvectionEquations<DIM>
  {
  public:
    /// Constructor: Call constructors for QElement and
    /// Advection Diffusion equations
    QScalarAdvectionElement()
      : QElement<DIM, NNODE_1D>(), ScalarAdvectionEquations<DIM>()
    {
    }
 
    /// Broken copy constructor
    QScalarAdvectionElement(
      const QScalarAdvectionElement<DIM, NNODE_1D>& dummy) = delete;
 
    /// Broken assignment operator
    /*void operator=(
     const QScalarAdvectionElement<DIM,NNODE_1D>&) = delete;*/
 
    ///  Required  # of `values' (pinned or dofs)
    /// at node n
    inline unsigned required_nvalue(const unsigned& n) const
    {
      return 1;
    }
 
    /// Output function:
    ///  x,y,u   or    x,y,z,u
    void output(std::ostream& outfile)
    {
      ScalarAdvectionEquations<DIM>::output(outfile);
    }
 
    /// Output function:
    ///  x,y,u   or    x,y,z,u at n_plot^DIM plot points
    void output(std::ostream& outfile, const unsigned& n_plot)
    {
      ScalarAdvectionEquations<DIM>::output(outfile, n_plot);
    }
 
 
    /*/// C-style output function:
   ///  x,y,u   or    x,y,z,u
   void output(FILE* file_pt)
    {
     ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt);
    }
 
   ///  C-style output function:
   ///   x,y,u   or    x,y,z,u at n_plot^DIM plot points
   void output(FILE* file_pt, const unsigned &n_plot)
    {
     ScalarAdvectionEquations<NFLUX,DIM>::output(file_pt,n_plot);
    }
 
   /// Output function for an exact solution:
   ///  x,y,u_exact   or    x,y,z,u_exact at n_plot^DIM plot points
   void output_fct(std::ostream &outfile, const unsigned &n_plot,
                   FiniteElement::SteadyExactSolutionFctPt
                   exact_soln_pt)
    {
     ScalarAdvectionEquations<NFLUX,DIM>::
      output_fct(outfile,n_plot,exact_soln_pt);}
 
 
   /// Output function for a time-dependent exact solution.
   ///  x,y,u_exact   or    x,y,z,u_exact at n_plot^DIM plot points
   /// (Calls the steady version)
   void output_fct(std::ostream &outfile, const unsigned &n_plot,
                   const double& time,
                   FiniteElement::UnsteadyExactSolutionFctPt
                   exact_soln_pt)
    {
     ScalarAdvectionEquations<NFLUX,DIM>::
      output_fct(outfile,n_plot,time,exact_soln_pt);
      }*/
 
 
  protected:
    /// Shape, test functions & derivs. w.r.t. to global coords. Return
    /// Jacobian.
    inline 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. at
    /// integration point ipt. Return Jacobian.
    inline double dshape_and_dtest_eulerian_at_knot_flux_transport(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsidx,
      Shape& test,
      DShape& dtestdx) const;
  };
 
  // Inline functions:
 
 
  //======================================================================
  /// Define the shape functions and test functions and derivatives
  /// w.r.t. global coordinates and return Jacobian of mapping.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double QScalarAdvectionElement<DIM, NNODE_1D>::
    dshape_and_dtest_eulerian_flux_transport(const Vector<double>& s,
                                             Shape& psi,
                                             DShape& dpsidx,
                                             Shape& test,
                                             DShape& dtestdx) const
  {
    // Call the geometrical shape functions and derivatives
    double J = this->dshape_eulerian(s, psi, dpsidx);
 
    // Loop over the test functions and derivatives and set them equal to the
    // shape functions
    for (unsigned i = 0; i < NNODE_1D; i++)
    {
      test[i] = psi[i];
      for (unsigned j = 0; j < DIM; j++)
      {
        dtestdx(i, j) = dpsidx(i, j);
      }
    }
 
    // Return the jacobian
    return J;
  }
 
 
  //======================================================================
  /// Define the shape functions and test functions and derivatives
  /// w.r.t. global coordinates and return Jacobian of mapping.
  ///
  /// Galerkin: Test functions = shape functions
  //======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  double QScalarAdvectionElement<DIM, NNODE_1D>::
    dshape_and_dtest_eulerian_at_knot_flux_transport(const unsigned& ipt,
                                                     Shape& psi,
                                                     DShape& dpsidx,
                                                     Shape& test,
                                                     DShape& dtestdx) const
  {
    // Call the geometrical shape functions and derivatives
    double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
 
    // Set the test functions equal to the shape functions (pointer copy)
    test = psi;
    dtestdx = dpsidx;
 
    // Return the jacobian
    return J;
  }
 
 
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// Face geometry for the QScalarAdvectionElement elements:
  /// The spatial dimension of the face elements is one lower than that
  /// of the bulk element but they have the same number of points along
  /// their 1D edges.
  //=======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class FaceGeometry<QScalarAdvectionElement<DIM, NNODE_1D>>
    : public virtual QElement<DIM - 1, NNODE_1D>
  {
  public:
    /// Constructor: Call the constructor for the
    /// appropriate lower-dimensional QElement
    FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
  };
 
 
  //=================================================================
  /// General DGScalarAdvectionClass. Establish the template parameters
  //===================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class DGScalarAdvectionElement
  {
  };
 
 
  //==================================================================
  // Specialization for 1D DG Advection element
  //==================================================================
  template<unsigned NNODE_1D>
  class DGScalarAdvectionElement<1, NNODE_1D>
    : public QScalarAdvectionElement<1, NNODE_1D>, public DGElement
  {
    friend class DGScalarAdvectionFaceElement<
      DGScalarAdvectionElement<1, NNODE_1D>>;
 
  public:
    // There is a single required n_flux
    unsigned required_nflux()
    {
      return 1;
    }
 
    // Calculate averages
    void calculate_element_averages(double*& average_value)
    {
      FluxTransportEquations<1>::calculate_element_averages(average_value);
    }
 
    // Constructor
    DGScalarAdvectionElement()
      : QScalarAdvectionElement<1, NNODE_1D>(), DGElement()
    {
    }
 
    ~DGScalarAdvectionElement() {}
 
    void build_all_faces()
    {
      // Make the two faces
      Face_element_pt.resize(2);
      // Make the face on the left
      Face_element_pt[0] =
        new DGScalarAdvectionFaceElement<DGScalarAdvectionElement<1, NNODE_1D>>(
          this, -1);
      // Make the face on the right
      Face_element_pt[1] =
        new DGScalarAdvectionFaceElement<DGScalarAdvectionElement<1, NNODE_1D>>(
          this, +1);
    }
 
 
    /// Compute the residuals for the Navier--Stokes equations;
    /// flag=1(or 0): do (or don't) compute the Jacobian as well.
    void fill_in_generic_residual_contribution_flux_transport(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix,
      unsigned flag)
    {
      QScalarAdvectionElement<1, NNODE_1D>::
        fill_in_generic_residual_contribution_flux_transport(
          residuals, jacobian, mass_matrix, flag);
 
      this->add_flux_contributions_to_residuals(residuals, jacobian, flag);
    }
  };
 
 
  //=======================================================================
  /// Face geometry of the 1D  DG elements
  //=======================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<DGScalarAdvectionElement<1, NNODE_1D>>
    : public virtual PointElement
  {
  public:
    FaceGeometry() : PointElement() {}
  };
 
 
  //==================================================================
  /// Specialisation for 2D DG Elements
  //==================================================================
  template<unsigned NNODE_1D>
  class DGScalarAdvectionElement<2, NNODE_1D>
    : public QScalarAdvectionElement<2, NNODE_1D>, public DGElement
  {
    friend class DGScalarAdvectionFaceElement<
      DGScalarAdvectionElement<2, NNODE_1D>>;
 
  public:
    // There is a single required n_flux
    unsigned required_nflux()
    {
      return 1;
    }
 
    // Calculate averages
    void calculate_element_averages(double*& average_value)
    {
      FluxTransportEquations<2>::calculate_element_averages(average_value);
    }
 
    // Constructor
    DGScalarAdvectionElement()
      : QScalarAdvectionElement<2, NNODE_1D>(), DGElement()
    {
    }
 
    ~DGScalarAdvectionElement() {}
 
    void build_all_faces()
    {
      Face_element_pt.resize(4);
      Face_element_pt[0] =
        new DGScalarAdvectionFaceElement<DGScalarAdvectionElement<2, NNODE_1D>>(
          this, 2);
      Face_element_pt[1] =
        new DGScalarAdvectionFaceElement<DGScalarAdvectionElement<2, NNODE_1D>>(
          this, 1);
      Face_element_pt[2] =
        new DGScalarAdvectionFaceElement<DGScalarAdvectionElement<2, NNODE_1D>>(
          this, -2);
      Face_element_pt[3] =
        new DGScalarAdvectionFaceElement<DGScalarAdvectionElement<2, NNODE_1D>>(
          this, -1);
    }
 
 
    /// Compute the residuals for the Navier--Stokes equations;
    /// flag=1(or 0): do (or don't) compute the Jacobian as well.
    void fill_in_generic_residual_contribution_flux_transport(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix,
      unsigned flag)
    {
      QScalarAdvectionElement<2, NNODE_1D>::
        fill_in_generic_residual_contribution_flux_transport(
          residuals, jacobian, mass_matrix, flag);
 
      this->add_flux_contributions_to_residuals(residuals, jacobian, flag);
    }
  };
 
 
  //=======================================================================
  /// Face geometry of the DG elements
  //=======================================================================
  template<unsigned NNODE_1D>
  class FaceGeometry<DGScalarAdvectionElement<2, NNODE_1D>>
    : public virtual QElement<1, NNODE_1D>
  {
  public:
    FaceGeometry() : QElement<1, NNODE_1D>() {}
  };
 
 
} // namespace oomph
 
#endif