electrodynamic.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_PDELAB_LOCALOPERATOR_ELECTRODYNAMIC_HH
#define DUNE_PDELAB_LOCALOPERATOR_ELECTRODYNAMIC_HH
 
#include <cstddef>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
#include <dune/common/deprecated.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
 
#include <dune/pdelab/common/quadraturerules.hh>
#include <dune/pdelab/localoperator/numericalresidual.hh>
#include <dune/pdelab/localoperator/pattern.hh>
#include <dune/pdelab/localoperator/flags.hh>
 
namespace Dune {
  namespace PDELab {
 
    namespace ElectrodynamicImpl {
 
 
      // Curl manipulation
 
      // dimension of the curl for a given space dimension
      constexpr std::size_t dimOfCurl(std::size_t dimOfSpace)
      {
        return
          dimOfSpace == 1 ? 2 :
          dimOfSpace == 2 ? 1 :
          dimOfSpace == 3 ? 3 :
          // Dune exceptions are difficult to construct in constexpr functions
          throw std::invalid_argument("Only applicable for dimensions 1-3");
      }
 
      template<typename RF>
      void jacobianToCurl(FieldVector<RF, 1> &curl,
                          const FieldMatrix<RF, 2, 2> &jacobian)
      {
        curl[0] = jacobian[1][0] - jacobian[0][1];
      }
      template<typename RF>
      void jacobianToCurl(FieldVector<RF, 3> &curl,
                          const FieldMatrix<RF, 3, 3> &jacobian)
      {
        for(unsigned i = 0; i < 3; ++i)
          curl[i] = jacobian[(i+2)%3][(i+1)%3] - jacobian[(i+1)%3][(i+2)%3];
      }
 
    } // namespace ElectrodynamicImpl
 
 
    template<typename Eps>
    class Electrodynamic_T
      : public FullVolumePattern
      , public LocalOperatorDefaultFlags
      , public JacobianBasedAlphaVolume<Electrodynamic_T<Eps> >
    {
 
    public:
 
      // pattern assembly flags
      static constexpr bool doPatternVolume = true;
      static constexpr bool doAlphaVolume = true;
 
 
      Electrodynamic_T(const Eps &eps, int qorder = 2) :
        eps_(eps), qorder_(qorder)
      {}
 
      Electrodynamic_T(Eps &&eps, int qorder = 2) :
        eps_(std::move(eps)), qorder_(qorder)
      {}
 
      template<typename EG, typename LFS, typename X, typename M>
      void jacobian_volume (const EG& eg, const LFS& lfsu, const X& x,
                            const LFS& lfsv, M& mat) const
      {
        using BasisTraits =
          typename LFS::Traits::FiniteElementType::Traits::Basis::Traits;
 
        // static checks
        static constexpr unsigned dimR = BasisTraits::dimRange;
        static_assert(dimR == 3 || dimR == 2, "Works only in 2D or 3D");
 
        using ctype = typename EG::Geometry::ctype;
        using DF = typename BasisTraits::DomainField;
        static_assert(std::is_same<ctype, DF>::value, "Grids ctype and "
                      "Finite Elements DomainFieldType must match");
 
        using Range = typename BasisTraits::Range;
        std::vector<Range> phi(lfsu.size());
 
        // loop over quadrature points
        for(const auto &qp : quadratureRule(eg.geometry(), qorder_))
        {
          // values of basefunctions
          lfsu.finiteElement().basis().evaluateFunction(qp.position(),phi);
 
          // calculate T
          auto factor = qp.weight()
            * eg.geometry().integrationElement(qp.position())
            * eps_(eg.entity(), qp.position());
 
          for(unsigned i = 0; i < lfsu.size(); ++i)
            for(unsigned j = 0; j < lfsu.size(); ++j)
              mat.accumulate(lfsv,i,lfsu,j,factor * (phi[i] * phi[j]));
        }
      }
 
    private:
      Eps eps_;
      int qorder_;
    };
 
 
    template<class Eps>
    Electrodynamic_T<std::decay_t<Eps> >
    makeLocalOperatorEdynT(Eps &&eps, int qorder = 2)
    {
      return { std::forward<Eps>(eps), qorder };
    }
 
 
    template<typename Mu>
    class Electrodynamic_S
      : public FullVolumePattern
      , public LocalOperatorDefaultFlags
      , public JacobianBasedAlphaVolume<Electrodynamic_S<Mu> >
    {
 
    public:
 
      // pattern assembly flags
      static constexpr bool doPatternVolume = true;
      static constexpr bool doAlphaVolume = true;
 
 
      Electrodynamic_S(const Mu &mu, int qorder = 2) :
        mu_(mu), qorder_(qorder)
      {}
 
      Electrodynamic_S(Mu &&mu, int qorder = 2) :
        mu_(std::move(mu)), qorder_(qorder)
      {}
 
      template<typename EG, typename LFS, typename X, typename M>
      void jacobian_volume (const EG& eg, const LFS& lfsu, const X& x,
                            const LFS& lfsv, M& mat) const
      {
        using ElectrodynamicImpl::dimOfCurl;
        using ElectrodynamicImpl::jacobianToCurl;
 
        using BasisTraits =
          typename LFS::Traits::FiniteElementType::Traits::Basis::Traits;
 
        // static checks
        static constexpr unsigned dimR = BasisTraits::dimRange;
        static_assert(dimR == 3 || dimR == 2, "Works only in 2D or 3D");
 
        using ctype = typename EG::Geometry::ctype;
        using DF = typename BasisTraits::DomainField;
        static_assert(std::is_same<ctype, DF>::value, "Grids ctype and "
                      "Finite Elements DomainFieldType must match");
 
        using Jacobian = typename BasisTraits::Jacobian;
        std::vector<Jacobian> J(lfsu.size());
 
        using RF = typename BasisTraits::RangeField;
        using Curl = FieldVector<RF, dimOfCurl(dimR)>;
        std::vector<Curl> rotphi(lfsu.size());
 
        // loop over quadrature points
        for(const auto &qp : quadratureRule(eg.geometry(), qorder_))
        {
          // curl of the basefunctions
          lfsu.finiteElement().basis().evaluateJacobian(qp.position(),J);
          for(unsigned i = 0; i < lfsu.size(); ++i)
            jacobianToCurl(rotphi[i], J[i]);
 
          // calculate S
          auto factor = qp.weight()
            * eg.geometry().integrationElement(qp.position())
            / mu_(eg.entity(), qp.position());
 
          for(unsigned i = 0; i < lfsu.size(); ++i)
            for(unsigned j = 0; j < lfsu.size(); ++j)
              mat.accumulate(lfsv,i,lfsu,j,factor * (rotphi[i] * rotphi[j]));
 
        }
      }
 
    private:
      Mu mu_;
      int qorder_;
    };
 
 
    template<class Mu>
    Electrodynamic_S<std::decay_t<Mu> >
    makeLocalOperatorEdynS(Mu &&mu, int qorder = 2)
    {
      return { std::forward<Mu>(mu), qorder };
    }
 
  } // namespace PDELab
} // namespace Dune
 
#endif // DUNE_PDELAB_LOCALOPERATOR_ELECTRODYNAMIC_HH