dg_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//    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//====================================================================
//Header functions for classes that define Discontinuous Galerkin elements

//Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_DG_ELEMENT_HEADER
#define OOMPH_DG_ELEMENT_HEADER

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

//Oomph-lib header
#include "elements.h"
#include "mesh.h"

namespace oomph
{
//=============================================================
/// \short Base class for Discontinuous Galerkin Faces.
/// These are responsible for calculating the normal fluxes
/// that provide the communication between the discontinuous 
/// elements.
//===============================================================
class DGFaceElement : public virtual FaceElement
 {
  /// Vector of neighbouring face elements at the integration points
  Vector<FaceElement*> Neighbour_face_pt;
  
  /// Vector of neighbouring local coordinates at the integration points
  Vector<Vector<double> > Neighbour_local_coordinate;

  /// Vector of the vectors that will store the number of the 
  /// bulk external data that correspond to the dofs in the neighbouring face
  /// This is only used if we are using implict timestepping and wish
  /// to assemble the jacobian matrix. In order that the data be correctly
  /// set up DGMesh::setup_face_neighbour_info() must be called with the
  /// boolean flag set to true.
  Vector<Vector<unsigned> > Neighbour_external_data;

 protected:

  /// \short Return the index at which the i-th unknown flux is stored. 
  // The default return is suitable for single-physics problem
  virtual inline unsigned flux_index(const unsigned &i) const {return i;}

  /// Set the number of flux components
  virtual unsigned required_nflux() {return 0;}

 public:

  /// Empty Constructor
  DGFaceElement() : FaceElement()  { }
 
  /// Empty Destructor 
  virtual ~DGFaceElement() {}

  /// Access function for neighbouring face information
  FaceElement* neighbour_face_pt(const unsigned &i)
   {return Neighbour_face_pt[i];}
  
  /// Setup information from the neighbouring face, return a set
  /// of nodes (as data) in the neighour. The boolean flag
  /// is used to determine whether the data in the neighbour are
  /// added as external data to the bulk element --- required when
  /// computing the jacobian of the system
  void setup_neighbour_info(const bool &add_neighbour_data_to_bulk);

  ///Output information about the present element and its neighbour
  void report_info();

  //Get the value of the unknowns
  virtual void interpolated_u(const Vector<double> &s, Vector<double> &f);
 
  ///\short Get the data that are used to interpolate the unkowns
  ///in the element. These must be returned in order.
  virtual void get_interpolation_data(Vector<Data*> &interpolation_data);

  
  ///\short Calculate the normal numerical flux at the integration point.
  ///This is the most general interface that can be overloaded if desired
  virtual void 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);

  ///\short Calculate the normal flux, which is the dot product of our
  ///approximation to the flux with the outer unit normal
  virtual void numerical_flux(const Vector<double> &n_out,
                              const Vector<double> &u_int,
                              const Vector<double> &u_ext,
                              Vector<double> &flux)
   {
    std::ostringstream error_stream;
    error_stream 
     << "Empty numerical flux function called\n"
     << "This function should be overloaded with a specific flux\n"
     << "that is appropriate to the equations being solved.\n";
    throw OomphLibError(error_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }

  ///\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
  virtual void 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);

  
  ///\short Add the contribution from integrating the numerical flux
  //over the face to the residuals
  void add_flux_contributions(Vector<double> &residuals,
                              DenseMatrix<double> &jacobian,
                              unsigned flag);

 };

class DGMesh;
class SlopeLimiter;

//==================================================================
/// A Base class for DGElements
//=================================================================
class DGElement : public virtual FiniteElement
{
protected:

 ///\short The DGFaceElement requires access to the nodal local equation
 ///information, so it's a friend
 friend class DGFaceElement;
 
 ///\short Vector of pointers to faces of the element
 Vector<FaceElement*> Face_element_pt;
 
 ///Pointer to Mesh, which will be responsible for the neighbour finding
 DGMesh *DG_mesh_pt;

 ///\short Pointer to storage for a mass matrix that can be recycled if
 ///desired
 DenseDoubleMatrix *M_pt;

 /// \short Pointer to storage for the average values of the of the 
 /// variables over the element
 double *Average_value;

 ///\short Boolean flag to indicate whether to reuse the mass matrix
 bool Mass_matrix_reuse_is_enabled;

 ///\short Boolean flag to indicate whether the mass matrix has been
 ///computed
 bool Mass_matrix_has_been_computed;

 ///\short Boolean flag to indicate whether the mass matrix can be
 ///deleted (i.e. was it created by this element)
 bool Can_delete_mass_matrix;

 /// Set the number of flux components
 virtual unsigned required_nflux() {return 0;}

  public:

 ///Constructor, initialise the pointers to zero
 DGElement() : DG_mesh_pt(0), M_pt(0), Average_value(0),
  Mass_matrix_reuse_is_enabled(false),
  Mass_matrix_has_been_computed(false),
  Can_delete_mass_matrix(true) { }

 ///\short Virtual destructor, destroy the mass matrix, if we created it
 ///Clean-up storage associated with average values
 virtual ~DGElement()
  {
   if((M_pt!=0) && Can_delete_mass_matrix) {delete M_pt;}
   if(this->Average_value!=0) {delete[] Average_value; Average_value=0;}
  }
 
 ///\short Access function for the boolean to indicate whether the
 ///mass matrix has been computed
 bool mass_matrix_has_been_computed() {return Mass_matrix_has_been_computed;}

 ///Function that allows the reuse of the mass matrix
 void enable_mass_matrix_reuse() 
  {
   Mass_matrix_reuse_is_enabled=true;
   Mass_matrix_has_been_computed=false;
  }

 ///Function that disables the reuse of the mass matrix
 void disable_mass_matrix_reuse()
  {
   //If we are using another element's mass matrix then reset our pointer
   //to zero
   if(!Can_delete_mass_matrix) {M_pt = 0;}
   //Otherwise we do not reuse the mass matrix
   Mass_matrix_reuse_is_enabled = false;
   //Recalculate the mass matrix
   Mass_matrix_has_been_computed=false;
  }

                               
 ///Set the mass matrix to point to one in another element
 virtual void set_mass_matrix_from_element(DGElement* const &element_pt)
  {
   //If the element's mass matrix has not been computed, compute it!
   if(!element_pt->mass_matrix_has_been_computed())
    {
     element_pt->pre_compute_mass_matrix();
    }
   
   //Now set the mass matrix in this element to address that
   //of element_pt
   this->M_pt = element_pt->M_pt;
   //We must reuse the mass matrix, or there will be trouble
   //Because we will recalculate it in the original element
   Mass_matrix_reuse_is_enabled=true;
   Mass_matrix_has_been_computed=true;
   //We cannot delete the mass matrix
   Can_delete_mass_matrix = false;
  }

 ///\short Function that computes and stores the (inverse) mass matrix
 void pre_compute_mass_matrix();

 //Function that is used to construct all the faces of the DGElement
 virtual void build_all_faces()=0;

 ///\short 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 saving tricks can be applied)
 virtual void get_inverse_mass_matrix_times_residuals(
  Vector<double> &minv_res);

 ///\short Construct all nodes and faces of the element. 
 ///The vector of booleans boundary should be the same size
 ///as the number of nodes and if any entries are true
 ///that node will be constructed as a boundary node.
 void construct_boundary_nodes_and_faces(DGMesh* const &mesh_pt,
                                         std::vector<bool> &boundary_flag,
                                         TimeStepper* const &time_stepper_pt)
  {
   //Construct the nodes (This should not be used in a base class)
   const unsigned n_node = this->nnode();
#ifdef PARANOID
   if(boundary_flag.size() != n_node)
    {
     std::ostringstream error_stream;
     error_stream 
      << "Size of boundary_flag vector is " 
      << boundary_flag.size() << "\n"
      << "It must be the same as the number of nodes in the element "
      << n_node << "\n";
     
     throw OomphLibError(error_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   //Loop over the nodes
   for(unsigned n=0;n<n_node;n++)
    {
     //If the node is marked on a boundary construct a boundary node
     if(boundary_flag[n])
      {
       (void)this->construct_boundary_node(n,time_stepper_pt);
      }
     //Otherwise, construct a normal node
     else
      {
       (void)this->construct_node(n,time_stepper_pt);
      }
    }

   //Make the faces
   this->build_all_faces();
   
   //Set the Mesh pointer
   DG_mesh_pt = mesh_pt;
  }
   

 ///Construct the nodes and faces of the element
 void construct_nodes_and_faces(DGMesh* const &mesh_pt,
                                TimeStepper* const &time_stepper_pt)
  {
   //Loop over the nodes
   const unsigned n_node = this->nnode();
   for(unsigned n=0;n<n_node;n++)
    {
     //Construct the node and ignore the return value
     (void) this->construct_node(n,time_stepper_pt);
    }
  
   //Make the faces
   this->build_all_faces();

   //Set the Mesh pointer
   DG_mesh_pt = mesh_pt;
  }

 //Set the mesh pointer of the element
 void set_mesh_pt(DGMesh* &mesh_pt) {DG_mesh_pt = mesh_pt;}

 /// Return the number of faces
 unsigned nface() const {return Face_element_pt.size();}

 ///Access function for the faces
 DGFaceElement* face_element_pt(const unsigned &i) 
  {return dynamic_cast<DGFaceElement*>(Face_element_pt[i]);}

 ///Output the faces of the element
 void output_faces(std::ostream &outfile)
  {
   //Loop over the faces
   unsigned n_face = nface();
   for(unsigned f=0;f<n_face;f++)
    {
     face_element_pt(f)->output(outfile);
    }
  }
 
 ///Return the neighbour info
 void get_neighbouring_face_and_local_coordinate(const int &face_index,
                                                 const Vector<double> &s, 
                                                 FaceElement* &face_element_pt,
                                                 Vector<double> &s_face);

 //Setup the face information
 ///The boolean flag determines whether the data from the neighbouring
 ///elements is added as external data to the element (required for
 ///correct computation of the jacobian)
 void setup_face_neighbour_info(const bool &add_face_data_as_external)
  {
   unsigned n_face = this->nface();
   for(unsigned f=0;f<n_face;f++)
    {
     face_element_pt(f)->setup_neighbour_info(add_face_data_as_external);
     face_element_pt(f)->report_info();
    }
  }


///Loop over all faces and add their integrated numerical fluxes
///to the residuals
void add_flux_contributions_to_residuals(Vector<double> &residuals,
                                         DenseMatrix<double> &jacobian,
                                         unsigned flag)
{
 //Add up the contributions from each face
 unsigned n_face = this->nface();
 for(unsigned f=0;f<n_face;f++)
  {
   face_element_pt(f)->add_flux_contributions(residuals,jacobian,flag);
  }
}

///Limit the slope within the element
 void slope_limit(SlopeLimiter* const &slope_limiter_pt); 
 
 ///Calculate the averages in the element
 virtual void calculate_element_averages(double* &average_values) 
  {
   throw OomphLibError("Default (empty) version called",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
 
 ///Calculate the elemental averages
 void calculate_averages()
  {this->calculate_element_averages(this->Average_value);}

  /// \short Return the average values
 double &average_value(const unsigned &i)
  {
   if(Average_value==0)
    {
     throw OomphLibError("Averages not calculated yet",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   return Average_value[i];
  }


 /// \short Return the average values
 const double &average_value(const unsigned &i) const
  {
   if(Average_value==0)
    {
     throw OomphLibError("Averages not calculated yet",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   return Average_value[i];
  }

};

class DGMesh : public Mesh
{


public:
 
 static double FaceTolerance;

 DGMesh() : Mesh() { }

 virtual ~DGMesh() {}

 virtual void neighbour_finder(FiniteElement* const &bulk_element_pt,
                               const int &face_index,
                               const Vector<double> &s_bulk,
                               FaceElement* &face_element_pt,
                               Vector<double> &s_face)
  {
  std::string error_message =
    "Empty neighbour_finder() has been called.\n";
   error_message +=
    "This function is implemented in the base class of a DGMesh.\n";
   error_message += 
    "It must be overloaded in a specific DGMesh\n";
   
   throw OomphLibError(error_message,
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  } 


 //Setup the face information for all the elements
 ///If the boolean flag is set to true then the data from neighbouring
 ///faces will be added as external data to the bulk element. 
 void setup_face_neighbour_info(const bool &add_face_data_as_external=false)
  {
   //Loop over all the elements and setup their face neighbour information
   const unsigned n_element = this->nelement();
   for(unsigned e=0;e<n_element;e++)
    {
     dynamic_cast<DGElement*>(this->element_pt(e))
      ->setup_face_neighbour_info(add_face_data_as_external);
    }
  }

 //Limit the slopes on the entire mesh
 void limit_slopes(SlopeLimiter* const &slope_limiter_pt)
  {
   //Loop over all the elements and calculate the averages
   const unsigned n_element = this->nelement();
   for(unsigned e=0;e<n_element;e++)
    {
     dynamic_cast<DGElement*>(this->element_pt(e))
      ->calculate_averages();
    }

   //Now loop over again and limit the values
   for(unsigned e=0;e<n_element;e++)
    {
     dynamic_cast<DGElement*>(this->element_pt(e))
      ->slope_limit(slope_limiter_pt);
    }
  }

};


//======================================================
/// \short Base class for slope limiters
//=====================================================
class SlopeLimiter
{
  public:
   
 ///Empty constructor
 SlopeLimiter() {}

 ///virtual destructor
 virtual ~SlopeLimiter() {}

 ///Basic function
 virtual void limit(const unsigned &i,
                    const Vector<DGElement*> &required_element_pt)
  {
   throw OomphLibError("Calling default empty limiter\n",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
};

class MinModLimiter : public SlopeLimiter
{
 ///\short Constant that is used in the modified min-mod function to
 ///provide better properties at extrema
 double M;

 ///Boolean flag to indicate a MUSCL or straight min-mod limiter
 bool MUSCL;

  public:
 
 ///\short Constructor takes a value for the modification parameter M 
 ///(default to zero --- classic min mod) and a flag to indicate whether
 ///we use MUSCL limiting or not --- default false
 MinModLimiter(const double &m=0.0, const bool &muscl=false) : 
  SlopeLimiter(), M(m), MUSCL(muscl) {}

 ///Empty destructor
 virtual ~MinModLimiter() {}

 ///The basic minmod function
 double minmod(Vector<double> &args);


 ///The modified  minmod function
 double minmodB(Vector<double> &args, const double &h);

 ///The limit function
 void limit(const unsigned &i, 
            const Vector<DGElement*> &required_element_pt);
};



}
#endif