<commit_before>/* $Id$ */
/* Author: <NAME>, University of Heidelberg, 2001 */

/*    $Id$       */
/*    Version: $Name$                                          */
/*                                                                */
/*    Copyright (C) 2001, 2002, 2003 by the deal.II authors */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyright and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

				 // The first few files have already
				 // been covered in previous examples
				 // and will thus not be further
				 // commented on.
#include <base/quadrature_lib.h>
#include <base/function.h>
#include <lac/vector.h>
#include <lac/sparse_matrix.h>
#include <lac/vector_memory.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/grid_out.h>
#include <grid/grid_refinement.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <fe/fe_values.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <numerics/data_out.h>

				 // This is the first new file. It
				 // declares the ``MappingQ1'' class
				 // that gives the standard bilinear
				 // mapping. For bilinear mappings use
				 // an object of this class rather
				 // than an object of the
				 // ``MappingQ(1)'' class, as the
				 // ``MappingQ1'' class is optimized
				 // due to the pre-knowledge of the
				 // actual polynomial degree 1.
#include <fe/mapping_q1.h>
				 // Here the discontinuous finite
				 // elements are defined. They are
				 // used as all other finite elements.
#include <fe/fe_dgq.h>
				 // We are going to use the simplest
				 // possible solver, called Richardson
				 // iteration, that represents a simple
				 // defect correction. This, in
				 // combination with a block SSOR
				 // preconditioner (defined in
				 // precondition_block.h), that uses
				 // the special block matrix structur
				 // of system matrices arising from DG
				 // discretizations.
#include <lac/solver_richardson.h>
#include <lac/precondition_block.h>
				 // We are going to use gradients as
				 // refinement indicator.
#include <numerics/derivative_approximation.h>
				 // Finally we do some time comparison
				 // using the ``Timer'' class.
#include <base/timer.h>

				 // And this again is C++:
#include <iostream>
#include <fstream>


				 // @sect3{Equation data}
				 //
				 // First we define the classes
				 // representing the equation-specific
				 // functions. Both classes, ``RHS''
				 // and ``BoundaryValues'', are
				 // derived from the ``Function''
				 // class. Only the ``value_list''
				 // function are implemented because
				 // only lists of function values are
				 // computed rather than single
				 // values.
template <int dim>
class RHS:  public Function<dim>
{
  public:
    virtual void value_list (const std::vector<Point<dim> > &points,
			     std::vector<double> &values,
			     const unsigned int component=0) const;
};


template <int dim>
class BoundaryValues:  public Function<dim>
{
  public:
    virtual void value_list (const std::vector<Point<dim> > &points,
			     std::vector<double> &values,
			     const unsigned int component=0) const;
};


				 // The class ``Beta'' represents the
				 // vector valued flow field of the
				 // linear transport equation and is
				 // not derived from the ``Function''
				 // class as we prefer to get function
				 // values of type ``Point'' rather
				 // than of type
				 // ``Vector<double>''. This, because
				 // there exist scalar products
				 // between ``Point'' and ``Point'' as
				 // well as between ``Point'' and
				 // ``Tensor'', simplifying terms like
				 // $\beta\cdot n$ and
				 // $\beta\cdot\nabla v$.
                                 //
                                 // An unnecessary empty constructor
                                 // is added to the class to work
                                 // around a bug in Compaq's cxx
                                 // compiler which otherwise reports
                                 // an error about an omitted
                                 // initializer for an object of
                                 // this class further down.
template <int dim>
class Beta
{
  public:
    Beta () {};
    void value_list (const std::vector<Point<dim> > &points,
		     std::vector<Point<dim> > &values) const;
};


				 // The implementation of the
				 // ``value_list'' functions of these
				 // classes are rather simple.  For
				 // simplicity the right hand side is
				 // set to be zero but will be
				 // assembled anyway.
template <int dim>
void RHS<dim>::value_list(const std::vector<Point<dim> > &points,
			  std::vector<double> &values,
			  const unsigned int) const
{
				   // Usually we check whether input
				   // parameter have the right sizes.
  Assert(values.size()==points.size(),
	 ExcDimensionMismatch(values.size(),points.size()));

  for (unsigned int i=0; i<values.size(); ++i)
    values[i]=0;
}


				 // The flow field is chosen to be
				 // circular, counterclockwise, and with
				 // the origin as midpoint.
template <int dim>
void Beta<dim>::value_list(const std::vector<Point<dim> > &points,
			   std::vector<Point<dim> > &values) const
{
  Assert(values.size()==points.size(),
	 ExcDimensionMismatch(values.size(),points.size()));

  for (unsigned int i=0; i<points.size(); ++i)
    {
      const Point<dim> &p=points[i];
      Point<dim> &beta=values[i];

      beta(0) = -p(1);
      beta(1) = p(0);
      beta /= std::sqrt(beta.square());
    }
}


				 // Hence the inflow boundary of the
				 // unit square [0,1]^2 are the right
				 // and the lower boundaries. We
				 // prescribe discontinuous boundary
				 // values 1 and 0 on the x-axis and
				 // value 0 on the right boundary. The
				 // values of this function on the
				 // outflow boundaries will not be
				 // used within the DG scheme.
template <int dim>
void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
				       std::vector<double> &values,
				       const unsigned int) const
{
  Assert(values.size()==points.size(),
	 ExcDimensionMismatch(values.size(),points.size()));

  for (unsigned int i=0; i<values.size(); ++i)
    {
      if (points[i](0)<0.5)
	values[i]=1.;
      else
	values[i]=0.;
    }
}


				 // @sect3{Class: DGTransportEquation}
				 //
				 // Next we define the
				 // equation-dependent and
				 // DG-method-dependent class
				 // ``DGTransportEquation''. Its
				 // member functions were already
				 // mentioned in the Introduction and
				 // will be explained
				 // below. Furthermore it includes
				 // objects of the previously defined
				 // ``Beta'', ``RHS'' and
				 // ``BoundaryValues'' function
				 // classes.
template <int dim>
class DGTransportEquation
{
  public:
    DGTransportEquation();

    void assemble_cell_term(const FEValues<dim>& fe_v,
			    FullMatrix<double> &u_v_matrix,
			    Vector<double> &cell_vector) const;
    
    void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
				FullMatrix<double> &u_v_matrix,
				Vector<double> &cell_vector) const;
    
    void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
			     const FEFaceValuesBase<dim>& fe_v_neighbor,
			     FullMatrix<double> &u_v_matrix,
			     FullMatrix<double> &un_v_matrix) const;

    void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
			     const FEFaceValuesBase<dim>& fe_v_neighbor,
			     FullMatrix<double> &u_v_matrix,
			     FullMatrix<double> &un_v_matrix,
			     FullMatrix<double> &u_vn_matrix,
			     FullMatrix<double> &un_vn_matrix) const;
  private:
    const Beta<dim> beta_function;
    const RHS<dim> rhs_function;
    const BoundaryValues<dim> boundary_function;
};


template <int dim>
DGTransportEquation<dim>::DGTransportEquation ()
		:
		beta_function (),
		rhs_function (),
		boundary_function ()
{}


				 // @sect4{Function: assemble_cell_term}
				 //
				 // The ``assemble_cell_term''
				 // function assembles the cell terms
				 // of the discretization.
				 // ``u_v_matrix'' is a cell matrix,
				 // i.e. for a DG method of degree 1,
				 // it is of size 4 times 4, and
				 // ``cell_vector'' is of size 4.
				 // When this function is invoked,
				 // ``fe_v'' is already reinit'ed with the
				 // current cell before and includes
				 // all shape values needed.
template <int dim>
void DGTransportEquation<dim>::assemble_cell_term(
  const FEValues<dim> &fe_v,
  FullMatrix<double> &u_v_matrix,
  Vector<double> &cell_vector) const
{
				   // First we ask ``fe_v'' for the
				   // shape gradients, shape values and
				   // quadrature weights,
  const std::vector<double> &JxW = fe_v.get_JxW_values ();

				   // Then the flow field beta and the
				   // ``rhs_function'' are evaluated at
				   // the quadrature points,
  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
  std::vector<double> rhs (fe_v.n_quadrature_points);
  
  beta_function.value_list (fe_v.get_quadrature_points(), beta);
  rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
  
				   // and the cell matrix and cell
				   // vector are assembled due to the
				   // terms $-(u,\beta\cdot\nabla
				   // v)_K$ and $(f,v)_K$.
  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
    for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
      {
	for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
	  u_v_matrix(i,j) -= beta[point]*fe_v.shape_grad(i,point)*
			      fe_v.shape_value(j,point) *
			      JxW[point];
	
	cell_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point];
      }
}


				 // @sect4{Function: assemble_boundary_term}
				 //
				 // The ``assemble_boundary_term''
				 // function assembles the face terms
				 // at boundary faces.  When this
				 // function is invoked, ``fe_v'' is
				 // already reinit'ed with the current
				 // cell and current face. Hence it
				 // provides the shape values on that
				 // boundary face.
template <int dim>
void DGTransportEquation<dim>::assemble_boundary_term(
  const FEFaceValues<dim>& fe_v,    
  FullMatrix<double> &u_v_matrix,
  Vector<double> &cell_vector) const
{
				   // First we check whether the
				   // current face is really at the
				   // boundary.
  Assert(fe_v.get_face()->at_boundary(), ExcInternalError());
  
				   // Again, as in the previous
				   // function, we ask the ``FEValues''
				   // object for the shape values and
				   // the quadrature weights
  const std::vector<double> &JxW = fe_v.get_JxW_values ();
				   // but here also for the normals.
  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();

				   // We evaluate the flow field
				   // and the boundary values at the
				   // quadrature points.
  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
  std::vector<double> g(fe_v.n_quadrature_points);
  
  beta_function.value_list (fe_v.get_quadrature_points(), beta);
  boundary_function.value_list (fe_v.get_quadrature_points(), g);

				   // Then we assemble cell vector and
				   // cell matrix according to the DG
				   // method given in the
				   // introduction.
  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
    {
      const double beta_n=beta[point] * normals[point];      
					 // We assemble the term
					 // $(\beta\cdot n
					 // u,v)_{\partial K_+}$,
      if (beta_n>0)
	for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
	  for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
	    u_v_matrix(i,j) += beta_n *
			       fe_v.shape_value(j,point) *
			       fe_v.shape_value(i,point) *
			       JxW[point];
      else
					 // and the term $(\beta\cdot
					 // n g,v)_{\partial
					 // K_-\cap\partial\Omega}$,
	for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
	  cell_vector(i) -= beta_n *
			    g[point] *
			    fe_v.shape_value(i,point) *
			    JxW[point];
    }
}


				 // @sect4{Function: assemble_face_term1}
				 //
				 // The ``assemble_face_term1''
				 // function assembles the face terms
				 // corresponding to the first version
				 // of the DG method, cf. above. For
				 // that case, the face terms are
				 // given as a sum of integrals over
				 // all cell boundaries.
				 //
				 // When this function is invoked,
				 // ``fe_v'' and ``fe_v_neighbor'' are
				 // already reinit'ed with the current
				 // cell and the neighoring cell,
				 // respectively, as well as with the
				 // current face. Hence they provide
				 // the inner and outer shape values
				 // on the face.
				 //
				 // In addition to the cell matrix
				 // ``u_v_matrix'' this function has
				 // got a new argument
				 // ``un_v_matrix'', that stores
				 // contributions to the system matrix
				 // that are based on outer values of
				 // u, see $\hat u_h$ in the
				 // introduction, and inner values of
				 // v, see $v_h$. Here we note that
				 // ``un'' is the short notation for
				 // ``u_neighbor'' and represents
				 // $\hat u_h$.
template <int dim>
void DGTransportEquation<dim>::assemble_face_term1(
  const FEFaceValuesBase<dim>& fe_v,
  const FEFaceValuesBase<dim>& fe_v_neighbor,      
  FullMatrix<double> &u_v_matrix,
  FullMatrix<double> &un_v_matrix) const
{
				   // First we check that the current
				   // face is not at the boundary by
				   // accident.
  Assert(!fe_v.get_face()->at_boundary(), ExcInternalError());
  
				   // Again, as in the previous
				   // function, we ask the FEValues
				   // objects for the shape values,
				   // the quadrature weights and the
				   // normals
  const std::vector<double> &JxW = fe_v.get_JxW_values ();
  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();

				   // and we evaluate the flow field
				   // at the quadrature points.
  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
  beta_function.value_list (fe_v.get_quadrature_points(), beta);

				   // Then we assemble the cell
				   // matrices according to the DG
				   // method given in the
				   // introduction.
  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
    {
      const double beta_n=beta[point] * normals[point];
					 // We assemble the term
					 // $(\beta\cdot n
					 // u,v)_{\partial K_+}$,
      if (beta_n>0)
	for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
	  for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
	    u_v_matrix(i,j) += beta_n *
			       fe_v.shape_value(j,point) *
			       fe_v.shape_value(i,point) *
			       JxW[point];<commit_msg>Don't use get_face. Previous use may have been honorable, but that function is now gone and we can either add another parameter to the functions that use this, or leave out the check. Since the functions that use this are called in only one place, it seems acceptable to remove the check, although this is arguably not the best way...<commit_after>/* $Id$ */
/* Author: <NAME>, University of Heidelberg, 2001 */

/*    $Id$       */
/*    Version: $Name$                                          */
/*                                                                */
/*    Copyright (C) 2001, 2002, 2003 by the deal.II authors */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyright and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

				 // The first few files have already
				 // been covered in previous examples
				 // and will thus not be further
				 // commented on.
#include <base/quadrature_lib.h>
#include <base/function.h>
#include <lac/vector.h>
#include <lac/sparse_matrix.h>
#include <lac/vector_memory.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/grid_out.h>
#include <grid/grid_refinement.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <fe/fe_values.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <numerics/data_out.h>

				 // This is the first new file. It
				 // declares the ``MappingQ1'' class
				 // that gives the standard bilinear
				 // mapping. For bilinear mappings use
				 // an object of this class rather
				 // than an object of the
				 // ``MappingQ(1)'' class, as the
				 // ``MappingQ1'' class is optimized
				 // due to the pre-knowledge of the
				 // actual polynomial degree 1.
#include <fe/mapping_q1.h>
				 // Here the discontinuous finite
				 // elements are defined. They are
				 // used as all other finite elements.
#include <fe/fe_dgq.h>
				 // We are going to use the simplest
				 // possible solver, called Richardson
				 // iteration, that represents a simple
				 // defect correction. This, in
				 // combination with a block SSOR
				 // preconditioner (defined in
				 // precondition_block.h), that uses
				 // the special block matrix structur
				 // of system matrices arising from DG
				 // discretizations.
#include <lac/solver_richardson.h>
#include <lac/precondition_block.h>
				 // We are going to use gradients as
				 // refinement indicator.
#include <numerics/derivative_approximation.h>
				 // Finally we do some time comparison
				 // using the ``Timer'' class.
#include <base/timer.h>

				 // And this again is C++:
#include <iostream>
#include <fstream>


				 // @sect3{Equation data}
				 //
				 // First we define the classes
				 // representing the equation-specific
				 // functions. Both classes, ``RHS''
				 // and ``BoundaryValues'', are
				 // derived from the ``Function''
				 // class. Only the ``value_list''
				 // function are implemented because
				 // only lists of function values are
				 // computed rather than single
				 // values.
template <int dim>
class RHS:  public Function<dim>
{
  public:
    virtual void value_list (const std::vector<Point<dim> > &points,
			     std::vector<double> &values,
			     const unsigned int component=0) const;
};


template <int dim>
class BoundaryValues:  public Function<dim>
{
  public:
    virtual void value_list (const std::vector<Point<dim> > &points,
			     std::vector<double> &values,
			     const unsigned int component=0) const;
};


				 // The class ``Beta'' represents the
				 // vector valued flow field of the
				 // linear transport equation and is
				 // not derived from the ``Function''
				 // class as we prefer to get function
				 // values of type ``Point'' rather
				 // than of type
				 // ``Vector<double>''. This, because
				 // there exist scalar products
				 // between ``Point'' and ``Point'' as
				 // well as between ``Point'' and
				 // ``Tensor'', simplifying terms like
				 // $\beta\cdot n$ and
				 // $\beta\cdot\nabla v$.
                                 //
                                 // An unnecessary empty constructor
                                 // is added to the class to work
                                 // around a bug in Compaq's cxx
                                 // compiler which otherwise reports
                                 // an error about an omitted
                                 // initializer for an object of
                                 // this class further down.
template <int dim>
class Beta
{
  public:
    Beta () {};
    void value_list (const std::vector<Point<dim> > &points,
		     std::vector<Point<dim> > &values) const;
};


				 // The implementation of the
				 // ``value_list'' functions of these
				 // classes are rather simple.  For
				 // simplicity the right hand side is
				 // set to be zero but will be
				 // assembled anyway.
template <int dim>
void RHS<dim>::value_list(const std::vector<Point<dim> > &points,
			  std::vector<double> &values,
			  const unsigned int) const
{
				   // Usually we check whether input
				   // parameter have the right sizes.
  Assert(values.size()==points.size(),
	 ExcDimensionMismatch(values.size(),points.size()));

  for (unsigned int i=0; i<values.size(); ++i)
    values[i]=0;
}


				 // The flow field is chosen to be
				 // circular, counterclockwise, and with
				 // the origin as midpoint.
template <int dim>
void Beta<dim>::value_list(const std::vector<Point<dim> > &points,
			   std::vector<Point<dim> > &values) const
{
  Assert(values.size()==points.size(),
	 ExcDimensionMismatch(values.size(),points.size()));

  for (unsigned int i=0; i<points.size(); ++i)
    {
      const Point<dim> &p=points[i];
      Point<dim> &beta=values[i];

      beta(0) = -p(1);
      beta(1) = p(0);
      beta /= std::sqrt(beta.square());
    }
}


				 // Hence the inflow boundary of the
				 // unit square [0,1]^2 are the right
				 // and the lower boundaries. We
				 // prescribe discontinuous boundary
				 // values 1 and 0 on the x-axis and
				 // value 0 on the right boundary. The
				 // values of this function on the
				 // outflow boundaries will not be
				 // used within the DG scheme.
template <int dim>
void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
				       std::vector<double> &values,
				       const unsigned int) const
{
  Assert(values.size()==points.size(),
	 ExcDimensionMismatch(values.size(),points.size()));

  for (unsigned int i=0; i<values.size(); ++i)
    {
      if (points[i](0)<0.5)
	values[i]=1.;
      else
	values[i]=0.;
    }
}


				 // @sect3{Class: DGTransportEquation}
				 //
				 // Next we define the
				 // equation-dependent and
				 // DG-method-dependent class
				 // ``DGTransportEquation''. Its
				 // member functions were already
				 // mentioned in the Introduction and
				 // will be explained
				 // below. Furthermore it includes
				 // objects of the previously defined
				 // ``Beta'', ``RHS'' and
				 // ``BoundaryValues'' function
				 // classes.
template <int dim>
class DGTransportEquation
{
  public:
    DGTransportEquation();

    void assemble_cell_term(const FEValues<dim>& fe_v,
			    FullMatrix<double> &u_v_matrix,
			    Vector<double> &cell_vector) const;
    
    void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
				FullMatrix<double> &u_v_matrix,
				Vector<double> &cell_vector) const;
    
    void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
			     const FEFaceValuesBase<dim>& fe_v_neighbor,
			     FullMatrix<double> &u_v_matrix,
			     FullMatrix<double> &un_v_matrix) const;

    void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
			     const FEFaceValuesBase<dim>& fe_v_neighbor,
			     FullMatrix<double> &u_v_matrix,
			     FullMatrix<double> &un_v_matrix,
			     FullMatrix<double> &u_vn_matrix,
			     FullMatrix<double> &un_vn_matrix) const;
  private:
    const Beta<dim> beta_function;
    const RHS<dim> rhs_function;
    const BoundaryValues<dim> boundary_function;
};


template <int dim>
DGTransportEquation<dim>::DGTransportEquation ()
		:
		beta_function (),
		rhs_function (),
		boundary_function ()
{}


				 // @sect4{Function: assemble_cell_term}
				 //
				 // The ``assemble_cell_term''
				 // function assembles the cell terms
				 // of the discretization.
				 // ``u_v_matrix'' is a cell matrix,
				 // i.e. for a DG method of degree 1,
				 // it is of size 4 times 4, and
				 // ``cell_vector'' is of size 4.
				 // When this function is invoked,
				 // ``fe_v'' is already reinit'ed with the
				 // current cell before and includes
				 // all shape values needed.
template <int dim>
void DGTransportEquation<dim>::assemble_cell_term(
  const FEValues<dim> &fe_v,
  FullMatrix<double> &u_v_matrix,
  Vector<double> &cell_vector) const
{
				   // First we ask ``fe_v'' for the
				   // shape gradients, shape values and
				   // quadrature weights,
  const std::vector<double> &JxW = fe_v.get_JxW_values ();

				   // Then the flow field beta and the
				   // ``rhs_function'' are evaluated at
				   // the quadrature points,
  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
  std::vector<double> rhs (fe_v.n_quadrature_points);
  
  beta_function.value_list (fe_v.get_quadrature_points(), beta);
  rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
  
				   // and the cell matrix and cell
				   // vector are assembled due to the
				   // terms $-(u,\beta\cdot\nabla
				   // v)_K$ and $(f,v)_K$.
  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
    for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
      {
	for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
	  u_v_matrix(i,j) -= beta[point]*fe_v.shape_grad(i,point)*
			      fe_v.shape_value(j,point) *
			      JxW[point];
	
	cell_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point];
      }
}


				 // @sect4{Function: assemble_boundary_term}
				 //
				 // The ``assemble_boundary_term''
				 // function assembles the face terms
				 // at boundary faces.  When this
				 // function is invoked, ``fe_v'' is
				 // already reinit'ed with the current
				 // cell and current face. Hence it
				 // provides the shape values on that
				 // boundary face.
template <int dim>
void DGTransportEquation<dim>::assemble_boundary_term(
  const FEFaceValues<dim>& fe_v,    
  FullMatrix<double> &u_v_matrix,
  Vector<double> &cell_vector) const
{
				   // Again, as in the previous
				   // function, we ask the ``FEValues''
				   // object for the shape values and
				   // the quadrature weights
  const std::vector<double> &JxW = fe_v.get_JxW_values ();
				   // but here also for the normals.
  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();

				   // We evaluate the flow field
				   // and the boundary values at the
				   // quadrature points.
  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
  std::vector<double> g(fe_v.n_quadrature_points);
  
  beta_function.value_list (fe_v.get_quadrature_points(), beta);
  boundary_function.value_list (fe_v.get_quadrature_points(), g);

				   // Then we assemble cell vector and
				   // cell matrix according to the DG
				   // method given in the
				   // introduction.
  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
    {
      const double beta_n=beta[point] * normals[point];      
					 // We assemble the term
					 // $(\beta\cdot n
					 // u,v)_{\partial K_+}$,
      if (beta_n>0)
	for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
	  for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
	    u_v_matrix(i,j) += beta_n *
			       fe_v.shape_value(j,point) *
			       fe_v.shape_value(i,point) *
			       JxW[point];
      else
					 // and the term $(\beta\cdot
					 // n g,v)_{\partial
					 // K_-\cap\partial\Omega}$,
	for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
	  cell_vector(i) -= beta_n *
			    g[point] *
			    fe_v.shape_value(i,point) *
			    JxW[point];
    }
}


				 // @sect4{Function: assemble_face_term1}
				 //
				 // The ``assemble_face_term1''
				 // function assembles the face terms
				 // corresponding to the first version
				 // of the DG method, cf. above. For
				 // that case, the face terms are
				 // given as a sum of integrals over
				 // all cell boundaries.
				 //
				 // When this function is invoked,
				 // ``fe_v'' and ``fe_v_neighbor'' are
				 // already reinit'ed with the current
				 // cell and the neighoring cell,
				 // respectively, as well as with the
				 // current face. Hence they provide
				 // the inner and outer shape values
				 // on the face.
				 //
				 // In addition to the cell matrix
				 // ``u_v_matrix'' this function has
				 // got a new argument
				 // ``un_v_matrix'', that stores
				 // contributions to the system matrix
				 // that are based on outer values of
				 // u, see $\hat u_h$ in the
				 // introduction, and inner values of
				 // v, see $v_h$. Here we note that
				 // ``un'' is the short notation for
				 // ``u_neighbor'' and represents
				 // $\hat u_h$.
template <int dim>
void DGTransportEquation<dim>::assemble_face_term1(
  const FEFaceValuesBase<dim>& fe_v,
  const FEFaceValuesBase<dim>& fe_v_neighbor,      
  FullMatrix<double> &u_v_matrix,
  FullMatrix<double> &un_v_matrix) const
{
				   // Again, as in the previous
				   // function, we ask the FEValues
				   // objects for the shape values,
				   // the quadrature weights and the
				   // normals
  const std::vector<double> &JxW = fe_v.get_JxW_values ();
  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();

				   // and we evaluate the flow field
				   // at the quadrature points.
  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
  beta_function.value_list (fe_v.get_quadrature_points(), beta);

				   // Then we assemble the cell
				   // matrices according to the DG
				   // method given in the
				   // introduction.
  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
    {
      const double beta_n=beta[point] * normals[point];
					 // We assemble the term
					 // $(\beta\cdot n
					 // u,v)_{\partial K_+}$,
      if (beta_n>0)
	for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
	  for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
	    u_v_matrix(i,j) += beta_n *
			       fe_v.shape_value(j,point) *
			       fe_v.shape_value(i,point) *
			       JxW[point];
{'fraction_non_alphanumeric': 0.12060794800622843, 'fraction_numerical': 0.004603615191930133, 'mean_word_length': 2.319065273564768, 'pattern_counts': {'":': 0, '<': 255, '<?xml version=': 0, '>': 233, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4993073', 'n_tokens_mistral': 10775, 'n_tokens_neox': 8717, 'n_words': 2726}