obstacle_constraint_interface.h

Go to the documentation of this file.

/*
 * Copyright (c) 2013-2015:  G-CSC, Goethe University Frankfurt
 * Author: Raphael Prohl
 * 
 * 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__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__OBSTACLE_CONSTRAINT_INTERFACE__
#define __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__OBSTACLE_CONSTRAINT_INTERFACE__
 
#include "lib_disc/common/multi_index.h"
#include "lib_disc/function_spaces/grid_function.h"
#include "lib_disc/spatial_disc/user_data/const_user_data.h"
 
#include "lib_disc/spatial_disc/constraints/constraint_interface.h"
 
using namespace std;
 
namespace ug{
 
 
template <typename TDomain, typename TAlgebra>
class IObstacleConstraint:
  public IDomainConstraint<TDomain, TAlgebra>
  //  TODO: think about, restructuring the IConstraint-Interface in order to distinguish between
  //  a constraint having an effect on the assembling process and a constraint, which is considered in the
  //  solver (e.g. by means of 'adjust_restriction') only!
{
  public:
    typedef IDomainConstraint<TDomain, TAlgebra> base_type;
 
    typedef TAlgebra algebra_type;
 
    typedef typename algebra_type::matrix_type matrix_type;
 
    typedef typename algebra_type::vector_type vector_type;
 
    typedef typename vector_type::value_type value_type;
 
    typedef TDomain domain_type;
 
    static const int dim = domain_type::dim;
 
    typedef typename domain_type::position_type position_type;
 
  public:
    IObstacleConstraint(const GridFunction<TDomain, TAlgebra>& u)
    {
      clear();
 
      m_spDD = u.dof_distribution();
      m_spDomain = u.domain();
    };
 
    IObstacleConstraint(){
      //clear();
      UG_THROW("In 'IObstacleConstraint()': A constructor with a GridFunction as parameter"
          "is needed here!");
    };
 
  #ifdef UG_FOR_LUA
    void add(const char* name, const char* function);
    void add(const char* name, const char* function, const char* subsets);
  #endif
 
    void add(SmartPtr<UserData<number, dim, bool> > func, const char* function);
    void add(SmartPtr<UserData<number, dim, bool> > func, const char* function, const char* subsets);
 
    void add(SmartPtr<UserData<number, dim> > func, const char* function);
    void add(SmartPtr<UserData<number, dim> > func, const char* function, const char* subsets);
 
    void add(number value, const char* function);
    void add(number value, const char* function, const char* subsets);
 
    void add(SmartPtr<UserData<MathVector<dim>, dim> > func, const char* functions);
    void add(SmartPtr<UserData<MathVector<dim>, dim> > func, const char* functions, const char* subsets);
 
 
    void init();
 
    bool is_obs_dof(const DoFIndex& dof);
 
    void reset_active_dofs(){m_vActiveDofs.resize(0);}
 
    void active_dofs(vector<DoFIndex>& vActiveDoFs)
    {vActiveDoFs = m_vActiveDofs;}
 
 
    virtual void preprocess(){};
 
    virtual void adjust_sol_and_cor(value_type& sol_i, value_type& c_i, bool& dofIsAdmissible,
        const DoFIndex& dof) = 0;
 
    virtual void adjust_defect_to_constraint(vector_type& d) = 0;
 
 
    virtual void restrict_obs_values() = 0;
 
 
    virtual ~IObstacleConstraint(){};
 
  public:
  //  Implement 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 = NULL,
               const number s_a0 = 1.0){
      UG_LOG("IObstacleConstraint::adjust_jacobian() \n");
    };
 
    void adjust_defect(vector_type& d, const vector_type& u,
                       ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0,
                           ConstSmartPtr<VectorTimeSeries<vector_type> > vSol = NULL,
               const std::vector<number>* vScaleMass = NULL,
               const std::vector<number>* vScaleStiff = NULL){
      UG_LOG("IObstacleConstraint::adjust_defect() \n");
    };
 
    void adjust_solution(vector_type& u,
                         ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0){
      UG_LOG("IObstacleConstraint::adjust_solution() \n");
    };
 
    void adjust_linear(matrix_type& A, vector_type& b,
                       ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0){
      UG_LOG("IObstacleConstraint::adjust_linear() \n");
    };
 
    void adjust_rhs(vector_type& b, const vector_type& u,
                    ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0){
      UG_LOG("IObstacleConstraint::adjust_rhs() \n");
    };
 
    virtual void adjust_prolongation(matrix_type& P,
                     ConstSmartPtr<DoFDistribution> ddFine,
                     ConstSmartPtr<DoFDistribution> ddCoarse,
                     int type,
                     number time = 0.0){
      UG_LOG("IObstacleConstraint::adjust_prolongationP() \n");
    };
 
    virtual void adjust_restriction(matrix_type& R,
                    ConstSmartPtr<DoFDistribution> ddCoarse,
                    ConstSmartPtr<DoFDistribution> ddFine,
                    int type,
                    number time = 0.0);
 
    //TODO: does another type make more sense?
    virtual int type() const {return CT_CONSTRAINTS;}
 
  private:
    void extract_data();
 
    template <typename TUserData, typename TScheduledUserData>
    void extract_data(std::map<int, std::vector<TUserData*> >& mvUserDataObsSegment,
                      std::vector<TScheduledUserData>& vUserData);
 
    void clear();
 
    void check_functions_and_subsets(FunctionGroup& functionGroup, SubsetGroup& subsetGroup,
        size_t numFct) const;
 
    void init_obstacle_dofs_with_values(number time);
 
    template <typename TUserData>
    void init_obstacle_dofs_with_values(const std::map<int, std::vector<TUserData*> >& mvUserData, number time);
 
    template <typename TBaseElem, typename TUserData>
    void init_obstacle_dofs_with_values(const std::vector<TUserData*>& vUserData, int si, number time);
 
  protected:
    map<DoFIndex, number> m_mObstacleValues;
 
    vector<int> m_vObsSubsets;
 
    vector<DoFIndex> m_vActiveDofs;
 
  private:
    ConstSmartPtr<TDomain> m_spDomain;
 
    ConstSmartPtr<DoFDistribution> m_spDD;
 
    struct NumberData
    {
      const static bool isConditional = false;
      const static size_t numFct = 1;
      typedef MathVector<1> value_type;
      NumberData(SmartPtr<UserData<number, dim> > functor_,
             std::string fctName_)
        : spFunctor(functor_), fctName(fctName_), bWholeDomain(true)
      {}
      NumberData(SmartPtr<UserData<number, dim> > functor_,
             std::string fctName_, std::string ssName_)
        : spFunctor(functor_), fctName(fctName_), ssName(ssName_),
          bWholeDomain(false)
      {}
      bool operator()(MathVector<1>& val, const MathVector<dim> x,
              number time, int si) const
      {
        (*spFunctor)(val[0], x, time, si); return true;
      }
 
      SmartPtr<UserData<number, dim> > spFunctor;
      std::string fctName;
      std::string ssName;
      size_t fct[numFct];
      SubsetGroup ssGrp;
      bool bWholeDomain;
    };
 
    struct CondNumberData
    {
      const static bool isConditional = true;
      const static size_t numFct = 1;
      typedef MathVector<1> value_type;
      CondNumberData(SmartPtr<UserData<number, dim, bool> > functor_,
              std::string fctName_)
        : spFunctor(functor_), fctName(fctName_), bWholeDomain(true)
      {}
      CondNumberData(SmartPtr<UserData<number, dim, bool> > functor_,
              std::string fctName_, std::string ssName_)
        : spFunctor(functor_), fctName(fctName_), ssName(ssName_),
          bWholeDomain(true)
      {}
      bool operator()(MathVector<1>& val, const MathVector<dim> x,
              number time, int si) const
      {
        return (*spFunctor)(val[0], x, time, si);
      }
 
      SmartPtr<UserData<number, dim, bool> > spFunctor;
      std::string fctName;
      std::string ssName;
      size_t fct[numFct];
      SubsetGroup ssGrp;
      bool bWholeDomain;
    };
 
    struct ConstNumberData
    {
      const static bool isConditional = false;
      const static size_t numFct = 1;
      typedef MathVector<1> value_type;
      ConstNumberData(number value_,
              std::string fctName_)
        : functor(value_), fctName(fctName_), bWholeDomain(true)
      {}
      ConstNumberData(number value_,
              std::string fctName_, std::string ssName_)
        : functor(value_), fctName(fctName_), ssName(ssName_),
          bWholeDomain(false)
      {}
      inline bool operator()(MathVector<1>& val, const MathVector<dim> x,
                   number time, int si) const
      {
        val[0] = functor; return true;
      }
 
      number functor;
      std::string fctName;
      std::string ssName;
      size_t fct[numFct];
      SubsetGroup ssGrp;
      bool bWholeDomain;
    };
 
    struct VectorData
    {
      const static bool isConditional = false;
      const static size_t numFct = dim;
      typedef MathVector<dim> value_type;
      VectorData(SmartPtr<UserData<MathVector<dim>, dim> > value_,
             std::string fctName_)
        : spFunctor(value_), fctName(fctName_), bWholeDomain(true)
      {}
      VectorData(SmartPtr<UserData<MathVector<dim>, dim> > value_,
             std::string fctName_, std::string ssName_)
        : spFunctor(value_), fctName(fctName_), ssName(ssName_),
          bWholeDomain(false)
      {}
      bool operator()(MathVector<dim>& val, const MathVector<dim> x,
              number time, int si) const
      {
        (*spFunctor)(val, x, time, si); return true;
      }
 
      SmartPtr<UserData<MathVector<dim>, dim> > spFunctor;
      std::string fctName;
      std::string ssName;
      size_t fct[numFct];
      SubsetGroup ssGrp;
      bool bWholeDomain;
    };
 
    std::vector<CondNumberData> m_vCondNumberData;
    std::vector<NumberData> m_vNumberData;
    std::vector<ConstNumberData> m_vConstNumberData;
 
    std::vector<VectorData> m_vVectorData;
 
    std::map<int, std::vector<NumberData*> > m_mNumberObsSegment;
 
    std::map<int, std::vector<ConstNumberData*> > m_mConstNumberObsSegment;
 
    std::map<int, std::vector<CondNumberData*> > m_mCondNumberObsSegment;
 
    std::map<int, std::vector<VectorData*> > m_mVectorObsSegment;
};
 
 
} // end namespace ug
 
// include implementation
#include "obstacle_constraint_interface_impl.h"
 
#endif /* __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__OBSTACLE_CONSTRAINT_INTERFACE__ */