maxwell.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_maxwell_h
#define dealii_integrators_maxwell_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 Maxwell
  {
    template <int dim>
    Tensor<1,dim>
    curl_curl (
      const Tensor<2,dim> &h0,
      const Tensor<2,dim> &h1,
      const Tensor<2,dim> &h2)
    {
      Tensor<1,dim> result;
      switch (dim)
        {
        case 2:
          result[0] = h1[0][1]-h0[1][1];
          result[1] = h0[0][1]-h1[0][0];
          break;
        case 3:
          result[0] = h1[0][1]+h2[0][2]-h0[1][1]-h0[2][2];
          result[1] = h2[1][2]+h0[1][0]-h1[2][2]-h1[0][0];
          result[2] = h0[2][0]+h1[2][1]-h2[0][0]-h2[1][1];
          break;
        default:
          Assert(false, ExcNotImplemented());
        }
      return result;
    }

    template <int dim>
    Tensor<1,dim>
    tangential_curl (
      const Tensor<1,dim> &g0,
      const Tensor<1,dim> &g1,
      const Tensor<1,dim> &g2,
      const Tensor<1,dim> &normal)
    {
      Tensor<1,dim> result;

      switch (dim)
        {
        case 2:
          result[0] = normal[1] * (g1[0]-g0[1]);
          result[1] =-normal[0] * (g1[0]-g0[1]);
          break;
        case 3:
          result[0] = normal[2]*(g2[1]-g0[2])+normal[1]*(g1[0]-g0[1]);
          result[1] = normal[0]*(g0[2]-g1[0])+normal[2]*(g2[1]-g1[2]);
          result[2] = normal[1]*(g1[0]-g2[1])+normal[0]*(g0[2]-g2[0]);
          break;
        default:
          Assert(false, ExcNotImplemented());
        }
      return result;
    }

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

      AssertDimension(fe.get_fe().n_components(), dim);
      AssertDimension(M.m(), n_dofs);
      AssertDimension(M.n(), n_dofs);

      // Depending on the dimension,
      // the cross product is either
      // a scalar (2d) or a vector
      // (3d). Accordingly, in the
      // latter case we have to sum
      // up three bilinear forms, but
      // in 2d, we don't. Thus, we
      // need to adapt the loop over
      // all dimensions
      const unsigned int d_max = (dim==2) ? 1 : dim;

      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
        {
          const double dx = factor * fe.JxW(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<d_max; ++d)
                {
                  const unsigned int d1 = (d+1)%dim;
                  const unsigned int d2 = (d+2)%dim;

                  const double cv = fe.shape_grad_component(i,k,d2)[d1] - fe.shape_grad_component(i,k,d1)[d2];
                  const double cu = fe.shape_grad_component(j,k,d2)[d1] - fe.shape_grad_component(j,k,d1)[d2];

                  M(i,j) += dx * cu * cv;
                }
        }
    }

    template <int dim>
    void curl_matrix (
      FullMatrix<double> &M,
      const FEValuesBase<dim> &fe,
      const FEValuesBase<dim> &fetest,
      double factor = 1.)
    {
      const unsigned int n_dofs = fe.dofs_per_cell;
      const unsigned int t_dofs = fetest.dofs_per_cell;
      AssertDimension(fe.get_fe().n_components(), dim);
      // There should be the right number of components (3 in 3D, otherwise 1)
      // for the curl.
      AssertDimension(fetest.get_fe().n_components(), (dim == 3) ? dim : 1);
      AssertDimension(M.m(), t_dofs);
      AssertDimension(M.n(), n_dofs);

      const unsigned int d_max = (dim==2) ? 1 : dim;

      for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
        {
          const double dx = fe.JxW(k) * factor;
          for (unsigned int i=0; i<t_dofs; ++i)
            for (unsigned int j=0; j<n_dofs; ++j)
              for (unsigned int d=0; d<d_max; ++d)
                {
                  const unsigned int d1 = (d+1)%dim;
                  const unsigned int d2 = (d+2)%dim;

                  const double vv = fetest.shape_value_component(i,k,d);
                  const double cu = fe.shape_grad_component(j,k,d2)[d1] - fe.shape_grad_component(j,k,d1)[d2];
                  M(i,j) += dx * cu * vv;
                }
        }
    }

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

      AssertDimension(fe.get_fe().n_components(), dim);
      AssertDimension(M.m(), n_dofs);
      AssertDimension(M.n(), n_dofs);

      // Depending on the
      // dimension, the cross
      // product is either a scalar
      // (2d) or a vector
      // (3d). Accordingly, in the
      // latter case we have to sum
      // up three bilinear forms,
      // but in 2d, we don't. Thus,
      // we need to adapt the loop
      // over all dimensions
      const unsigned int d_max = (dim==2) ? 1 : dim;

      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 j=0; j<n_dofs; ++j)
              if (fe.get_fe().has_support_on_face(i, face_no) && fe.get_fe().has_support_on_face(j, face_no))
                {
                  for (unsigned int d=0; d<d_max; ++d)
                    {
                      const unsigned int d1 = (d+1)%dim;
                      const unsigned int d2 = (d+2)%dim;

                      const double cv = fe.shape_grad_component(i,k,d2)[d1] - fe.shape_grad_component(i,k,d1)[d2];
                      const double cu = fe.shape_grad_component(j,k,d2)[d1] - fe.shape_grad_component(j,k,d1)[d2];
                      const double v= fe.shape_value_component(i,k,d1)*n[d2] - fe.shape_value_component(i,k,d2)*n[d1];
                      const double u= fe.shape_value_component(j,k,d1)*n[d2] - fe.shape_value_component(j,k,d2)*n[d1];

                      M(i,j) += dx*(2.*penalty*u*v - cv*u - cu*v);
                    }
                }
        }
    }
    template <int dim>
    void tangential_trace_matrix (
      FullMatrix<double> &M,
      const FEValuesBase<dim> &fe,
      double factor = 1.)
    {
      const unsigned int n_dofs = fe.dofs_per_cell;

      AssertDimension(fe.get_fe().n_components(), dim);
      AssertDimension(M.m(), n_dofs);
      AssertDimension(M.n(), n_dofs);

      // Depending on the
      // dimension, the cross
      // product is either a scalar
      // (2d) or a vector
      // (3d). Accordingly, in the
      // latter case we have to sum
      // up three bilinear forms,
      // but in 2d, we don't. Thus,
      // we need to adapt the loop
      // over all dimensions
      const unsigned int d_max = (dim==2) ? 1 : dim;

      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 j=0; j<n_dofs; ++j)
              for (unsigned int d=0; d<d_max; ++d)
                {
                  const unsigned int d1 = (d+1)%dim;
                  const unsigned int d2 = (d+2)%dim;

                  const double v= fe.shape_value_component(i,k,d1)*n(d2) - fe.shape_value_component(i,k,d2)*n(d1);
                  const double u= fe.shape_value_component(j,k,d1)*n(d2) - fe.shape_value_component(j,k,d2)*n(d1);

                  M(i,j) += dx*u*v;
                }
        }
    }

    template <int dim>
    inline void ip_curl_matrix (
      FullMatrix<double> &M11,
      FullMatrix<double> &M12,
      FullMatrix<double> &M21,
      FullMatrix<double> &M22,
      const FEValuesBase<dim> &fe1,
      const FEValuesBase<dim> &fe2,
      const double pen,
      const double factor1 = 1.,
      const 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.m(), n_dofs);
      AssertDimension(M11.n(), n_dofs);
      AssertDimension(M12.m(), n_dofs);
      AssertDimension(M12.n(), n_dofs);
      AssertDimension(M21.m(), n_dofs);
      AssertDimension(M21.n(), n_dofs);
      AssertDimension(M22.m(), n_dofs);
      AssertDimension(M22.n(), n_dofs);

      const double nu1 = factor1;
      const double nu2 = (factor2 < 0) ? factor1 : factor2;
      const double penalty = .5 * pen * (nu1 + nu2);

      // Depending on the
      // dimension, the cross
      // product is either a scalar
      // (2d) or a vector
      // (3d). Accordingly, in the
      // latter case we have to sum
      // up three bilinear forms,
      // but in 2d, we don't. Thus,
      // we need to adapt the loop
      // over all dimensions
      const unsigned int d_max = (dim==2) ? 1 : dim;

      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)
            for (unsigned int j=0; j<n_dofs; ++j)
              for (unsigned int d=0; d<d_max; ++d)
                {
                  const unsigned int d1 = (d+1)%dim;
                  const unsigned int d2 = (d+2)%dim;
                  // curl u, curl v
                  const double cv1 = nu1*fe1.shape_grad_component(i,k,d2)[d1] - fe1.shape_grad_component(i,k,d1)[d2];
                  const double cv2 = nu2*fe2.shape_grad_component(i,k,d2)[d1] - fe2.shape_grad_component(i,k,d1)[d2];
                  const double cu1 = nu1*fe1.shape_grad_component(j,k,d2)[d1] - fe1.shape_grad_component(j,k,d1)[d2];
                  const double cu2 = nu2*fe2.shape_grad_component(j,k,d2)[d1] - fe2.shape_grad_component(j,k,d1)[d2];

                  // u x n, v x n
                  const double u1= fe1.shape_value_component(j,k,d1)*n(d2) - fe1.shape_value_component(j,k,d2)*n(d1);
                  const double u2=-fe2.shape_value_component(j,k,d1)*n(d2) + fe2.shape_value_component(j,k,d2)*n(d1);
                  const double v1= fe1.shape_value_component(i,k,d1)*n(d2) - fe1.shape_value_component(i,k,d2)*n(d1);
                  const double v2=-fe2.shape_value_component(i,k,d1)*n(d2) + fe2.shape_value_component(i,k,d2)*n(d1);

                  M11(i,j) += .5*dx*(2.*penalty*u1*v1 - cv1*u1 - cu1*v1);
                  M12(i,j) += .5*dx*(2.*penalty*v1*u2 - cv1*u2 - cu2*v1);
                  M21(i,j) += .5*dx*(2.*penalty*u1*v2 - cv2*u1 - cu1*v2);
                  M22(i,j) += .5*dx*(2.*penalty*u2*v2 - cv2*u2 - cu2*v2);
                }
        }
    }


  }
}


DEAL_II_NAMESPACE_CLOSE

#endif