From 85af76e941c5e11aafe11136f99d5b06fecf805b Mon Sep 17 00:00:00 2001
From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Tue, 4 Sep 2012 08:42:26 +0000
Subject: [PATCH] Implement advection operators in weak form

git-svn-id: https://svn.dealii.org/trunk@26230 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../include/deal.II/integrators/advection.h   | 390 ++++++++++++++++++
 1 file changed, 390 insertions(+)
 create mode 100644 deal.II/include/deal.II/integrators/advection.h

diff --git a/deal.II/include/deal.II/integrators/advection.h b/deal.II/include/deal.II/integrators/advection.h
new file mode 100644
index 0000000000..56cdbbcdd4
--- /dev/null
+++ b/deal.II/include/deal.II/integrators/advection.h
@@ -0,0 +1,390 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2010, 2011, 2012 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__integrators_advection_h
+#define __deal2__integrators_advection_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/fe/mapping.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/meshworker/dof_info.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LocalIntegrators
+{
+/**
+ * @brief Local integrators related to advection along a vector field and its DG formulations
+ *
+ * All advection operators depend on an advection velocity denoted by
+ * <b>w</b> in the formulas below. It is denoted as <tt>velocity</tt>
+ * in the parameter lists.
+ *
+ * The functions cell_matrix() and both upwind_value_matrix() are
+ * taking the equation in weak form, that is, the directional
+ * derivative is on the test function.
+ *
+ * @ingroup Integrators
+ * @author Guido Kanschat
+ * @date 2012
+ */
+  namespace Advection
+  {
+/**
+ * Advection along the direction <b>w</b> in weak form
+ * with derivative on the test function
+ * \f[
+ * m_{ij} = \int_Z u_j\,(\mathbf w \cdot \nabla) v_i \, dx.
+ * \f]
+ *
+ * The FiniteElement in <tt>fe</tt> may be scalar or vector valued. In
+ * the latter case, the advection operator is applied to each component
+ * separately.
+ *
+ * @param M: The advection matrix obtained as result
+ * @param fe: The FEValues object describing the local trial function
+ * space. #update_values and #update_gradients, and #update_JxW_values
+ * must be set.
+ * @param fetest: The FEValues object describing the local test
+ * function space. #update_values and #update_gradients must be set.
+ * @param velocity: The advection velocity, a vector of dimension
+ * <tt>dim</tt>. Each component may either contain a vector of length
+ * one, in which case a constant velocity is assumed, or a vector with
+ * as many entries as quadrature points if the velocity is not constant.
+ * @param factor is an optional multiplication factor for the result.
+ *
+ * @ingroup Integrators
+ * @author Guido Kanschat
+ * @date 2012
+ */
+    template<int dim>
+      void cell_matrix (
+	FullMatrix<double>& M,
+	const FEValuesBase<dim>& fe,
+	const FEValuesBase<dim>& fetest,
+	const VectorSlice<const std::vector<std::vector<double> > >& velocity,
+	const double factor = 1.)
+      {
+	const unsigned int n_dofs = fe.dofs_per_cell;
+	const unsigned int t_dofs = fetest.dofs_per_cell;
+	const unsigned int n_components = fe.get_fe().n_components();
+	
+	AssertDimension(velocity.size(), dim);
+					 // If the size of the
+					 // velocity vectors is one,
+					 // then do not increment
+					 // between quadrature points.
+	const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
+	
+	if (v_increment == 1)
+	  {
+	    AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
+	  }
+	
+	AssertDimension(M.n(), n_dofs);
+	AssertDimension(M.m(), t_dofs);
+      
+	for (unsigned k=0;k<fe.n_quadrature_points;++k)
+	  {
+	    const double dx = factor * fe.JxW(k);
+	    const unsigned int vindex = k * v_increment;
+	      
+	    for (unsigned j=0;j<n_dofs;++j)
+	      for (unsigned i=0;i<t_dofs;++i)
+		for (unsigned int c=0;c<n_components;++c)
+		  {
+		    double wgradv = velocity[0][vindex]
+				    * fe.shape_grad_component(i,k,c)[0];
+		    for (unsigned int d=1;d<dim;++d)
+		      wgradv += velocity[d][vindex]
+				* fe.shape_grad_component(i,k,c)[d];
+		    M(i,j) -= dx * wgradv * fe.shape_value_component(j,k,c);
+		  }
+	  }
+      }
+
+/**
+ * Advection residual operator in weak form
+ *
+ * \f[
+ * r_i = \int_Z  u\,(\mathbf w \cdot \nabla) v_i \, dx.
+ * \f]
+ */
+    template <int dim>
+      inline void
+      cell_residual  (
+	Vector<double>& result,
+	const FEValuesBase<dim>& fe,
+	const std::vector<Tensor<1,dim> >& input,
+	const VectorSlice<const std::vector<std::vector<double> > >& velocity,
+	double factor = 1.)
+      {
+	const unsigned int nq = fe.n_quadrature_points;
+	const unsigned int n_dofs = fe.dofs_per_cell;
+	Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq));
+	Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs));
+      
+	AssertDimension(velocity.size(), dim);
+	const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
+	if (v_increment == 1)
+	  {
+	  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
+	  }
+	
+	for (unsigned k=0;k<nq;++k)
+	  {
+	    const double dx = factor * fe.JxW(k);	      
+	    for (unsigned i=0;i<n_dofs;++i)
+	      for (unsigned int d=0;d<dim;++d)
+		result(i) += dx * input[k]
+			     * fe.shape_grad(i,k)[d] * velocity[d][k * v_increment];
+	  }
+      }
+    
+
+/**
+ * Vector-valued advection residual operator in weak form
+ *
+ *
+ * \f[
+ * r_i = \int_Z  \mathbf u\cdot\bigl((\mathbf w \cdot \nabla) \mathbf v_i\bigr) \, dx.
+ * \f]
+ */
+    template <int dim>
+      inline void
+      cell_residual  (
+	Vector<double>& result,
+	const FEValuesBase<dim>& fe,
+	const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > >& input,
+	const VectorSlice<const std::vector<std::vector<double> > >& velocity,
+	double factor = 1.)
+      {
+	const unsigned int nq = fe.n_quadrature_points;
+	const unsigned int n_dofs = fe.dofs_per_cell;      
+	const unsigned int n_comp = fe.get_fe().n_components();
+      
+	AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points);
+	Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs));
+      
+	AssertDimension(velocity.size(), dim);
+	const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
+	if (v_increment == 1)
+	  {
+	  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
+	  }
+	
+	for (unsigned k=0;k<nq;++k)
+	  {
+	    const double dx = factor * fe.JxW(k);
+	    for (unsigned i=0;i<n_dofs;++i)
+	      for (unsigned int c=0;c<n_comp;++c)
+		for (unsigned int d=0;d<dim;++d)
+		  result(i) += dx * input[c][k]
+			       * fe.shape_grad_component(i,k,c)[d] * velocity[d][k * v_increment];
+	  }
+      }
+    
+				     /**
+				      * Upwind flux at the boundary
+				      * for weak advection
+				      * operator. This is the value of
+				      * the trial function at the
+				      * outflow boundary and zero
+				      * else:
+				      * @f[
+				      * a_{ij} = \int_{\partial\Omega}
+				      * [\mathbf w\cdot\mathbf n]_+
+				      * u_i v_j \, ds
+				      * @f]
+				      *
+				      * The <tt>velocity</tt> is
+				      * provided as a VectorSlice,
+				      * having <tt>dim</tt> vectors,
+				      * one for each velocity
+				      * component. Each of the
+				      * vectors must either have only
+				      * a single entry, if t he
+				      * advection velocity is
+				      * constant, or have an entry
+				      * for each quadrature point.
+				      *
+				      * The finite element can have
+				      * several components, in which
+				      * case each component is
+				      * advected by the same velocity.
+				      */
+    template <int dim>
+      void upwind_value_matrix(
+	FullMatrix<double>& M,
+	const FEValuesBase<dim>& fe,
+	const FEValuesBase<dim>& fetest,
+	const VectorSlice<const std::vector<std::vector<double> > >& velocity,
+	double factor = 1.)
+      {
+	const unsigned int n_dofs = fe.dofs_per_cell; 
+	const unsigned int t_dofs = fetest.dofs_per_cell; 
+	unsigned int n_components = fe.get_fe().n_components();
+	AssertDimension (M.m(), n_dofs);
+	AssertDimension (M.n(), n_dofs);
+	
+	AssertDimension(velocity.size(), dim);
+	const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
+	if (v_increment == 1)
+	  {
+	  AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points);
+	  }	
+	
+	for (unsigned k=0;k<fe.n_quadrature_points;++k)
+	  {
+	    const double dx = factor * fe.JxW(k);
+	      
+	    double nv = 0.;
+	    for (unsigned int d=0;d<dim;++d)
+	      nv += fe.normal_vector(k)[d] * velocity[d][k * v_increment];
+	      
+	    if (nv > 0)
+	      {
+		for (unsigned i=0;i<t_dofs;++i)
+		  for (unsigned j=0;j<n_dofs;++j)
+		    {
+		      if (fe.get_fe().is_primitive())
+			M(i,j) += dx * nv * fe.shape_value(i,k) * fe.shape_value(j,k);
+		      else
+			for (unsigned int c=0;c<n_components;++c)
+			  M(i,j) += dx * nv * fetest.shape_value_component(i,k,c)
+				    * fe.shape_value_component(j,k,c);
+		    }
+	      }
+	  }
+      }
+
+    				     /**
+				      * Upwind flux in the interior
+				      * for weak advection
+				      * operator. Matrix entries
+				      * correspond to the upwind value
+				      * of the trial function, multiplied
+				      * by the jump of the test
+				      * functions
+				      * @f[
+				      * a_{ij} = \int_F \left|\mathbf w
+				      * \cdot \mathbf n\right|
+				      * u^\uparrow
+				      * (v^\uparrow-v^\downarrow)
+				      * \,ds
+				      * @f]
+				      *
+				      * The <tt>velocity</tt> is
+				      * provided as a VectorSlice,
+				      * having <tt>dim</tt> vectors,
+				      * one for each velocity
+				      * component. Each of the
+				      * vectors must either have only
+				      * a single entry, if t he
+				      * advection velocity is
+				      * constant, or have an entry
+				      * for each quadrature point.
+				      *
+				      * The finite element can have
+				      * several components, in which
+				      * case each component is
+				      * advected the same way.
+				      */
+    template <int dim>
+      void upwind_value_matrix (
+	FullMatrix<double>& M11,
+	FullMatrix<double>& M12,
+	FullMatrix<double>& M21,
+	FullMatrix<double>& M22,
+	const FEValuesBase<dim>& fe1,
+	const FEValuesBase<dim>& fe2,
+	const FEValuesBase<dim>& fetest1,
+	const FEValuesBase<dim>& fetest2,
+	const VectorSlice<const std::vector<std::vector<double> > >& velocity,
+	const double factor = 1.)
+      {
+	const unsigned int n1 = fe1.dofs_per_cell;
+					 // Multiply the quadrature point
+					 // index below with this factor to
+					 // have simpler data for constant
+					 // velocities.
+	AssertDimension(velocity.size(), dim);
+	const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1;
+	if (v_increment == 1)
+	  {
+	    AssertVectorVectorDimension(velocity, dim, fe1.n_quadrature_points);
+	  }
+	
+	for (unsigned k=0;k<fe1.n_quadrature_points;++k)
+	  {
+	    double nbeta = fe1.normal_vector(k)[0] * velocity[0][k * v_increment];
+	    for (unsigned int d=1;d<dim;++d)
+	      nbeta += fe1.normal_vector(k)[d] * velocity[d][k * v_increment];
+	    const double dx_nbeta = factor * nbeta * fe1.JxW(k);
+	      
+	    for (unsigned i=0;i<n1;++i)
+	      for (unsigned j=0;j<n1;++j)
+		if (fe1.get_fe().is_primitive())
+		  {
+		    if (nbeta > 0)
+		      {
+			M11(i,j) += dx_nbeta
+				    * fe1.shape_value(j,k)
+				    * fetest1.shape_value(i,k);
+			M21(i,j) -= dx_nbeta
+				    * fe1.shape_value(j,k)
+				    * fetest2.shape_value(i,k);
+		      }
+		    else
+		      {
+			M22(i,j) -= dx_nbeta
+				    * fe2.shape_value(j,k)
+				    * fetest2.shape_value(i,k);
+			M12(i,j) += dx_nbeta
+				    * fe2.shape_value(j,k)
+				    * fetest1.shape_value(i,k);
+		      }
+		  }
+		else
+		  {
+		    for (unsigned int d=0;d<fe1.get_fe().n_components();++d)
+		      if (nbeta > 0)
+			{
+			  M11(i,j) += dx_nbeta
+				      * fe1.shape_value_component(j,k,d)
+				      * fetest1.shape_value_component(i,k,d);
+			  M21(i,j) -= dx_nbeta
+				      * fe1.shape_value_component(j,k,d)
+				      * fetest2.shape_value_component(i,k,d);
+			}
+		      else
+			{
+			  M22(i,j) -= dx_nbeta
+				      * fe2.shape_value_component(j,k,d)
+				      * fetest2.shape_value_component(i,k,d);
+			  M12(i,j) += dx_nbeta
+				      * fe2.shape_value_component(j,k,d)
+				      * fetest1.shape_value_component(i,k,d);
+			}	      
+		  }
+	  }
+      }
+  }
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
-- 
2.39.5