nedelec_dirichlet.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.
 */
 
/*
 * Dirichlet boundary conditions for the rot-rot-discretizations based on the
 * Nedelec-type-1 elements (the Whitney-1 forms).
 */
#ifndef __H__UG__PLUGINS__ELECTROMAGNETISM__NEDELEC_DIRICHLET__
#define __H__UG__PLUGINS__ELECTROMAGNETISM__NEDELEC_DIRICHLET__
 
#include <map>
#include <vector>
 
// ug4 headers
#include "common/common.h"
 
// library-specific headers
#include "lib_grid/lg_base.h"
#include "lib_disc/common/function_group.h"
#include "lib_disc/common/groups_util.h"
#include "lib_disc/spatial_disc/domain_disc_interface.h"
#include "lib_disc/function_spaces/approximation_space.h"
#include "lib_disc/spatial_disc/user_data/const_user_data.h"
#include "lib_grid/tools/subset_handler_interface.h"
 
#include "em_material.h"
 
namespace ug{
namespace Electromagnetism{
 
 
template <typename TDomain, typename TAlgebra>
class NedelecDirichletBC
  : public EMDirichlet<TDomain, TAlgebra>
{
public:
    typedef IDomainConstraint<TDomain, TAlgebra> base_type;
    
    typedef NedelecDirichletBC<TDomain, TAlgebra> this_type;
 
    typedef TDomain domain_type;
 
    static const int dim = domain_type::dim;
 
    typedef typename domain_type::position_type position_type;
 
    typedef TAlgebra algebra_type;
 
    typedef typename algebra_type::matrix_type matrix_type;
 
    typedef typename algebra_type::vector_type vector_type;
  
private:
    typedef DoFDistribution::traits<Edge>::const_iterator t_edge_iterator;
  
public:
    NedelecDirichletBC
    (
      const char * funcNames 
    )
    : m_vDirichletFunc (TokenizeString (funcNames))
    {};
    
    void add_0 (const char * subsets);
  
    void add (MathVector<dim> & value, const char * function, const char * subsets);
    void add (std::vector<number> vValue, const char * function, const char * subsets);
  
    void add (SmartPtr<UserData<MathVector<dim>, dim> > & func, const char * function, const char * subsets);
#ifdef UG_FOR_LUA
    void add (const char * name, const char * function, const char * subsets);
#endif
 
    virtual void get_dirichlet_subsets
    (
      SubsetGroup & dirichlet_ssgrp 
    ) const;
  
private:
// Data, auxiliary types and auxiliary functions:
  
  typename domain_type::position_accessor_type m_aaPos;
  
  std::vector<std::string> m_vDirichletFunc;
  
 
  std::map<std::string, std::vector<std::string> > m_mDirichletSS;
  
  struct TConstBC
  {
    MathVector<dim> bcValue; 
    std::string fctName; 
    size_t fct; 
    std::string ssName; 
    SubsetGroup ssGrp; 
    
    TConstBC
    (
      MathVector<dim> & the_bcValue,
      std::string the_fctName,
      std::string the_ssName
    )
    : bcValue (the_bcValue), fctName (the_fctName), ssName (the_ssName)
    {}
    
    inline number operator ()
    (
      Edge * edge, 
      int si, 
      typename domain_type::position_accessor_type & aaPos, 
      number time 
    )
    {
      // Get the vector pointing along the edge from its beginning to its end:
      position_type edgeVector = aaPos [edge->vertex(1)];
      edgeVector -= aaPos [edge->vertex(0)];
      
      // Return the scalar product:
      return bcValue * edgeVector;
    }
  };
  
  std::vector<TConstBC> m_vConstBCData;
  std::map<int, std::vector<TConstBC *> > m_mConstBC;
  
  struct TUserDataBC
  {
    SmartPtr<UserData<MathVector<dim>, dim> > spBCFunc; 
    std::string fctName; 
    size_t fct; 
    std::string ssName; 
    SubsetGroup ssGrp; 
    
    TUserDataBC
    (
      SmartPtr<UserData<MathVector<dim>, dim> > the_spFunc,
      std::string the_fctName,
      std::string the_ssName
    )
    : spBCFunc (the_spFunc), fctName (the_fctName), ssName (the_ssName)
    {}
    
    inline number operator ()
    (
      Edge * edge, 
      int si, 
      typename domain_type::position_accessor_type & aaPos, 
      number time 
    )
    {
      number dofValue;
      MathVector<dim> func_value;
      
      // Get the vector pointing along the edge from its beginning to its end:
      position_type edgeVector = aaPos [edge->vertex(1)];
      edgeVector -= aaPos [edge->vertex(0)];
      
      // Evaluate the function
      (*spBCFunc) (func_value, aaPos [edge->vertex(0)], time, si);
      dofValue = func_value * edgeVector;
      
      (*spBCFunc) (func_value, aaPos [edge->vertex(1)], time, si);
      dofValue += func_value * edgeVector;
      
      return dofValue / 2;
    }
  };
 
  std::vector<TUserDataBC> m_vUserDataBCData;
  std::map<int, std::vector<TUserDataBC *> > m_mUserDataBC;
  
  std::map<int, FunctionGroup> m_mZeroBC;
  
  void add_subset
  (
    const char* f_name, 
    const char* ss_names 
  );
  
  void set_approximation_space
  (
    SmartPtr<ApproximationSpace<TDomain> > approxSpace
  )
  {
    base_type::set_approximation_space(approxSpace);
    m_aaPos = base_type::m_spApproxSpace->domain()->position_accessor();
  }
 
  void check_functions_and_subsets
  (
    FunctionGroup & functionGroup,
    SubsetGroup & subsetGroup
  );
  
  void extract_implicit ();
  
  template <typename TUserData>
  void extract_data
  (
    std::map<int, std::vector<TUserData*> >& mvUserDataBndSegment,
    std::vector<TUserData>& vUserData
  );
  void extract_data ()
  {
    if (! base_type::m_spApproxSpace.valid ())
      UG_THROW ("NedelecDirichletBC: Approximation space not set.");
    extract_data (m_mConstBC, m_vConstBCData);
    extract_data (m_mUserDataBC, m_vUserDataBCData);
    extract_implicit ();
  }
  
  template <typename TUserData>
  void adjust_solution
  (
    const std::vector<TUserData *> & vUserData,
    int si,
    vector_type& u,
    ConstSmartPtr<DoFDistribution> dd,
    number time
  );
  
  void adjust_solution_implicit
  (
    vector_type& u,
    ConstSmartPtr<DoFDistribution> dd,
    number time
  );
 
public:
// The interface:
  
  void adjust_jacobian
  (
    matrix_type & J,
    const vector_type & u,
    ConstSmartPtr<DoFDistribution> dd,
    int type,
    number time = 0.0,
    ConstSmartPtr<VectorTimeSeries<vector_type> > vSol = SPNULL,
    const number s_a0 = 1.0
  );
 
  void adjust_defect
  (
    vector_type & d,
    const vector_type & u,
    ConstSmartPtr<DoFDistribution> dd,
    int type,
    number time = 0.0,
    ConstSmartPtr<VectorTimeSeries<vector_type> > vSol = SPNULL,
    const std::vector<number> * vScaleMass = NULL,
    const std::vector<number> * vScaleStiff = NULL
  );
 
  void adjust_solution
  (
    vector_type & u,
    ConstSmartPtr<DoFDistribution> dd,
    int type,
    number time = 0.0
  );
 
  void adjust_linear
  (
    matrix_type & A,
    vector_type & b,
    ConstSmartPtr<DoFDistribution> dd,
    int type,
    number time = 0.0
  )
  {
    this_type::adjust_jacobian (A, b, dd, type, time); // the 2nd arg. is dummy: A does not depend on u
    this_type::adjust_solution (b, dd, type, time);
  };
 
  void adjust_rhs
  (
    vector_type & b,
    const vector_type & u,
    ConstSmartPtr<DoFDistribution> dd,
    int type,
    number time = 0.0
  )
  {
    this_type::adjust_solution (b, dd, type, time);
  };
 
  int type () const {return CT_DIRICHLET;}
  
};
 
} // end namespace Electromagnetism
} // end namespace ug
 
// Implementation of the functions:
#include "nedelec_dirichlet_impl.h"
 
#endif // __H__UG__PLUGINS__ELECTROMAGNETISM__NEDELEC_DIRICHLET__
 
/* End of File */