libdeal.ii-dev 8.4.2-2+b1 / usr / include / deal.II / integrators / advection.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2010 - 2016 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__integrators_advection_h
#define dealii__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.
     *
     * @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);
                }
        }
    }



    /**
     * Scalar advection residual operator in strong form
     *
     * \f[ r_i = \int_Z  (\mathbf w \cdot \nabla)u\, v_i \, dx. \f]
     *
     * \warning This is not the residual consistent with cell_matrix(), but
     * with its transpose.
     */
    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][d]
                           * fe.shape_value(i,k) * velocity[d][k * v_increment];
        }
    }



    /**
     * Vector-valued advection residual operator in strong form
     *
     *
     * \f[ r_i = \int_Z \bigl((\mathbf w \cdot \nabla) \mathbf u\bigr)
     * \cdot\mathbf v_i \, dx. \f]
     *
     * \warning This is not the residual consistent with cell_matrix(), but
     * with its transpose.
     */
    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][d]
                             * fe.shape_value_component(i,k,c) * velocity[d][k * v_increment];
        }
    }



    /**
     * Scalar advection residual operator in weak form
     *
     * \f[ r_i = \int_Z  (\mathbf w \cdot \nabla)v\, u_i \, dx. \f]
     */
    template <int dim>
    inline void
    cell_residual  (
      Vector<double> &result,
      const FEValuesBase<dim> &fe,
      const std::vector<double> &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 \bigl((\mathbf w \cdot \nabla) \mathbf v\bigr)
     * \cdot\mathbf u_i \, dx. \f]
     */
    template <int dim>
    inline void
    cell_residual  (
      Vector<double> &result,
      const FEValuesBase<dim> &fe,
      const VectorSlice<const std::vector<std::vector<double> > > &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 the 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);
                  }
            }
        }
    }



    /**
     * Scalar case: Residual for upwind flux at the boundary for weak
     * advection operator. This is the value of the trial function at the
     * outflow boundary and the value of the incoming boundary condition on
     * the inflow boundary:
     * @f[
     * a_{ij} = \int_{\partial\Omega}
     * (\mathbf w\cdot\mathbf n)
     * \widehat u v_j \, ds
     * @f]
     *
     * Here, the numerical flux $\widehat u$ is the upwind value at the face,
     * namely the finite element function whose values are given in the
     * argument `input` on the outflow boundary. On the inflow boundary, it is
     * the inhomogenous boundary value in the argument `data`.
     *
     * 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 the 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>
    inline void
    upwind_value_residual(
      Vector<double> &result,
      const FEValuesBase<dim> &fe,
      const std::vector<double> &input,
      const std::vector<double> &data,
      const VectorSlice<const std::vector<std::vector<double> > > &velocity,
      double factor = 1.)
    {
      const unsigned int n_dofs = fe.dofs_per_cell;

      AssertDimension(input.size(), fe.n_quadrature_points);
      AssertDimension(data.size(), fe.n_quadrature_points);

      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];

          // Always use the upwind value
          const double val = (nv > 0.) ? input[k] : -data[k];

          for (unsigned i=0; i<n_dofs; ++i)
            {
              const double v= fe.shape_value(i,k);
              result(i) += dx * nv * val *v;
            }
        }
    }



    /**
     * Vector-valued case: Residual for upwind flux at the boundary for weak
     * advection operator. This is the value of the trial function at the
     * outflow boundary and the value of the incoming boundary condition on
     * the inflow boundary:
     * @f[
     * a_{ij} = \int_{\partial\Omega}
     * (\mathbf w\cdot\mathbf n)
     * \widehat u v_j \, ds
     * @f]
     *
     * Here, the numerical flux $\widehat u$ is the upwind value at the face,
     * namely the finite element function whose values are given in the
     * argument `input` on the outflow boundary. On the inflow boundary, it is
     * the inhomogenous boundary value in the argument `data`.
     *
     * 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 the 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>
    inline void
    upwind_value_residual(
      Vector<double> &result,
      const FEValuesBase<dim> &fe,
      const VectorSlice<const std::vector<std::vector<double> > > &input,
      const VectorSlice<const std::vector<std::vector<double> > > &data,
      const VectorSlice<const std::vector<std::vector<double> > > &velocity,
      double factor = 1.)
    {
      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);
      AssertVectorVectorDimension(data, n_comp, fe.n_quadrature_points);

      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];

          std::vector<double> val(n_comp);

          for (unsigned int d=0; d<n_comp; ++d)
            {
              val[d] = (nv > 0.) ? input[d][k] : -data[d][k];
              for (unsigned i=0; i<n_dofs; ++i)
                {
                  const double v= fe.shape_value_component(i,k,d);
                  result(i) += dx * nv * val[d] *v;
                }
            }
        }
    }



    /**
     * 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 the 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 * std::abs(nbeta) * fe1.JxW(k);
          FullMatrix<double> &M1 = nbeta > 0. ? M11 : M22;
          FullMatrix<double> &M2 = nbeta > 0. ? M21 : M12;
          const FEValuesBase<dim> &fe = nbeta > 0.  ? fe1 : fe2;
          const FEValuesBase<dim> &fetest = nbeta > 0.  ? fetest1 : fetest2;
          const FEValuesBase<dim> &fetestn = nbeta > 0.  ? fetest2 : fetest1;
          for (unsigned i=0; i<n1; ++i)
            for (unsigned j=0; j<n1; ++j)
              {
                if (fe1.get_fe().is_primitive())
                  {
                    M1(i,j) += dx_nbeta*fe.shape_value(j,k)*fetest.shape_value(i,k);
                    M2(i,j) -= dx_nbeta*fe.shape_value(j,k)*fetestn.shape_value(i,k);
                  }
                else
                  {
                    for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
                      {
                        M1(i,j) += dx_nbeta*fe.shape_value_component(j,k,d)*fetest.shape_value_component(i,k,d);
                        M2(i,j) -= dx_nbeta*fe.shape_value_component(j,k,d)*fetestn.shape_value_component(i,k,d);
                      }
                  }
              }
        }
    }



    /**
     * Scalar case: 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 the 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_face_residual (
      Vector<double> &result1,
      Vector<double> &result2,
      const FEValuesBase<dim> &fe1,
      const FEValuesBase<dim> &fe2,
      const std::vector<double> &input1,
      const std::vector<double> &input2,
      const VectorSlice<const std::vector<std::vector<double> > > &velocity,
      const double factor = 1.)
    {
      Assert(fe1.get_fe().n_components() == 1,
             ExcDimensionMismatch(fe1.get_fe().n_components(), 1));
      Assert(fe2.get_fe().n_components() == 1,
             ExcDimensionMismatch(fe2.get_fe().n_components(), 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)
            {
              const double v1 = fe1.shape_value(i,k);
              const double v2 = fe2.shape_value(i,k);
              const double u1 = input1[k];
              const double u2 = input2[k];
              if (nbeta > 0)
                {
                  result1(i) += dx_nbeta*u1*v1;
                  result2(i) -= dx_nbeta*u1*v2;
                }
              else
                {

                  result1(i) += dx_nbeta*u2*v1;
                  result2(i) -= dx_nbeta*u2*v2;
                }
            }
        }
    }



    /**
     * Vector-valued case: 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 the 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_face_residual (
      Vector<double> &result1,
      Vector<double> &result2,
      const FEValuesBase<dim> &fe1,
      const FEValuesBase<dim> &fe2,
      const VectorSlice<const std::vector<std::vector<double> > > &input1,
      const VectorSlice<const std::vector<std::vector<double> > > &input2,
      const VectorSlice<const std::vector<std::vector<double> > > &velocity,
      const double factor = 1.)
    {
      const unsigned int n_comp = fe1.get_fe().n_components();
      const unsigned int n1 = fe1.dofs_per_cell;
      AssertVectorVectorDimension(input1, n_comp, fe1.n_quadrature_points);
      AssertVectorVectorDimension(input2, n_comp, fe2.n_quadrature_points);

      // 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 int d=0; d<n_comp; ++d)
              {
                const double v1 = fe1.shape_value_component(i,k,d);
                const double v2 = fe2.shape_value_component(i,k,d);
                const double u1 = input1[d][k];
                const double u2 = input2[d][k];
                if (nbeta > 0)
                  {
                    result1(i) += dx_nbeta*u1*v1;
                    result2(i) -= dx_nbeta*u1*v2;
                  }
                else
                  {

                    result1(i) += dx_nbeta*u2*v1;
                    result2(i) -= dx_nbeta*u2*v2;
                  }
              }
        }
    }

  }
}


DEAL_II_NAMESPACE_CLOSE

#endif