nedelec_source.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.
 */
 
#ifndef __H__UG__PLUGINS__ELECTROMAGNETISM__NEDELEC_SOURCE__
#define __H__UG__PLUGINS__ELECTROMAGNETISM__NEDELEC_SOURCE__
 
#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_disc/operator/linear_operator/transfer_interface.h"
#include "lib_algebra/operator/interface/linear_operator_inverse.h"
#include "lib_algebra/operator/debug_writer.h"
 
#ifdef UG_PARALLEL
#include "lib_grid/parallelization/util/attachment_operations.hpp"
#include "lib_disc/parallelization/parallelization_util.h"
#endif
 
namespace ug{
namespace Electromagnetism{
 
 
template <typename TDomain, typename TAlgebra>
class NedelecLoopCurrent
{
  typedef NedelecLoopCurrent<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 AChar a_edge_flag_type;
  typedef Grid::EdgeAttachmentAccessor<a_edge_flag_type> aa_edge_flag_type;
  typedef ABool a_vert_flag_type;
  typedef Grid::VertexAttachmentAccessor<a_vert_flag_type> aa_vert_flag_type;
  
public:
 
  NedelecLoopCurrent
  (
    const char * ssNames, 
    const char * posSsNames, 
    const char * cutSsNames, 
    SmartPtr<ApproximationSpace<TDomain> > vertApproxSpace, 
    SmartPtr<ILinearOperatorInverse<pot_vector_type> > potSolver 
  );
  
  void set
  (
    const char * fctNames, 
    number I 
  );
 
  void compute
  (
    SmartPtr<GridFunction<TDomain, TAlgebra> > sp_u 
  );
  
  std::string subsets () {return m_allSsNames;};
  
  SmartPtr<ITransferPostProcess<TDomain, TPotAlgebra> > zero_average ()
  {
    return m_zeroAverage;
  };
  
private:
  
  template <typename TElem>
  void mark_source_edges
  (
    const DoFDistribution & edgeDD, 
    aa_edge_flag_type & in_source 
  );
  struct MarkSourceEdges
  {
    this_type * m_pThis;
    const DoFDistribution & m_edgeDD;
    aa_edge_flag_type & m_in_source;
    
    MarkSourceEdges
    (
      this_type * pThis, 
      const DoFDistribution & edgeDD, 
      aa_edge_flag_type & in_source 
    )
    : m_pThis (pThis), m_edgeDD (edgeDD), m_in_source (in_source) {}
    
    template <typename TElem> void operator() (TElem &)
    {
      m_pThis->template mark_source_edges<TElem> (m_edgeDD, m_in_source);
    }
  };
  
  void compute_potential
  (
    pot_gf_type & pot_u 
  );
  
  void distribute_source_potential
  (
    const DoFDistribution & vertDD, 
    pot_vector_type & src_pot, 
    DoFDistribution & edgeDD, 
    size_t func, 
    number value, 
    vector_type & src_field 
  );
  
  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] 
    );
    
    static bool bnd_corners
    (
      GridObject * elem, 
      const SubsetGroup & bndGrp, 
      bool bnd_flag [numCorners] 
    );
  };
  
 
  class AuxLaplaceLocAss : public IElemDisc<TDomain>
  {
  private:
    static const size_t _C_ = 0;
    
    static const size_t maxCorners = domain_traits<WDim>::MaxNumVerticesOfElem;
    
  public:
    AuxLaplaceLocAss
    (
      NedelecLoopCurrent & 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 ass_rhs_elem (LocalVector & d, 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_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:
    
    NedelecLoopCurrent & m_master;
    
    SubsetGroup m_posSsGrp;
    SubsetGroup m_cutSsGrp;
  
    bool m_in_pos_subset;
  };
  
 
  class OutOfSource : public IDomainConstraint<TDomain, TPotAlgebra>
  {
  private:
    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;
  
  private:
  
    template <typename TElem>
    void mark_source_vertices_elem_type
    (
      const TDomain * dom 
    );
    struct MarkSourceVertices
    {
      OutOfSource * m_pThis;
      const TDomain * m_pDom;
      
      MarkSourceVertices
      (
        OutOfSource * pThis, 
        const TDomain * pDom 
      )
      : m_pThis (pThis), m_pDom (pDom) {}
      
      template <typename TElem> void operator() (TElem &)
      {
        m_pThis->template mark_source_vertices_elem_type<TElem> (m_pDom);
      }
    };
    void mark_source_vertices
    (
      const TDomain * dom 
    );
  
    void adjust_matrix
    (
      const DoFDistribution & vertDD, 
      pot_matrix_type & A 
    );
    
    void adjust_vector
    (
      const DoFDistribution & vertDD, 
      pot_vector_type & u 
    );
    
  public:
    void set_zero_average
    (
      const DoFDistribution & vertDD, 
      pot_vector_type & u 
    );
    
  public:
    OutOfSource
    (
      NedelecLoopCurrent & master 
    );
    
    void init
    (
      domain_type * dom, 
      a_vert_flag_type & a_in_source 
    )
    {
      MultiGrid * mg = dom->grid().get ();
      m_in_source.access (*mg, a_in_source);
      mark_source_vertices (dom);
#     ifdef UG_PARALLEL
      AttachmentAllReduce<Vertex> (*mg, a_in_source, PCL_RO_LOR);
#     endif
    }
    
    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
    )
    {
      adjust_matrix (* dd.get (), J);
    }
  
    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
    )
    {
      adjust_vector (* dd.get (), d);
    }
  
    void adjust_solution
    (
      pot_vector_type & u,
      ConstSmartPtr<DoFDistribution> dd,
      int type,
      number time = 0.0
    )
    {
      adjust_vector (* dd.get (), u);
    }
  
    void adjust_linear
    (
      pot_matrix_type & A,
      pot_vector_type & b,
      ConstSmartPtr<DoFDistribution> dd,
      int type,
      number time = 0.0
    )
    {
      adjust_matrix (* dd.get (), A);
      adjust_vector (* dd.get (), b);
    }
  
    void adjust_rhs
    (
      pot_vector_type & b,
      const pot_vector_type & u,
      ConstSmartPtr<DoFDistribution> dd,
      int type,
      number time = 0.0
    )
    {
      adjust_vector (* dd.get (), b);
    }
 
    int type () const {return CT_DIRICHLET;}
  
  private:
    
    NedelecLoopCurrent & m_master;
    
    aa_vert_flag_type m_in_source;
  };
  
 
  class ZeroAverage : public ITransferPostProcess<TDomain, TPotAlgebra>
  {
  public:
  
    ZeroAverage (SmartPtr<OutOfSource> spOoS) : m_spOoS (spOoS) {};
  
    virtual void post_process (SmartPtr<pot_gf_type> sp_gf)
    {
      m_spOoS->set_zero_average (* sp_gf->dd (), * sp_gf);
    }
    
  private:
  
    SmartPtr<OutOfSource> m_spOoS;
  };
  
  template <typename TElem>
  void get_flux_of_pot
  (
    const TDomain & domain, 
    const pot_vector_type & pot, 
    const DoFDistribution & vertDD, 
    number & flux 
  );
  struct GetFluxOfPotential
  {
    this_type * m_pThis;
    const domain_type & m_domain;
    const pot_vector_type & m_pot;
    const DoFDistribution & m_vertDD;
    number & m_flux;
    
    GetFluxOfPotential
    (
      this_type * pThis, 
      const TDomain & domain, 
      const pot_vector_type & pot, 
      const DoFDistribution & vertDD, 
      number & flux 
    )
    : m_pThis (pThis), m_domain (domain), m_pot (pot), m_vertDD (vertDD), m_flux (flux)
    {
      m_flux = 0;
    }
    
    template <typename TElem> void operator() (TElem &)
    {
      m_pThis->template get_flux_of_pot<TElem> (m_domain, m_pot, m_vertDD, m_flux);
    }
  };
  
private:
  
  struct TSrcData
  {
    std::string fctNames; 
    number I; 
    
    TSrcData (const char * the_fctNames, number the_I)
      : fctNames (the_fctNames), I (the_I) {};
  };
  
  std::vector<TSrcData> m_vSrcData;
 
  std::string m_allSsNames;
  std::string m_posSsNames;
  std::string m_cutSsNames;
  
  SubsetGroup m_allSsGrp;
  SubsetGroup m_posSsGrp;
  SubsetGroup m_cutSsGrp;
  
  SmartPtr<ApproximationSpace<TDomain> > m_spVertApproxSpace;
  
  SmartPtr<AuxLaplaceLocAss> m_auxLocLaplace;
  SmartPtr<OutOfSource> m_outOfSource;
  SmartPtr<ZeroAverage> m_zeroAverage;
 
  SmartPtr<DomainDiscretization<TDomain, TPotAlgebra> > m_auxLaplaceAss;
  SmartPtr<AssembledLinearOperator<TPotAlgebra> > m_auxLaplaceOp;
  SmartPtr<ILinearOperatorInverse<pot_vector_type> > m_potSolver;
};
 
} // end namespace Electromagnetism
} // end namespace ug
 
#include "nedelec_source_impl.h"
 
#endif // __H__UG__PLUGINS__ELECTROMAGNETISM__NEDELEC_PROJECT__
 
/* End of File */