ddr-magnetostatics.hpp

Go to the documentation of this file.

// Implementation of the discrete de Rham sequence for the magnetostatic problem.
//
// Author: Daniele Di Pietro (daniele.di-pietro@umontpellier.fr)
//
 
/*
 * The DDR scheme is described in
 *
 *  An arbitrary-order method for magnetostatics on polyhedral meshes based on a discrete de Rham
 sequence. 
 *   D. A. Di Pietro and J. Droniou, 31p, 2020. url: https://arxiv.org/abs/2005.06890.
 *
 * If you use this code in a scientific publication, please mention the above article.
 *
 */
 
#ifndef MAGNETOSTATICS_HPP
#define MAGNETOSTATICS_HPP
 
#include <iostream>
 
#include <boost/math/constants/constants.hpp>
 
#include <Eigen/Sparse>
#include <unsupported/Eigen/SparseExtra>
 
#include <mesh.hpp>
#include <mesh_builder.hpp>
#include <local_static_condensation.hpp>
#include <linearsolver.hpp>
 
#include <xdiv.hpp>
#include <xcurl.hpp>
 
namespace HArDCore3D
{
 
  struct Magnetostatics
  {
    typedef Eigen::SparseMatrix<double> SystemMatrixType;
    
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> ForcingTermType;
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> SolutionPotentialType;
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> SolutionCurlType;
    typedef IntegralWeight PermeabilityType;
 
    Magnetostatics(
                   const DDRCore & ddrcore,           
                   bool use_threads,                  
                   std::ostream & output = std::cout  
                   );
 
    void assembleLinearSystem(
                              const ForcingTermType & f,      
                              const PermeabilityType & mu,    
                              const SolutionPotentialType & u 
                              );
 
 
    inline size_t dimensionSpace() const
    {
      return m_xcurl.dimension() + m_xdiv.dimension();
    }
 
    inline size_t nbSCDOFs() const
    {
      return m_ddrcore.mesh().n_cells() * m_nloc_sc_H;
    }
 
    inline size_t sizeSystem() const
    {
      return dimensionSpace() - nbSCDOFs();
    }
 
    inline const XCurl & xCurl() const
    {
      return m_xcurl;
    }
 
    inline const XDiv & xDiv() const
    {
      return m_xdiv;
    }
 
    inline const SystemMatrixType & systemMatrix() const {
      return m_A;
    }
 
    inline SystemMatrixType & systemMatrix() {
      return m_A;
    }
 
    inline const Eigen::VectorXd & systemVector() const {
      return m_b;
    }
 
    inline Eigen::VectorXd & systemVector() {
      return m_b;
    }
 
    inline const double & stabilizationParameter() const {
      return m_stab_par;
    }
 
    inline double & stabilizationParameter() {
      return m_stab_par;
    }
 
    inline const SystemMatrixType & scMatrix() const {
      return m_sc_A;
    }
 
    inline Eigen::VectorXd & scVector() {
      return m_sc_b;
    }
 
    std::vector<double> computeNorms(
                       const std::vector<Eigen::VectorXd> & list_dofs 
                      ) const;
 
  private:
    std::pair<Eigen::MatrixXd, Eigen::VectorXd>
    _compute_local_contribution(
                                size_t iT,                      
                                const ForcingTermType & f,      
                                const PermeabilityType & mu     
                                );
 
    LocalStaticCondensation _compute_static_condensation(const size_t & iT) const;
 
    void _assemble_local_contribution(
                                      size_t iT,                                               
                                      const std::pair<Eigen::MatrixXd, Eigen::VectorXd> & lsT, 
                                      std::list<Eigen::Triplet<double> > & A1,                 
                                      Eigen::VectorXd & b1,                                    
                                      std::list<Eigen::Triplet<double> > & A2,                 
                                      Eigen::VectorXd & b2                                     
                                      );
    
    const DDRCore & m_ddrcore;
    bool m_use_threads;
    std::ostream & m_output;
    XCurl m_xcurl;
    XDiv m_xdiv;
    const size_t m_nloc_sc_H;   // Number of statically condensed DOFs of H in each cell
    SystemMatrixType m_A;   // Matrix and RHS of statically condensed system
    Eigen::VectorXd m_b;
    SystemMatrixType m_sc_A; // Static condensation operator and RHS (to recover statically condensed DOFs)
    Eigen::VectorXd m_sc_b;
    double m_stab_par;
  };
 
  //------------------------------------------------------------------------------
  // Exact solutions
  //------------------------------------------------------------------------------
 
  static const double PI = boost::math::constants::pi<double>();
  using std::sin;
  using std::cos;
 
  //------------------------------------------------------------------------------
  
  static Magnetostatics::SolutionPotentialType
  constant_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                 return Eigen::Vector3d(1., 2., 3.);
               };
 
  static Magnetostatics::SolutionCurlType  
  constant_sigma = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                     return Eigen::Vector3d::Zero();
                   };
 
  static Magnetostatics::ForcingTermType
  constant_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                 return Eigen::Vector3d::Zero();
               };
 
  static Magnetostatics::PermeabilityType
  constant_mu = Magnetostatics::PermeabilityType(1.);
  
  //------------------------------------------------------------------------------
 
  static Magnetostatics::SolutionPotentialType
  linear_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
               return Eigen::Vector3d(
                                      x(0) + x(1) + x(2),
                                      2.*(x(0)-x(1)+x(2)),
                                      -x(0) - x(1) + x(2)
                                      );
             };
 
  static Magnetostatics::SolutionCurlType  
  linear_sigma = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                   return Eigen::Vector3d(-3., 2., 1.);
                 };
 
  static Magnetostatics::ForcingTermType
  linear_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
               return Eigen::Vector3d::Zero();
             };
 
  static Magnetostatics::PermeabilityType
  linear_mu = Magnetostatics::PermeabilityType(1.);
  
  //------------------------------------------------------------------------------
 
  static Magnetostatics::SolutionPotentialType
  trigonometric_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return Eigen::Vector3d(
                                             cos(PI*x(0)) * sin(PI*x(1)) * sin(PI*x(2)),
                                             -2. * sin(PI*x(0)) * cos(PI*x(1)) * sin(PI*x(2)),
                                             sin(PI*x(0)) * sin(PI*x(1)) * cos(PI*x(2))
                                             );
                    };
 
  static Magnetostatics::SolutionCurlType  
  trigonometric_sigma = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                          return 3. * PI * Eigen::Vector3d(
                                                           sin(PI*x(0)) * cos(PI*x(1)) * cos(PI*x(2)),
                                                           0.,
                                                           -cos(PI*x(0)) * cos(PI*x(1)) * sin(PI*x(2))
                                                           );
                        };
 
  static Magnetostatics::ForcingTermType
  trigonometric_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return 3. * std::pow(PI, 2) * Eigen::Vector3d(
                                                                    cos(PI*x(0)) * sin(PI*x(1)) * sin(PI*x(2)),
                                                                    -2. * sin(PI*x(0)) * cos(PI*x(1)) * sin(PI*x(2)),
                                                                    sin(PI*x(0)) * sin(PI*x(1)) * cos(PI*x(2))
                                                                    );
                    };
 
  static Magnetostatics::PermeabilityType
  trigonometric_mu = Magnetostatics::PermeabilityType(1.);
 
  //------------------------------------------------------------------------------
  
 
  static Magnetostatics::SolutionPotentialType
  variable_permeability_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                              return Eigen::Vector3d(
                                                     cos(PI*x(0)) * sin(PI*x(1)) * sin(PI*x(2)),
                                                     -2. * sin(PI*x(0)) * cos(PI*x(1)) * sin(PI*x(2)),
                                                     sin(PI*x(0)) * sin(PI*x(1)) * cos(PI*x(2))
                                                     );
                            };
 
  static Magnetostatics::SolutionCurlType  
  variable_permeability_sigma = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                                  return 3 * PI / (1. + x(0) + x(1) + x(2))
                                    * Eigen::Vector3d(
                                                      sin(PI*x(0)) * cos(PI*x(1)) * cos(PI*x(2)),
                                                      0.,
                                                      -cos(PI*x(0)) * cos(PI*x(1)) * sin(PI*x(2))
                                                      );
                                };
 
  static Magnetostatics::ForcingTermType
  variable_permeability_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                              double a = 3 * PI / std::pow(1. + x(0) + x(1) + x(2), 2);
                              double b = 3 * std::pow(PI, 2) / (1. + x(0) + x(1) + x(2));
                              double dy_sigmax = -a * sin(PI*x(0)) * cos(PI*x(1)) * cos(PI*x(2))
                                - b * sin(PI*x(0)) * sin(PI*x(1)) * cos(PI*x(2));
                              double dz_sigmax = -a * sin(PI*x(0)) * cos(PI*x(1)) * cos(PI*x(2))
                                - b * sin(PI*x(0)) * cos(PI*x(1)) * sin(PI*x(2));
                              ;
                              double dx_sigmaz = a * cos(PI*x(0)) * cos(PI*x(1)) * sin(PI*x(2))
                                + b * sin(PI*x(0)) * cos(PI*x(1)) * sin(PI*x(2));
 
                              double dy_sigmaz = a * cos(PI*x(0)) * cos(PI*x(1)) * sin(PI*x(2))
                                + b * cos(PI*x(0)) * sin(PI*x(1)) * sin(PI*x(2));
 
                              return Eigen::Vector3d(
                                                     dy_sigmaz,
                                                     dz_sigmax - dx_sigmaz,
                                                     -dy_sigmax
                                                     );
                            };
 
  static Magnetostatics::PermeabilityType
  variable_permeability_mu = Magnetostatics::PermeabilityType(
                        [](const Cell & T, const Eigen::Vector3d & x) -> double { return 1. + x(0) + x(1) + x(2); },
                        [](const Cell & T) -> size_t { return 1; }
                        );
 
} // end of namespace HArDCore3D
 
#endif