laplace.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2010 - 2017 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_laplace_h
#define dealii_integrators_laplace_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
{
  namespace Laplace
  {
    template <int dim>
    void cell_matrix (
      FullMatrix<double> &M,
      const FEValuesBase<dim> &fe,
      const double factor = 1.)
    {
      const unsigned int n_dofs = fe.dofs_per_cell;
      const unsigned int n_components = fe.get_fe().n_components();

      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
        {
          const double dx = fe.JxW(k) * factor;
          for (unsigned int i=0; i<n_dofs; ++i)
            {
              double Mii = 0.0;
              for (unsigned int d=0; d<n_components; ++d)
                Mii += dx *
                       (fe.shape_grad_component(i,k,d) * fe.shape_grad_component(i,k,d));

              M(i,i) += Mii;

              for (unsigned int j=i+1; j<n_dofs; ++j)
                {
                  double Mij = 0.0;
                  for (unsigned int d=0; d<n_components; ++d)
                    Mij += dx *
                           (fe.shape_grad_component(j,k,d) * fe.shape_grad_component(i,k,d));

                  M(i,j) += Mij;
                  M(j,i) += Mij;
                }
            }
        }

    }

    template <int dim>
    inline void
    cell_residual  (
      Vector<double> &result,
      const FEValuesBase<dim> &fe,
      const std::vector<Tensor<1,dim> > &input,
      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));

      for (unsigned int k=0; k<nq; ++k)
        {
          const double dx = factor * fe.JxW(k);
          for (unsigned int i=0; i<n_dofs; ++i)
            result(i) += dx * (input[k] * fe.shape_grad(i,k));
        }
    }


    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,
      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));

      for (unsigned int k=0; k<nq; ++k)
        {
          const double dx = factor * fe.JxW(k);
          for (unsigned int i=0; i<n_dofs; ++i)
            for (unsigned int d=0; d<n_comp; ++d)
              {

                result(i) += dx * (input[d][k] * fe.shape_grad_component(i,k,d));
              }
        }
    }


    template <int dim>
    void nitsche_matrix (
      FullMatrix<double> &M,
      const FEValuesBase<dim> &fe,
      double penalty,
      double factor = 1.)
    {
      const unsigned int n_dofs = fe.dofs_per_cell;
      const unsigned int n_comp = fe.get_fe().n_components();

      Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs));
      Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs));

      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
        {
          const double dx = fe.JxW(k) * factor;
          const Tensor<1,dim> n = fe.normal_vector(k);
          for (unsigned int i=0; i<n_dofs; ++i)
            for (unsigned int j=0; j<n_dofs; ++j)
              for (unsigned int d=0; d<n_comp; ++d)
                M(i,j) += dx *
                          (2. * fe.shape_value_component(i,k,d) * penalty * fe.shape_value_component(j,k,d)
                           - (n * fe.shape_grad_component(i,k,d)) * fe.shape_value_component(j,k,d)
                           - (n * fe.shape_grad_component(j,k,d)) * fe.shape_value_component(i,k,d));
        }
    }

    template <int dim>
    void nitsche_tangential_matrix (
      FullMatrix<double> &M,
      const FEValuesBase<dim> &fe,
      double penalty,
      double factor = 1.)
    {
      const unsigned int n_dofs = fe.dofs_per_cell;
      AssertDimension(fe.get_fe().n_components(), dim);
      Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs));
      Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs));

      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
        {
          const double dx = fe.JxW(k) * factor;
          const Tensor<1,dim> n = fe.normal_vector(k);
          for (unsigned int i=0; i<n_dofs; ++i)
            for (unsigned int j=0; j<n_dofs; ++j)
              {
                double udotn = 0.;
                double vdotn = 0.;
                double ngradun = 0.;
                double ngradvn = 0.;

                for (unsigned int d=0; d<dim; ++d)
                  {
                    udotn += n[d]*fe.shape_value_component(j,k,d);
                    vdotn += n[d]*fe.shape_value_component(i,k,d);
                    ngradun += n*fe.shape_grad_component(j,k,d)*n[d];
                    ngradvn += n*fe.shape_grad_component(i,k,d)*n[d];
                  }

                for (unsigned int d=0; d<dim; ++d)
                  {
                    const double v_t = fe.shape_value_component(i,k,d)-vdotn*n[d];
                    const double dnv_t = n * fe.shape_grad_component(i,k,d)-ngradvn*n[d];
                    const double u_t = fe.shape_value_component(j,k,d)-udotn*n[d];
                    const double dnu_t = n * fe.shape_grad_component(j,k,d)-ngradun*n[d];

                    M(i,j) += dx *(2.*penalty*u_t*v_t-dnu_t*v_t-dnv_t*u_t);
                  }
              }
        }
    }

    template <int dim>
    void nitsche_residual (
      Vector<double> &result,
      const FEValuesBase<dim> &fe,
      const std::vector<double> &input,
      const std::vector<Tensor<1,dim> > &Dinput,
      const std::vector<double> &data,
      double penalty,
      double factor = 1.)
    {
      const unsigned int n_dofs = fe.dofs_per_cell;
      AssertDimension(input.size(), fe.n_quadrature_points);
      AssertDimension(Dinput.size(), fe.n_quadrature_points);
      AssertDimension(data.size(), fe.n_quadrature_points);

      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
        {
          const double dx = factor * fe.JxW(k);
          const Tensor<1,dim> n = fe.normal_vector(k);
          for (unsigned int i=0; i<n_dofs; ++i)
            {
              const double dnv = fe.shape_grad(i,k) * n;
              const double dnu = Dinput[k] * n;
              const double v= fe.shape_value(i,k);
              const double u= input[k];
              const double g= data[k];

              result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
            }
        }
    }

    template <int dim>
    void nitsche_residual (
      Vector<double> &result,
      const FEValuesBase<dim> &fe,
      const VectorSlice<const std::vector<std::vector<double> > > &input,
      const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput,
      const VectorSlice<const std::vector<std::vector<double> > > &data,
      double penalty,
      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(Dinput, n_comp, fe.n_quadrature_points);
      AssertVectorVectorDimension(data, n_comp, fe.n_quadrature_points);

      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
        {
          const double dx = factor * fe.JxW(k);
          const Tensor<1,dim> n = fe.normal_vector(k);
          for (unsigned int i=0; i<n_dofs; ++i)
            for (unsigned int d=0; d<n_comp; ++d)
              {
                const double dnv = fe.shape_grad_component(i,k,d) * n;
                const double dnu = Dinput[d][k] * n;
                const double v= fe.shape_value_component(i,k,d);
                const double u= input[d][k];
                const double g= data[d][k];

                result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
              }
        }
    }

    template <int dim>
    void ip_matrix (
      FullMatrix<double> &M11,
      FullMatrix<double> &M12,
      FullMatrix<double> &M21,
      FullMatrix<double> &M22,
      const FEValuesBase<dim> &fe1,
      const FEValuesBase<dim> &fe2,
      double penalty,
      double factor1 = 1.,
      double factor2 = -1.)
    {
      const unsigned int n_dofs = fe1.dofs_per_cell;
      AssertDimension(M11.n(), n_dofs);
      AssertDimension(M11.m(), n_dofs);
      AssertDimension(M12.n(), n_dofs);
      AssertDimension(M12.m(), n_dofs);
      AssertDimension(M21.n(), n_dofs);
      AssertDimension(M21.m(), n_dofs);
      AssertDimension(M22.n(), n_dofs);
      AssertDimension(M22.m(), n_dofs);

      const double nui = factor1;
      const double nue = (factor2 < 0) ? factor1 : factor2;
      const double nu = .5*(nui+nue);

      for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
        {
          const double dx = fe1.JxW(k);
          const Tensor<1,dim> n = fe1.normal_vector(k);
          for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
            {
              for (unsigned int i=0; i<n_dofs; ++i)
                {
                  for (unsigned int j=0; j<n_dofs; ++j)
                    {
                      const double vi = fe1.shape_value_component(i,k,d);
                      const double dnvi = n * fe1.shape_grad_component(i,k,d);
                      const double ve = fe2.shape_value_component(i,k,d);
                      const double dnve = n * fe2.shape_grad_component(i,k,d);
                      const double ui = fe1.shape_value_component(j,k,d);
                      const double dnui = n * fe1.shape_grad_component(j,k,d);
                      const double ue = fe2.shape_value_component(j,k,d);
                      const double dnue = n * fe2.shape_grad_component(j,k,d);
                      M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
                      M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-nu*penalty*vi*ue);
                      M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-nu*penalty*ui*ve);
                      M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+nu*penalty*ue*ve);
                    }
                }
            }
        }
    }

    template <int dim>
    void ip_tangential_matrix (
      FullMatrix<double> &M11,
      FullMatrix<double> &M12,
      FullMatrix<double> &M21,
      FullMatrix<double> &M22,
      const FEValuesBase<dim> &fe1,
      const FEValuesBase<dim> &fe2,
      double penalty,
      double factor1 = 1.,
      double factor2 = -1.)
    {
      const unsigned int n_dofs = fe1.dofs_per_cell;
      AssertDimension(fe1.get_fe().n_components(), dim);
      AssertDimension(fe2.get_fe().n_components(), dim);
      AssertDimension(M11.n(), n_dofs);
      AssertDimension(M11.m(), n_dofs);
      AssertDimension(M12.n(), n_dofs);
      AssertDimension(M12.m(), n_dofs);
      AssertDimension(M21.n(), n_dofs);
      AssertDimension(M21.m(), n_dofs);
      AssertDimension(M22.n(), n_dofs);
      AssertDimension(M22.m(), n_dofs);

      const double nui = factor1;
      const double nue = (factor2 < 0) ? factor1 : factor2;
      const double nu = .5*(nui+nue);

      for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
        {
          const double dx = fe1.JxW(k);
          const Tensor<1,dim> n = fe1.normal_vector(k);
          for (unsigned int i=0; i<n_dofs; ++i)
            {
              // We compute the tangential component by subtracting
              // the normal component from the total. Thus, compute
              // the normal component first
              for (unsigned int j=0; j<n_dofs; ++j)
                {
                  double u1dotn = 0.;
                  double v1dotn = 0.;
                  double u2dotn = 0.;
                  double v2dotn = 0.;

                  double ngradu1n = 0.;
                  double ngradv1n = 0.;
                  double ngradu2n = 0.;
                  double ngradv2n = 0.;

                  for (unsigned int d=0; d<dim; ++d)
                    {
                      u1dotn += n[d]*fe1.shape_value_component(j,k,d);
                      v1dotn += n[d]*fe1.shape_value_component(i,k,d);
                      u2dotn += n[d]*fe2.shape_value_component(j,k,d);
                      v2dotn += n[d]*fe2.shape_value_component(i,k,d);

                      ngradu1n += n*fe1.shape_grad_component(j,k,d)*n[d];
                      ngradv1n += n*fe1.shape_grad_component(i,k,d)*n[d];
                      ngradu2n += n*fe2.shape_grad_component(j,k,d)*n[d];
                      ngradv2n += n*fe2.shape_grad_component(i,k,d)*n[d];
                    }

                  // The following code is equal to ip_matrix() with
                  // the only exception that all variables introduced
                  // below denote tangential components.
                  for (unsigned int d=0; d<dim; ++d)
                    {
                      const double vi = fe1.shape_value_component(i,k,d)-v1dotn*n[d];
                      const double dnvi = n * fe1.shape_grad_component(i,k,d)-ngradv1n*n[d];

                      const double ve = fe2.shape_value_component(i,k,d)-v2dotn*n[d];
                      const double dnve = n * fe2.shape_grad_component(i,k,d)-ngradv2n*n[d];

                      const double ui = fe1.shape_value_component(j,k,d)-u1dotn*n[d];
                      const double dnui = n * fe1.shape_grad_component(j,k,d)-ngradu1n*n[d];

                      const double ue = fe2.shape_value_component(j,k,d)-u2dotn*n[d];
                      const double dnue = n * fe2.shape_grad_component(j,k,d)-ngradu2n*n[d];

                      M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
                      M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-nu*penalty*vi*ue);
                      M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-nu*penalty*ui*ve);
                      M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+nu*penalty*ue*ve);
                    }
                }
            }
        }
    }

    template <int dim>
    void
    ip_residual(
      Vector<double> &result1,
      Vector<double> &result2,
      const FEValuesBase<dim> &fe1,
      const FEValuesBase<dim> &fe2,
      const std::vector<double> &input1,
      const std::vector<Tensor<1,dim> > &Dinput1,
      const std::vector<double> &input2,
      const std::vector<Tensor<1,dim> > &Dinput2,
      double pen,
      double int_factor = 1.,
      double ext_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 double nui = int_factor;
      const double nue = (ext_factor < 0) ? int_factor : ext_factor;
      const double penalty = .5 * pen * (nui + nue);

      const unsigned int n_dofs = fe1.dofs_per_cell;

      for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
        {
          const double dx = fe1.JxW(k);
          const Tensor<1,dim> n = fe1.normal_vector(k);

          for (unsigned int i=0; i<n_dofs; ++i)
            {
              const double vi = fe1.shape_value(i,k);
              const Tensor<1,dim> &Dvi = fe1.shape_grad(i,k);
              const double dnvi = Dvi * n;
              const double ve = fe2.shape_value(i,k);
              const Tensor<1,dim> &Dve = fe2.shape_grad(i,k);
              const double dnve = Dve * n;

              const double ui = input1[k];
              const Tensor<1,dim> &Dui = Dinput1[k];
              const double dnui = Dui * n;
              const double ue = input2[k];
              const Tensor<1,dim> &Due = Dinput2[k];
              const double dnue = Due * n;

              result1(i) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi);
              result1(i) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue);
              result2(i) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve);
              result2(i) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve);
            }
        }
    }


    template <int dim>
    void
    ip_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<Tensor<1,dim> > > > &Dinput1,
      const VectorSlice<const std::vector<std::vector<double> > > &input2,
      const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput2,
      double pen,
      double int_factor = 1.,
      double ext_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(Dinput1, n_comp, fe1.n_quadrature_points);
      AssertVectorVectorDimension(input2, n_comp, fe2.n_quadrature_points);
      AssertVectorVectorDimension(Dinput2, n_comp, fe2.n_quadrature_points);

      const double nui = int_factor;
      const double nue = (ext_factor < 0) ? int_factor : ext_factor;
      const double penalty = .5 * pen * (nui + nue);


      for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
        {
          const double dx = fe1.JxW(k);
          const Tensor<1,dim> n = fe1.normal_vector(k);

          for (unsigned int i=0; i<n1; ++i)
            for (unsigned int d=0; d<n_comp; ++d)
              {
                const double vi = fe1.shape_value_component(i,k,d);
                const Tensor<1,dim> &Dvi = fe1.shape_grad_component(i,k,d);
                const double dnvi = Dvi * n;
                const double ve = fe2.shape_value_component(i,k,d);
                const Tensor<1,dim> &Dve = fe2.shape_grad_component(i,k,d);
                const double dnve = Dve * n;

                const double ui = input1[d][k];
                const Tensor<1,dim> &Dui = Dinput1[d][k];
                const double dnui = Dui * n;
                const double ue = input2[d][k];
                const Tensor<1,dim> &Due = Dinput2[d][k];
                const double dnue = Due * n;

                result1(i) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi);
                result1(i) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue);
                result2(i) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve);
                result2(i) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve);
              }
        }
    }



    template <int dim, int spacedim, typename number>
    double compute_penalty(
      const MeshWorker::DoFInfo<dim,spacedim,number> &dinfo1,
      const MeshWorker::DoFInfo<dim,spacedim,number> &dinfo2,
      unsigned int deg1,
      unsigned int deg2)
    {
      const unsigned int normal1 = GeometryInfo<dim>::unit_normal_direction[dinfo1.face_number];
      const unsigned int normal2 = GeometryInfo<dim>::unit_normal_direction[dinfo2.face_number];
      const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1+1);
      const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2+1);

      double penalty1 = deg1sq / dinfo1.cell->extent_in_direction(normal1);
      double penalty2 = deg2sq / dinfo2.cell->extent_in_direction(normal2);
      if (dinfo1.cell->has_children() ^ dinfo2.cell->has_children())
        {
          Assert (dinfo1.face == dinfo2.face, ExcInternalError());
          Assert (dinfo1.face->has_children(), ExcInternalError());
          penalty1 *= 2;
        }
      const double penalty = 0.5*(penalty1 + penalty2);
      return penalty;
    }
  }
}


DEAL_II_NAMESPACE_CLOSE

#endif