nedelec_project.h

Go to the documentation of this file.

/*
 * Copyright (c) 2013-2014:  G-CSC, Goethe University Frankfurt
 * Author: Dmitry Logashenko
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
/*
 * nedelec_project.h - classes for the projection of the functions based
 * on the Nedelec element to the space of the divergence-free functions.
 */
#ifndef __H__UG__PLUGINS__ELECTROMAGNETISM__NEDELEC_PROJECT__
#define __H__UG__PLUGINS__ELECTROMAGNETISM__NEDELEC_PROJECT__
 
#include "common/common.h"
#include "lib_grid/algorithms/attachment_util.h"
#include "lib_disc/common/function_group.h"
#include "lib_disc/common/groups_util.h"
#include "lib_disc/function_spaces/grid_function.h"
#include "lib_disc/spatial_disc/domain_disc.h"
#include "lib_disc/spatial_disc/domain_disc_interface.h"
#include "lib_disc/spatial_disc/elem_disc/elem_disc_interface.h"
#include "lib_disc/operator/linear_operator/transfer_interface.h"
#include "lib_disc/spatial_disc/constraints/constraint_interface.h"
#include "lib_disc/operator/linear_operator/assembled_linear_operator.h"
#include "lib_algebra/operator/interface/linear_operator_inverse.h"
#include "lib_algebra/operator/interface/pprocess.h"
 
#ifdef UG_PARALLEL
#include "lib_grid/parallelization/util/attachment_operations.hpp"
#include "lib_disc/parallelization/parallelization_util.h"
#endif
 
#include "em_material.h"
 
namespace ug{
namespace Electromagnetism{
 
 
template <typename TDomain, typename TAlgebra>
class NedelecProject : public IPProcessVector<typename TAlgebra::vector_type>
{
  typedef NedelecProject<TDomain, TAlgebra> this_type;
  
public:
  typedef TDomain domain_type;
  
  typedef typename TDomain::grid_type grid_type;
  
  typedef typename domain_type::subset_handler_type subset_handler_type;
  
  typedef TAlgebra algebra_type;
  typedef typename TAlgebra::vector_type vector_type;
  typedef typename TAlgebra::matrix_type matrix_type;
 
  typedef CPUAlgebra TPotAlgebra;
  typedef typename TPotAlgebra::vector_type pot_vector_type;
  typedef typename TPotAlgebra::matrix_type pot_matrix_type;
  typedef GridFunction<TDomain, TPotAlgebra> pot_gf_type;
 
  static const int WDim = TDomain::dim;
  
  typedef typename TDomain::position_type position_type;
 
private:
 
  typedef AInt a_vert_cond_type; // the insulators are marked with -2, not with -1!
  typedef Grid::VertexAttachmentAccessor<a_vert_cond_type> aa_vert_cond_type;
  
public:
 
  NedelecProject
  (
    SmartPtr<EMaterial<TDomain> > emMatherial, 
    SmartPtr<ApproximationSpace<TDomain> > vertApproxSpace, 
    SmartPtr<ILinearOperatorInverse<pot_vector_type> > potSolver 
  );
  
  ~NedelecProject ()
  {
    for (size_t i = 0; i < m_DVF_phi.size (); i++) delete m_DVF_phi [i];
  }
  
  void set_Dirichlet
  (
    SmartPtr<EMDirichlet<TDomain, TAlgebra> > spDirichlet 
  )
  {
    m_spDirichlet = spDirichlet;
  }
  
  void set_dampDVFs (bool val)
  {
    m_bDampDVFs = val;
  }
  
  void apply
  (
    vector_type & vec 
  );
  
  void apply
  (
    GridFunction<TDomain, TAlgebra> & u, 
    const char * fct_names 
  );
  
  void apply
  (
    GridFunction<TDomain, TAlgebra> & u, 
    const FunctionGroup & fctGrp 
  );
  
  void compute_div
  (
    SmartPtr<GridFunction<TDomain, TAlgebra> > sp_u, 
    const char * u_fct_name, 
    SmartPtr<GridFunction<TDomain, TPotAlgebra> > sp_div 
  );
  
private:
  
  inline void project_func
  (
    const SmartPtr<TDomain> & domain, 
    const SmartPtr<DoFDistribution> & edgeDD, 
    vector_type & u, 
    size_t fct, 
    const SmartPtr<DoFDistribution> & vertDD, 
    pot_gf_type & aux_rhs, 
    pot_gf_type & aux_cor 
  );
  
//  The weak divergence and gradient operators:
  
  template <typename TElem>
  void weak_div_elem_type
  (
    const TDomain & domain, 
    const DoFDistribution & edgeDD, 
    const vector_type & u, 
    size_t fct, 
    const DoFDistribution & vertDD, 
    pot_vector_type & div 
  );
  
  struct WeakDiv
  {
    this_type * m_pThis;
    const TDomain & m_domain;
    const DoFDistribution & m_edgeDD;
    const vector_type & m_u;
    size_t m_fct;
    const DoFDistribution & m_vertDD;
    pot_vector_type & m_div;
    
    WeakDiv
    (
      this_type * pThis,
      const TDomain & domain, 
      const DoFDistribution & edgeDD, 
      const vector_type & u, 
      size_t fct, 
      const DoFDistribution & vertDD, 
      pot_vector_type & div 
    )
    : m_pThis (pThis), m_domain (domain), m_edgeDD (edgeDD), m_u (u),
      m_fct (fct), m_vertDD (vertDD), m_div (div) {}
    
    template <typename TElem> void operator() (TElem &)
    {
      m_pThis->template weak_div_elem_type<TElem>
        (m_domain, m_edgeDD, m_u, m_fct, m_vertDD, m_div);
    }
  };
  
  template <typename TElem>
  void clear_div_in_conductors
  (
    const DoFDistribution & vertDD, 
    pot_vector_type & div, 
    DenseVector<VariableArray1<number> > & charge 
  );
  
  struct ClearDivInConductors
  {
    this_type * m_pThis;
    const DoFDistribution & m_vertDD;
    pot_vector_type & m_div;
    DenseVector<VariableArray1<number> > & m_charge;
    
    ClearDivInConductors
    (
      this_type * pThis,
      const DoFDistribution & vertDD, 
      pot_vector_type & div, 
      DenseVector<VariableArray1<number> > & charge 
    )
    : m_pThis (pThis), m_vertDD (vertDD), m_div (div), m_charge (charge) {}
    
    template <typename TElem> void operator() (TElem &)
    {
      m_pThis->template clear_div_in_conductors<TElem> (m_vertDD, m_div, m_charge);
    }
  };
  
  void assemble_div
  (
    const TDomain & domain, 
    const DoFDistribution & edgeDD, 
    const vector_type & u, 
    size_t fct, 
    const DoFDistribution & vertDD, 
    pot_vector_type & div, 
    DenseVector<VariableArray1<number> > & charge 
  );
  
  void distribute_cor
  (
    const TDomain & domain, 
    const DoFDistribution & edgeDD, 
    vector_type & u, 
    size_t fct, 
    const DoFDistribution & vertDD, 
    const pot_vector_type & cor 
  );
  
//  Computation of charges of the DVFs and elimination of the DVFs
  
  void alloc_DVFs
  (
    const TDomain & domain, 
    pot_gf_type & aux_rhs 
  );
  
  void compute_DVFs
  (
    pot_gf_type & aux_rhs 
  );
  
  void compute_DVF_potential_coeffs
  (
    const TDomain & domain, 
    DoFDistribution & vertDD 
  );
  
  void damp_DVFs
  (
    pot_vector_type & cor, 
    const DenseVector<VariableArray1<number> > & charge 
  );
  
  template <typename TElem>
  void mark_cond_vert_elem_type
  (
    DoFDistribution & vertDD, 
    aa_vert_cond_type & vert_base_cond 
  );
  
  struct MarkCondVert
  {
    this_type * m_pThis;
    DoFDistribution & m_vertDD;
    aa_vert_cond_type & m_vert_base_cond;
    
    MarkCondVert
    (
      this_type * pThis,
      DoFDistribution & vertDD, 
      aa_vert_cond_type & vert_base_cond 
    )
    : m_pThis (pThis), m_vertDD (vertDD), m_vert_base_cond (vert_base_cond) {}
    
    template <typename TElem> void operator() (TElem &)
    {
      m_pThis->template mark_cond_vert_elem_type<TElem> (m_vertDD, m_vert_base_cond);
    }
  };
  
  template <typename TElem>
  void integrate_div_DVF_elem_type
  (
    const TDomain & domain, 
    const DoFDistribution & vertDD, 
    const aa_vert_cond_type & vert_base_cond, 
    DenseMatrix<VariableArray2<number> > & C 
  );
  
  struct IntegrateDivDVF
  {
    this_type * m_pThis;
    const TDomain & m_domain;
    const DoFDistribution & m_vertDD;
    const aa_vert_cond_type & m_vert_base_cond;
    DenseMatrix<VariableArray2<number> > & m_C;
    
    IntegrateDivDVF
    (
      this_type * pThis,
      const TDomain & domain, 
      const DoFDistribution & vertDD, 
      const aa_vert_cond_type & vert_base_cond, 
      DenseMatrix<VariableArray2<number> > & C 
    )
    : m_pThis (pThis), m_domain (domain), m_vertDD (vertDD),
      m_vert_base_cond (vert_base_cond), m_C (C) {}
    
    template <typename TElem> void operator() (TElem &)
    {
      m_pThis->template integrate_div_DVF_elem_type<TElem>
        (m_domain, m_vertDD, m_vert_base_cond, m_C);
    }
  };
  
//  Laplace operators for the auxiliary problems
  
  template <typename TElem>
  class LocLaplaceA
  {
  public:
    typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
    static const int numCorners = ref_elem_type::numCorners;
    
    static void stiffness
    (
      GridObject * elem, 
      const position_type vCornerCoords [], 
      number loc_A [numCorners] [numCorners] 
    );
  };
  
  class AuxLaplaceLocAss : public IElemDisc<TDomain>
  {
  private:
    static const size_t _C_ = 0;
    
    static const size_t maxCorners = domain_traits<WDim>::MaxNumVerticesOfElem;
    
  public:
    AuxLaplaceLocAss
    (
      NedelecProject & master 
    );
    
    virtual void prepare_setting
    (
      const std::vector<LFEID> & vLfeID,
      bool bNonRegular
    );
  
  private:
    template <typename TElem>
    void prepare_element_loop (ReferenceObjectID roid, int si);
 
    template <typename TElem>
    void ass_JA_elem (LocalMatrix & J, const LocalVector & u, GridObject * elem, const position_type vCornerCoords []);
 
    template <typename TElem>
    void prepare_element (const LocalVector & u, GridObject * elem, ReferenceObjectID roid, const position_type vCornerCoords [])
      {};
    template <typename TElem>
    void finish_element_loop ()
      {};
    template <typename TElem>
    void ass_rhs_elem (LocalVector & d, GridObject * elem, const position_type vCornerCoords [])
      {};
    template <typename TElem>
    void ass_JM_elem (LocalMatrix & J, const LocalVector & u, GridObject * elem, const position_type vCornerCoords [])
      {
        UG_THROW ("NedelecProject: Attempt to use nonstationary discretization for the potential.");
      };
    template <typename TElem>
    void ass_dA_elem (LocalVector & d, const LocalVector & u, GridObject * elem, const position_type vCornerCoords [])
      {
        UG_THROW ("NedelecProject: Attempt to assemble non-linear system for the potential.");
      };
    template <typename TElem>
    void ass_dM_elem (LocalVector & d, const LocalVector & u, GridObject * elem, const position_type vCornerCoords [])
      {
        UG_THROW ("NedelecProject: Attempt to use nonstationary discretization for the potential.");
      };
  
    void register_all_loc_discr_funcs();
 
    struct RegisterLocalDiscr {
        RegisterLocalDiscr (AuxLaplaceLocAss* pThis) : m_pThis (pThis) {}
        AuxLaplaceLocAss* m_pThis;
        template< typename TElem > void operator () (TElem&)
        {m_pThis->register_loc_discr_func<TElem> ();}
    };
 
    template <typename TElem>
    void register_loc_discr_func();
 
  private:
    
    NedelecProject & m_master;
    
    bool m_do_assemble_here;
  };
  
  class AuxLaplaceRHS : public IDomainConstraint<TDomain, TPotAlgebra>
  {
  private:
    typedef DoFDistribution::traits<Edge>::const_iterator t_edge_iterator;
  
    static const size_t maxCorners = (size_t) element_list_traits<typename domain_traits<WDim>::DimElemList>::maxCorners;
    
    typedef DoFDistribution::traits<Vertex>::const_iterator t_vert_iterator;
  
  public:
    AuxLaplaceRHS
    (
      NedelecProject & master 
    );
    
    void set_base_conductor
    (
      int base_cond = -1 
    )
    {
      m_base_cond = base_cond;
    }
    
    void adjust_jacobian
    (
      pot_matrix_type & J,
      const pot_vector_type & u,
      ConstSmartPtr<DoFDistribution> dd,
      int type,
      number time = 0.0,
      ConstSmartPtr<VectorTimeSeries<pot_vector_type> > vSol = SPNULL,
      const number s_a0 = 1.0
    );
  
    void adjust_defect
    (
      pot_vector_type & d,
      const pot_vector_type & u,
      ConstSmartPtr<DoFDistribution> dd,
      int type,
      number time = 0.0,
      ConstSmartPtr<VectorTimeSeries<pot_vector_type> > vSol = SPNULL,
      const std::vector<number> * vScaleMass = NULL,
      const std::vector<number> * vScaleStiff = NULL
    );
  
    void adjust_solution
    (
      pot_vector_type & u,
      ConstSmartPtr<DoFDistribution> dd,
      int type,
      number time = 0.0
    );
  
    void adjust_linear
    (
      pot_matrix_type & A,
      pot_vector_type & b,
      ConstSmartPtr<DoFDistribution> dd,
      int type,
      number time = 0.0
    )
    { // Note that this function is not really used, so it needs not to be optimal.
      AuxLaplaceRHS::adjust_jacobian (A, b, dd, type, time); // the 2nd arg. is dummy: A does not depend on u
      AuxLaplaceRHS::adjust_solution (b, dd, type, time);
    }
  
    void adjust_rhs
    (
      pot_vector_type & b,
      const pot_vector_type & u,
      ConstSmartPtr<DoFDistribution> dd,
      int type,
      number time = 0.0
    )
    {
      AuxLaplaceRHS::adjust_solution (b, dd, type, time);
    }
 
    int type () const {return CT_DIRICHLET;}
  
  private:
    
    void set_zero_Dirichlet
    (
      pot_vector_type & u, 
      const DoFDistribution * dd 
    );
  
    void set_identity_Dirichlet
    (
      pot_matrix_type & A, 
      const DoFDistribution * dd 
    );
  
    template <typename TElem>
    void set_value_on_subset
    (
      int si, 
      number val, 
      pot_vector_type & u, 
      const DoFDistribution * dd 
    );
    
    struct SetValueOnSubset
    {
      AuxLaplaceRHS * m_pThis;
      int m_si;
      number m_val;
      pot_vector_type & m_u;
      const DoFDistribution * m_dd;
      
      SetValueOnSubset
      (
        AuxLaplaceRHS * pThis,
        int si, 
        number val, 
        pot_vector_type & u, 
        const DoFDistribution * dd 
      )
      : m_pThis (pThis), m_si (si), m_val (val), m_u (u), m_dd (dd) {};
      
      template <typename TElem> void operator() (TElem &)
      {
        m_pThis->template set_value_on_subset<TElem> (m_si, m_val, m_u, m_dd);
      }
    };
  
    template <typename TElem>
    void set_identity_on_subset
    (
      int si, 
      pot_matrix_type & A, 
      const DoFDistribution * dd 
    );
  
    struct SetIdentityOnSubset
    {
      AuxLaplaceRHS * m_pThis;
      int m_si;
      pot_matrix_type & m_A;
      const DoFDistribution * m_dd;
      
      SetIdentityOnSubset
      (
        AuxLaplaceRHS * pThis,
        int si, 
        pot_matrix_type & A, 
        const DoFDistribution * dd 
      )
      : m_pThis (pThis), m_si (si), m_A (A), m_dd (dd) {};
      
      template <typename TElem> void operator() (TElem &)
      {
        m_pThis->template set_identity_on_subset<TElem> (m_si, m_A, m_dd);
      }
    };
  
  private:
    
    NedelecProject & m_master;
  
    int m_base_cond;
  };
  
private:
 
  SmartPtr<EMaterial<TDomain> > m_spEmMaterial;
  
  SmartPtr<ApproximationSpace<TDomain> > m_spVertApproxSpace;
  
  bool m_bDampDVFs;
  
  SmartPtr<EMDirichlet<TDomain, TAlgebra> > m_spDirichlet;
  
  SmartPtr<AuxLaplaceLocAss> m_auxLocLaplace;
  SmartPtr<AuxLaplaceRHS> m_auxLaplaceRHS;
  SmartPtr<DomainDiscretization<TDomain, TPotAlgebra> > m_auxLaplaceAss;
  SmartPtr<AssembledLinearOperator<TPotAlgebra> > m_auxLaplaceOp;
  SmartPtr<ILinearOperatorInverse<pot_vector_type> > m_potSolver;
  std::vector<pot_gf_type *> m_DVF_phi;
  DenseMatrix<VariableArray2<number> > m_potCoeff;
};
 
} // end namespace Electromagnetism
} // end namespace ug
 
#include "nedelec_project_impl.h"
 
#endif // __H__UG__PLUGINS__ELECTROMAGNETISM__NEDELEC_PROJECT__
 
/* End of File */