dg_elements.cc

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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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//====================================================================
//Non-inline member functions for discontinuous galerkin elements

//oomph-lib includes
#include "dg_elements.h"
#include "shape.h"
#include <iomanip>

namespace oomph
{
//===================================================================
///Find pointers to neighbouring faces and the local coordinates
///in those faces that correspond to the integration points in the
///present face.
///This is achieved by moving up to the bulk element and thence
///the mesh which MUST have implemented a neighbour finding scheme
//==================================================================
void DGFaceElement::setup_neighbour_info(
 const bool &add_neighbour_data_to_bulk)
{
 //Cache the pointer to the bulk element
 DGElement* const bulk_element_pt = 
  dynamic_cast<DGElement*>(this->bulk_element_pt());
 
 //Find the number of points in the integration scheme
 const unsigned n_intpt = integral_pt()->nweight();
 //Resize the storage in the element
 Neighbour_face_pt.resize(n_intpt);
 Neighbour_local_coordinate.resize(n_intpt);

 //If we are adding the neighbour data to the bulk element
 //then resize this storage
 if(add_neighbour_data_to_bulk) {Neighbour_external_data.resize(n_intpt);}

 
 //Get the dimension of the present element
 const unsigned el_dim = this->dim();
 //Local coordinate in the face element
 Vector<double> s(el_dim);
 
 //Get the dimension of the bulk element
 const unsigned n_dim = bulk_element_pt->dim(); 
 //Local coordinate in the bulk element
 Vector<double> s_bulk(n_dim);

 //Storage for the interpolation data in the neighbour
 Vector<Data*> neighbour_data;
 
 //Now loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Get the local coordinate of the integration points
   for(unsigned i=0;i<el_dim;i++)
    {s[i] = integral_pt()->knot(ipt,i);}
   
   //Now get the bulk coordinate
   this->get_local_coordinate_in_bulk(s,s_bulk);

   //Now we find the neighbouring face via the bulk element's member
   //function which calls the Mesh's member function
   bulk_element_pt->
    get_neighbouring_face_and_local_coordinate(this->face_index(),
     s_bulk,Neighbour_face_pt[ipt],Neighbour_local_coordinate[ipt]);

   //If we are adding the external data to the bulk
   if(add_neighbour_data_to_bulk)
    {
     //Throw the appropriate data from the neighbour face to the set
     dynamic_cast<DGFaceElement*>(Neighbour_face_pt[ipt])
      ->get_interpolation_data(neighbour_data);
     
     //Find the number of data
     unsigned n_neighbour_data = neighbour_data.size();
     //Resize the storage accordingly
     Neighbour_external_data.resize(n_neighbour_data);
     
     //Add the data to the external data of the bulk element
     for(unsigned n=0;n<n_neighbour_data;n++)
      {
       Neighbour_external_data[ipt][n]
        = bulk_element_pt->add_external_data(neighbour_data[n]);
      }
    }
  }
}


//======================================================================
//Report the global coordinates corresponding to the integration points
//and the corresponding coordinates in the neighbouring faces
//======================================================================
void DGFaceElement::report_info()
{
 //Find the number of nodes
 const unsigned n_node = this->nnode();
 //Allocate storage for the shape functions
 Shape psi(n_node);
 
 //Find the dimension of the problem
 const unsigned dim = this->nodal_dimension();
 //Storage for local coordinates in this face and its neighbour
 Vector<double> x(dim), face_x(dim);
 
 //Find the dimension of the element
 const unsigned el_dim = this->dim();
 //Storage for local coordinate
 Vector<double> s(el_dim);
 
 //Calculate the number of integration points from the array
 const unsigned n_intpt = this->integral_pt()->nweight();
 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Find the local coordinate at the knot point
   for(unsigned i=0;i<el_dim;i++) 
    {s[i] = this->integral_pt()->knot(ipt,i);}
   //Get the value of the global coordinate in the present face
   this->interpolated_x(s,x);
   
   //Get the value of x in the neighbouring face
   Neighbour_face_pt[ipt]->interpolated_x(Neighbour_local_coordinate[ipt],
                                          face_x);
   
   //Let's compare
   oomph_info <<  "In Face                   In Neighbour\n";
   for(unsigned i=0;i<dim;i++)
    {
     if(i==0) {oomph_info << "(";}
     else {oomph_info << ", ";}
     oomph_info << std::setw(5) << std::left <<  x[i];
    }
   oomph_info << ")";
   
   oomph_info << "                   ";
   
   for(unsigned i=0;i<dim;i++)
    {
     if(i==0) {oomph_info << "(";}
     else {oomph_info << ", ";}
     oomph_info << std::setw(5) << std::left << face_x[i];
    }
   oomph_info << ")";
   oomph_info << std::endl;
  }
}


//=====================================================================
///Return the interpolated values of the unknown fluxes
//=====================================================================
void DGFaceElement::interpolated_u(const Vector<double> &s, 
                                   Vector<double> &u)
{
 //Find the number of nodes
 const unsigned n_node = nnode();
 //If there are no nodes then return immediately
 if(n_node==0) {return;}
 
 //Get the shape functions at the local coordinate
 Shape psi(n_node);
 this->shape(s,psi);
 
 //Find the number of fluxes
 const unsigned n_flux = this->required_nflux();
 
 //Find the indices at which the local fluxes are stored
 Vector<unsigned> flux_nodal_index(n_flux);
 for(unsigned i=0;i<n_flux;i++)
  {
   flux_nodal_index[i] = this->flux_index(i);
  }
 //Initialise the fluxes to zero
 for(unsigned i=0;i<n_flux;i++) {u[i] = 0.0;}
 
 //Loop over the nodes
 for(unsigned n=0;n<n_node;n++)
  {
   const double psi_ = psi[n];
   for(unsigned i=0;i<n_flux;i++)
    {
     u[i] += this->node_pt(n)->value(flux_nodal_index[i])*psi_;
    }
  }
}

//=====================================================================
///Add all the data that are used to interpolate the unknowns. This must
///be consistent with the interpolated_u function above and the default
///implementation assumes pure nodal interpolaton.
//=====================================================================
void DGFaceElement::get_interpolation_data(
 Vector<Data*> &interpolation_data)
{
 //Find the number of nodes
 const unsigned n_node = nnode();
 //Now resize the vector
 interpolation_data.resize(n_node);

 //If there are no nodes then return immediately
 if(n_node==0) {return;}

 //Otherwise loop over the nodes and add to the vector in order
 for(unsigned n=0;n<n_node;n++)
  {
   interpolation_data[n] = this->node_pt(n);
  }
}

//====================================================================
///Calculate the numerical flux at the knot point ipt. This is the 
///most general interface than can be overloaded if desired. The shape
///functions at the knot point will be passed into this function.
//====================================================================
void DGFaceElement::numerical_flux_at_knot(const unsigned &ipt,
                                           const Shape &psi,
                                           Vector<double> &flux,
                                           DenseMatrix<double> &dflux_du_int,
                                           DenseMatrix<double> &dflux_du_ext,
                                           unsigned flag)
{
 //Find the number of nodes
 const unsigned n_node = this->nnode();
 //Find the nodal dimension
 const unsigned nodal_dim = this->nodal_dimension();
 //Number of fluxes
 const unsigned n_flux = this->required_nflux();
 //Find the indices at which the local fluxes are stored
 Vector<unsigned> flux_nodal_index(n_flux);
 for(unsigned i=0;i<n_flux;i++)
  {
   flux_nodal_index[i] = this->flux_index(i);
  }
 
 //Calculate the local unknowns
 Vector<double> interpolated_u(n_flux,0.0);
 
 //Loop over the shape functions
 for(unsigned l=0;l<n_node;l++)
  {
   //Cache the shape functions
   const double psi_ = psi(l);
   //Loop over the fluxes
   for(unsigned i=0;i<n_flux;i++)
    {
     //Calculate the velocity from the most recent nodal values
     interpolated_u[i] += this->nodal_value(l,flux_nodal_index[i])*psi_;
    }
  }
 
 //Now calculate the outer unit normal Vector
 Vector<double> interpolated_n(nodal_dim);
 this->outer_unit_normal(ipt,interpolated_n);
 
 //Get the pointer to the neighbour
 DGFaceElement* neighbour_element_pt =
  dynamic_cast<DGFaceElement*>(Neighbour_face_pt[ipt]);
 
 //Get the neighbour's fluxes 
 Vector<double> interpolated_u_neigh(n_flux); 
 
 neighbour_element_pt->
  interpolated_u(Neighbour_local_coordinate[ipt],
                 interpolated_u_neigh);

 //Call the "standard" numerical flux function
 this->numerical_flux(interpolated_n,interpolated_u,
                      interpolated_u_neigh,flux);
 
 //If we are calculating the jacobian add this term
 if(flag && (flag < 3))
  {
   this->dnumerical_flux_du(interpolated_n,interpolated_u,
                            interpolated_u_neigh,dflux_du_int,
                            dflux_du_ext);
  }

}

//===============================================================
///\short Calculate the derivative of the 
///normal flux, which is the dot product of our
///approximation to the flux with the outer unit normal,
///with respect to the interior and exterior variables
/// Default is to use finite differences
//=============================================================
void DGFaceElement::dnumerical_flux_du(const Vector<double> &n_out,
                                       const Vector<double> &u_int,
                                       const Vector<double> &u_ext,
                                       DenseMatrix<double> &dflux_du_int,
                                       DenseMatrix<double> &dflux_du_ext)
{
 //Find the number of fluxes
 const unsigned n_flux = this->required_nflux();

 //Get a local copy of the unknowns
 Vector<double> u_int_local = u_int;
 Vector<double> u_ext_local = u_ext;
 
 //Storage for incremented and decremented fluxes
 Vector<double> flux_plus(n_flux), flux_minus(n_flux);
 
 const double fd_step = GeneralisedElement::Default_fd_jacobian_step;

 //Now loop over all the fluxes
 for(unsigned n=0;n<n_flux;n++)
  {
   //Increase internal value

   //Store the old value
   double old_var = u_int_local[n];
   //Increment the value
   u_int_local[n] += fd_step;
   //Get the new values
   this->numerical_flux(n_out,u_int_local,u_ext_local,flux_plus);
   
   //Reset the value
   u_int_local[n] = old_var;
   //Decrement the value
   u_int_local[n] -= fd_step;
   //Get the new values
   this->numerical_flux(n_out,u_int_local,u_ext_local,flux_minus);
         
   //Assemble the column of the jacobian
   for(unsigned m=0;m<n_flux;m++)
    {
     dflux_du_int(m,n) = (flux_plus[m] - flux_minus[m])/(2.0*fd_step);
    }
   
   //Reset the value
   u_int_local[n] = old_var;
   
   //Increase external value

   //Store the old value
   old_var = u_ext_local[n];
   //Increment the value
   u_ext_local[n] += fd_step;
   //Get the new values
   this->numerical_flux(n_out,u_int_local,u_ext_local,flux_plus);
   
   //Reset the value
   u_ext_local[n] = old_var;
   //Decrement the value
   u_ext_local[n] -= fd_step;
   //Get the new values
   this->numerical_flux(n_out,u_int_local,u_ext_local,flux_minus);
         
   //Assemble the column of the jacobian
   for(unsigned m=0;m<n_flux;m++)
    {
     dflux_du_ext(m,n) = (flux_plus[m] - flux_minus[m])/(2.0*fd_step);
    }
   
   //Reset the value
   u_ext_local[n] = old_var;
 }
}


//===================================================================
///Calculate the integrated (numerical) flux out of the face and add
///it to the residuals vector
//===================================================================
void DGFaceElement::add_flux_contributions(Vector<double> &residuals,
                                           DenseMatrix<double> &jacobian,
                                           unsigned flag)
{
//Check that we have set up the coupling data if we are computing the jacobin
#ifdef PARANOID
 if(flag && (flag < 3))
  {
   if(Neighbour_external_data.size()==0)
    {
     std::ostringstream error_stream;
     error_stream <<   
      "Coupling data between elements not included in jacobian\n"
                  << 
      "You should call DGMesh::setup_face_neighbour_info(true) to ensure\n"
                  << 
      "that this information is included in the jacobian\n";
      
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
  }
#endif


 //Find the number of nodes
 const unsigned n_node = nnode();
 //Dimension of the face
 const unsigned el_dim = dim();
 //Storage for the shape functions
 Shape psi(n_node);

 //Number of integration points
 const unsigned n_intpt = this->integral_pt()->nweight();
 //Number of fluxes
 const unsigned n_flux = this->required_nflux();
 //Find the indices at which the local fluxes are stored
 Vector<unsigned> flux_nodal_index(n_flux);
 for(unsigned i=0;i<n_flux;i++)
  {
   flux_nodal_index[i] = this->flux_index(i);
  }

 //Cache the bulk element
 DGElement* bulk_elem_pt =  dynamic_cast<DGElement*>(this->bulk_element_pt());
 
 //Storage for the flux and its derivatives 
 Vector<double> F(n_flux);
 DenseMatrix<double> dF_du_int(n_flux,n_flux);
 DenseMatrix<double> dF_du_ext(n_flux,n_flux);

 //Loop over the integration points
 for(unsigned ipt=0;ipt<n_intpt;ipt++)
  {
   //Get the integral weight
   double W = this->integral_pt()->weight(ipt);
   //Get the shape functions at the knot
   this->shape_at_knot(ipt,psi);
     
   //Calculate the Jacobian 
   //For a point element, it's one
   double J=W;
   //Otherwise calculate for the element
   if(el_dim != 0) {J *= this->J_eulerian_at_knot(ipt);}
   
   //Now calculate the numerical flux (and derivatives)
   this->numerical_flux_at_knot(ipt,psi,F,dF_du_int,dF_du_ext,flag);
   
   //Limit if desired here
 
   //Cache the pointer to the neighbour
   DGFaceElement* neighbour_element_pt =
    dynamic_cast<DGFaceElement*>(Neighbour_face_pt[ipt]);
   
   //Finally  we need to assemble the appropriate contributions
   //to the residuals
   for(unsigned l=0;l<n_node;l++)
    {
     //Loop over the fluxes
     for(unsigned i=0;i<n_flux;i++)
      {
       //Get the local equation number in the bulk
       int local_eqn = 
        bulk_elem_pt->nodal_local_eqn(bulk_node_number(l),flux_nodal_index[i]);
       
       //If it's not a boundary condition
       if(local_eqn >= 0)
        {
         //Add the flux multiplied by the shape function and the jacobian
         residuals[local_eqn] -= psi(l)*F[i]*J;

         //Add the jacobian contributions
         if(flag)
          {
           //If we are assembling the jacobian
           if(flag < 3)
            {
             
             //Loop over the internal nodes and fluxes again
             for(unsigned l2=0;l2<n_node;l2++)
              {
               for(unsigned i2=0;i2<n_flux;i2++)
                {
                 //Get the local unknown equation number in the bulk
                 int local_unknown = bulk_elem_pt
                  ->nodal_local_eqn(bulk_node_number(l2),flux_nodal_index[i2]);
                 
                 //If it's not a boundary condition
                 if(local_unknown >= 0)
                  {
                   //Add the flux multiplied by the shape function a
                   //nd the jacobian
                   jacobian(local_eqn,local_unknown) -= 
                    psi(l)*dF_du_int(i,i2)*psi(l2)*J;
                  }
                }
              }
             
             //How many nodes does it have
             unsigned neigh_n_node = neighbour_element_pt->nnode();
             
             //Get the flux indices from the neighbour
             Vector<unsigned> neigh_flux_index(n_flux);
             for(unsigned i2=0;i2<n_flux;i2++) 
              {
               neigh_flux_index[i2] = neighbour_element_pt->flux_index(i2);
              }
             
             //Loop over the neighbours nodes
             for(unsigned l2=0;l2<neigh_n_node;l2++)
              {
               //Loop over the fluxed
               for(unsigned i2=0;i2<n_flux;i2++)
                {
                 //Get the local unknown equation number in the bulk
                 int local_unknown = 
                  dynamic_cast<DGElement*>(this->bulk_element_pt())
                  ->external_local_eqn(
                   Neighbour_external_data[ipt][l2],flux_nodal_index[i2]);
                 
                 //If it's not a boundary condition
                 if(local_unknown >= 0)
                  {
                   //Add the flux multiplied by the shape function 
                   //and the jacobian
                   jacobian(local_eqn,local_unknown) -= 
                    psi(l)*dF_du_ext(i,i2)*psi(l2)*J;
                  }
                }
              }
            }
          } //End of jacobian calculation
         
        }
      }
    }
  }
}

//========================================================================
///\short Function that computes and stores the (inverse) mass matrix
//========================================================================
void DGElement::pre_compute_mass_matrix()
{
 //Now let's assemble stuff
 const unsigned n_dof = this->ndof();
 //Allocate storage for the local mass matrix (if required)
 if(M_pt==0) {M_pt = new DenseDoubleMatrix;}
 
 //Resize and initialise the vector that will holds the residuals
 Vector<double> dummy(n_dof,0.0);
 
 //Resize the mass matrix
 M_pt->resize(n_dof,n_dof);
 //Initialise the entries to zero
 M_pt->initialise(0.0);
 //Get the local mass matrix and residuals
 this->fill_in_contribution_to_mass_matrix(dummy,*M_pt);
 
 //Now invert the mass matrix it will always be small
 //This can possibly be streamlined (for example in spectral
 //elements the mass matrix is diagonal)
 M_pt->ludecompose();
 
 //The mass matrix has been computed
 Mass_matrix_has_been_computed=true;
}



//============================================================================
///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 DGElement::
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();
 //Allocate storage for the local mass matrix (if required)
 if(M_pt==0) {M_pt = new DenseDoubleMatrix;}
 
 //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
 {
  //Resize the mass matrix
  M_pt->resize(n_dof,n_dof);
  //Initialise the entries to zero
  M_pt->initialise(0.0);
  //Get the local mass matrix and residuals
  this->fill_in_contribution_to_mass_matrix(minv_res,*M_pt);
 
  //Now invert the mass matrix it will always be small
  //This can possibly be streamlined (for example in spectral
  //elements the mass matrix is diagonal)
  M_pt->ludecompose();

  //The mass matrix has been computed
  Mass_matrix_has_been_computed=true;
 }
 
 //Always do the backsubstitution
 M_pt->lubksub(minv_res);
}


void DGElement::get_neighbouring_face_and_local_coordinate(
 const int &face_index,
 const Vector<double> &s, FaceElement* &face_element_pt,
 Vector<double> &s_face)
{
 DG_mesh_pt->neighbour_finder(this,face_index,s,face_element_pt,s_face);
}


///Limit the slope within an element
void DGElement::slope_limit(SlopeLimiter* const &slope_limiter_pt)
{
 //Firstly find the dimension
 const unsigned n_dim = this->dim();
 //Find the number of fluxes
 const unsigned n_flux = this->required_nflux();

 switch(n_dim)
  {
   //One dimensional (easy-ish) case
  case 1:
  {
   //Storage for the element and its neighbours
   Vector<DGElement*> required_element_pt(3);
   required_element_pt[0] = this;
   
   //Get the pointer to the element on the left
   required_element_pt[1] = dynamic_cast<DGElement*>(
    dynamic_cast<DGFaceElement*>(this->face_element_pt(0))
    ->neighbour_face_pt(0)->bulk_element_pt());
   //Get the pointer to the element on the right
   required_element_pt[2] = dynamic_cast<DGElement*>(
    dynamic_cast<DGFaceElement*>(this->face_element_pt(1))
    ->neighbour_face_pt(0)->bulk_element_pt());

   //Loop over the fluxed
   for(unsigned i=0;i<n_flux;i++)
    {
     //Call our limiter, which will take as it's arguments, the current
     //element and the required neighbours
     slope_limiter_pt->limit(i,required_element_pt);
    }
  }
  break;
  
  default:
  {
   std::ostringstream error_stream;
   error_stream <<"Slope limiting is not implemented for this dimension: "
                << n_dim;
   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
  }
}

double DGMesh::FaceTolerance = 1.0e-10;


//====================================================
///Helper minmod function
//====================================================
double MinModLimiter::minmod(Vector<double> &args)
{
 const unsigned n_arg = args.size();
 //If no arguments return zero
 if(n_arg==0) {return 0.0;}

 //Initialise the sign from the sign of the first entry
 int sign = 0;
 if(args[0] < 0.0) {sign = -1;}
 else if(args[0] > 0.0) {sign = 1;}
 else {return 0.0;}
 
 //Initialise the minimum value
 double min = std::fabs(args[0]);
 
 //Now loop over the rest of the values
 for(unsigned i=1;i<n_arg;i++)
  {
   if(sign == 1)
    {
     if(args[i] < 0.0) {return 0.0;}
       else if(args[i] < min) {min = args[i];}
    }
   else if(sign == -1)
    {
     if(args[i] > 0.0) {return 0.0;}
     else if(std::fabs(args[i]) < min) {min = std::fabs(args[i]);}
    }
  }
 
 //Make sure to return the sign multiplied by the minimum value
 return sign*min;
}


///Modified minmod limiter to fix behaviour in smooth regions
double MinModLimiter::minmodB(Vector<double> &args, const double &h)
{
 const unsigned n_arg = args.size();
 //If no arguments return zero
 if(n_arg==0) {return 0.0;}
  
 //Modification to fix extrema
 if(std::fabs(args[0]) < this->M*h*h) {return args[0];}
 //Otherwise just return the usual minmod
 else {return minmod(args);}
}

///Implement the limiter function for the basic MinModlimiter
void MinModLimiter::limit(const unsigned &i, 
                          const Vector<DGElement*> &required_element_pt)
{
 //Set the tolerance
 const double tol = 1.0e-16;

 //Find the geometric parameters of the element
 const unsigned n_node = required_element_pt[0]->nnode();
 const double x_l = required_element_pt[0]->node_pt(0)->x(0);
 const double x_r = required_element_pt[0]->node_pt(n_node-1)->x(0); 
 const double h = x_r - x_l;
 const double x0 = 0.5*(x_l + x_r);                                           
 //Find the average value
 const double u_av = required_element_pt[0]->average_value(i);

 //Storage for the gradients to the minmod function
 Vector<double> arg; arg.reserve(3);
 
 //Add the approximation for the element's left flux
 arg.push_back(u_av - required_element_pt[0]->node_pt(0)->value(i));
 
 //If there is a left element add the approximate gradient 
 //to the argument vector
 if(required_element_pt[1] != required_element_pt[0]) 
  {arg.push_back(u_av - required_element_pt[1]->average_value(i));}
 //If there is a right element add the approximate gradient 
 //to the argument vector
 if(required_element_pt[2] != required_element_pt[0]) 
  {arg.push_back(required_element_pt[2]->average_value(i) - u_av);}
 
 //Set the left value
 const double u_l = u_av - this->minmod(arg);

 //Now replace the first term in the argument list with the 
 //approximation for the element's right flux
 arg.front() = required_element_pt[0]->node_pt(n_node-1)->value(i) - u_av;
 
 //Set the right value
 const double u_r = u_av + this->minmod(arg);

 //If the basic limited values are different from 
 //the unlimited values then limit
 if((std::fabs(u_l - required_element_pt[0]->node_pt(0)->value(i)) > tol) &&
    (std::fabs(u_r - required_element_pt[0]->node_pt(n_node-1)->value(i)) 
     > tol))
  {
   //Find the centre of the element on the left
   const double x0_l = 
    0.5*(required_element_pt[1]->node_pt(required_element_pt[1]->nnode()-1)
         ->x(0) + required_element_pt[1]->node_pt(0)->x(0));
   //Find the centre of the element on the right
   const double x0_r = 
    0.5*(required_element_pt[2]->node_pt(required_element_pt[2]->nnode()-1)
         ->x(0) + required_element_pt[2]->node_pt(0)->x(0));
   
   //Clear the argument list and reserve its size to 
   arg.clear(); arg.reserve(3);
   
   //Approximate the gradient over the whole cell
   arg.push_back((required_element_pt[0]->node_pt(n_node-1)->value(i) - 
                  required_element_pt[0]->node_pt(0)->value(i))/h);
   
   //Adjust the estimates for the gradient calculated from neighbouring
   //averages for the straight min-mod and MUSCL limiters
   double gradient_factor = 0.5;
   if(MUSCL) {gradient_factor = 1.0;}
   
   //If there is a left element, form the gradient
   if(required_element_pt[0]!=required_element_pt[1])
    {
     arg.push_back(
      (u_av - required_element_pt[1]->average_value(i))/
      (gradient_factor*(x0 - x0_l)));
    }

   //If there is a right element, form the gradient
   if(required_element_pt[0]!=required_element_pt[2])
    {
     arg.push_back(
      (required_element_pt[2]->average_value(i) - u_av)/
      (gradient_factor*(x0_r - x0)));
    }
   
   //Calculate the limited gradient of these three
   double limited_gradient = this->minmodB(arg,h);
   
   //Loop over the nodes and limit
   for(unsigned n=0;n<n_node;n++)
    {
     double x = required_element_pt[0]->node_pt(n)->x(0) - x0;
     required_element_pt[0]->node_pt(n)
      ->set_value(i,u_av  + x*limited_gradient);
    }
  }
}

} //End of namespace