level_set.h

Go to the documentation of this file.

/*
 * Copyright (c) 2011-2015:  G-CSC, Goethe University Frankfurt
 * Author: Christian Wehner
 * 
 * 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__LEVEL_SET__LEVEL_SET_H__
#define __H__UG__PLUGINS__LEVEL_SET__LEVEL_SET_H__
 
#include <vector>
 
#include "common/common.h"
#include "lib_disc/common/groups_util.h"
#include "lib_disc/spatial_disc/user_data/data_export.h"
#include "lib_disc/spatial_disc/user_data/data_import.h"
#include "lib_disc/spatial_disc/user_data/const_user_data.h"
#include "lib_disc/spatial_disc/user_data/user_data.h"
#include "lib_disc/function_spaces/approximation_space.h"
#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
#include "lib_disc/spatial_disc/constraints/constraint_interface.h"
#include <boost/function.hpp>
#ifdef UG_FOR_LUA
#include "bindings/lua/lua_user_data.h"
#endif
 
namespace ug{
namespace LevelSet{
 
template<typename TGridFunction>
class FV1LevelSetDisc
{
    typedef typename TGridFunction::domain_type domain_type;
    
    typedef typename TGridFunction::algebra_type algebra_type;
 
    static const int dim = domain_type::dim;
 
  //  for debug type of scv-size attachment
  //  typedef typename Grid::VertexAttachmentAccessor<Attachment<number> > aaDiv;
 
    typedef typename domain_type::grid_type grid_type;
  
    typedef typename domain_type::position_accessor_type position_accessor_type;
    
    typedef typename Grid::VertexAttachmentAccessor<Attachment<number> > aaVol;
 
    typedef typename Grid::VertexAttachmentAccessor<Attachment<MathVector<dim> > > aaGrad;
    
    typedef typename TGridFunction::template dim_traits<dim>::grid_base_object ElemType;
 
    typedef typename TGridFunction::template traits<Edge>::const_iterator EdgeConstIterator;
 
    typedef typename TGridFunction::template traits<Vertex>::const_iterator VertexConstIterator;
    
    typedef typename TGridFunction::template dim_traits<dim>::const_iterator ElemIterator;
 
  public:
    FV1LevelSetDisc() :
          m_dt(0.0), m_time(0), m_gamma(1.0), m_delta(0.0),
          m_divFree(false),
          m_nrOfSteps(1),
          m_maxCFL(0),m_print(false),m_timestep_nr(0),m_limiter(false),
            m_inside_elements_si(2),
            m_outside_elements_si(3),
            m_onls_elements_si(4),
            m_inside_nodes_si(5),
            m_outside_nodes_si(6),
            m_onls_nodes_si(7)
        {
        set_source(0.0);
        }
 
    ~FV1LevelSetDisc() {};
 
        void set_dt(number deltaT){ UG_LOG("Set dt="<<deltaT<<"\n"); m_dt=deltaT; };
 
    void set_nr_of_steps(size_t n){m_nrOfSteps = n;};
    
      void set_reinit(size_t n){ m_gamma=0;m_delta=1;m_nrOfSteps=n; };
    void set_divfree(bool b){m_divFree=b;};
    void set_gamma(number gamma){m_gamma =gamma;}
    void set_delta(number delta){m_delta =delta;}
    void set_time(number t){m_time = t;}
    number get_time(){return m_time;};
    void set_info(bool b){m_print=b;};
    void set_limiter(bool b){m_limiter=b;};
    void set_timestep_nr(size_t n){m_timestep_nr = n;};
    
    bool advect_lsf(TGridFunction& uNew,TGridFunction& u);
    
      bool compute_error(TGridFunction& numsol);
      
    void set_velocity(SmartPtr<CplUserData<MathVector<dim>, dim> > user) {m_imVelocity=user;}
 
    void set_velocity(number vel_x)
    {
      SmartPtr<ConstUserVector<dim> > vel(new ConstUserVector<dim>());
      vel->set_entry(0, vel_x);
      if (dim>1) vel->set_entry(1, vel_x);
      if (dim>2) vel->set_entry(2, vel_x);
      set_velocity(vel);
    }
 
    void set_velocity(number vel_x, number vel_y)
    {
      SmartPtr<ConstUserVector<dim> > vel(new ConstUserVector<dim>());
      vel->set_entry(0, vel_x);
      vel->set_entry(1, vel_y);
      if (dim>2){
        UG_THROW("ConvectionDiffusion: Setting velocity vector of dimension 2"
                  " to a Discretization for world dim " << dim);
      }
      set_velocity(vel);
    }
 
    void set_velocity(number vel_x, number vel_y, number vel_z)
    {
      SmartPtr<ConstUserVector<dim> > vel(new ConstUserVector<dim>());
      vel->set_entry(0, vel_x);
      vel->set_entry(1, vel_y);
      vel->set_entry(2, vel_z);
      if (dim<3){
        UG_THROW("ConvectionDiffusion: Setting velocity vector of dimension 3"
                          " to a Discretization for world dim " << dim);
      }
      set_velocity(vel);
    }
 
#ifdef UG_FOR_LUA
    void set_velocity(const char* fctName)
    {
      set_velocity(LuaUserDataFactory<MathVector<dim>,dim>::create(fctName));
    }
#endif
 
    void set_source(SmartPtr<CplUserData<number,dim> > user){m_imSource = user;};
    void set_source(number val){set_source(make_sp(new ConstUserNumber<dim>(val)));}
#ifdef UG_FOR_LUA
    void set_source(const char* fctName)
    {
      set_source(LuaUserDataFactory<number,dim>::create(fctName));
    }
#endif
 
    void set_dirichlet_data(SmartPtr<CplUserData<number,dim> > d){m_imDirichlet = d;};
    void set_dirichlet_data(number val){set_dirichlet_data(make_sp(new ConstUserNumber<dim>(val)));};
#ifdef UG_FOR_LUA
    void set_dirichlet_data(const char* fctName)
    {
      set_dirichlet_data(LuaUserDataFactory<number,dim>::create(fctName));
    }
#endif
 
    bool runtimetest (TGridFunction& u);
 
    bool compute_normal(TGridFunction& vx,TGridFunction& vy,TGridFunction& u);
    bool compute_dnormal(TGridFunction& dnormal,TGridFunction& vx,TGridFunction& vy,TGridFunction& phi,TGridFunction& u);
    bool compute_ddnormal(TGridFunction& ddnormal,TGridFunction& dnormal,TGridFunction& vx,TGridFunction& vy,TGridFunction& phi,TGridFunction& u);
    
    bool set_dirichlet_boundary(TGridFunction& uNew,const char* subsets){
      try{
        m_dirichlet_sg = uNew.subset_grp_by_name(subsets);
      }UG_CATCH_THROW("ERROR while parsing Subsets.");
 
      return true;
       };
 
     bool set_outflow_boundary(TGridFunction& uNew,const char* subsets){
      try{
        m_neumann_sg = uNew.subset_grp_by_name(subsets);
      }UG_CATCH_THROW("ERROR while parsing Subsets.");
 
      return true;
    }
 
 
     bool init_ls_subsets(TGridFunction& phi);
     void create_ls_subsets(TGridFunction& phi);
     bool update_ls_subsets(TGridFunction& phi);
 
     void set_elements_active(int sign){
          if (sign==-1) if (m_inactive_sg.contains(m_inside_elements_si)==true) m_inactive_sg.remove(m_inside_elements_si);
        if (sign==0) if (m_inactive_sg.contains(m_onls_elements_si)==true) m_inactive_sg.remove(m_onls_elements_si);
          if (sign==1) if (m_inactive_sg.contains(m_outside_elements_si)==true)  m_inactive_sg.remove(m_outside_elements_si);
      }
      void set_elements_active(int signi,int signj){
          set_elements_active(signi);
          set_elements_active(signj);
      }
      void set_elements_active(int signi,int signj,int signk){
          set_elements_active(signi,signj);
          set_elements_active(signk);
      }
      void set_elements_inactive(int sign){
          if (sign==-1) m_inactive_sg.add(m_inside_elements_si);
          if (sign==0) m_inactive_sg.add(m_onls_elements_si);
          if (sign==1) m_inactive_sg.add(m_outside_elements_si);
      }
      void set_elements_inactive(int signi,int signj){
        set_elements_inactive(signi);
          set_elements_inactive(signj);
      }
      void set_elements_inactive(int signi,int signj,int signk){
          set_elements_inactive(signi,signj);
          set_elements_inactive(signk);
      }
      void set_nodes_active(int sign){
          if (sign==-1) m_inactive_sg.remove(m_inside_nodes_si);
          if (sign==0) m_inactive_sg.remove(m_onls_nodes_si);
          if (sign==1) m_inactive_sg.remove(m_outside_nodes_si);
      }
      void set_nodes_active(int signi,int signj){
          set_nodes_active(signi);
          set_nodes_active(signj);
      }
      void set_nodes_active(int signi,int signj,int signk){
          set_nodes_active(signi,signj);
          set_nodes_active(signk);
      }
      void set_nodes_inactive(int sign){
          if (sign==-1) m_inactive_sg.add(m_inside_nodes_si);
          if (sign==0) m_inactive_sg.add(m_onls_nodes_si);
          if (sign==1) m_inactive_sg.add(m_outside_nodes_si);
      }
      void set_nodes_inactive(int signi,int signj){
          set_nodes_inactive(signi);
          set_nodes_inactive(signj);
      }
      void set_nodes_inactive(int signi,int signj,int signk){
          set_nodes_inactive(signi,signj);
          set_nodes_inactive(signk);
      }
      bool overwrite(TGridFunction&,TGridFunction&,TGridFunction&,int);
        bool overwrite(TGridFunction&,number,TGridFunction&,int);
 
   protected:
 
      bool calculate_vertex_vol(TGridFunction& u, aaVol& aaVolVolume);
 
      bool calculate_vertex_grad_vol(TGridFunction& u, aaGrad& aaGradient, aaVol& aaVolume );
 
      bool calculate_vertex_grad_vol_sign(TGridFunction&, aaGrad& ,aaVol& ,TGridFunction& ,int);
 
    template <typename TElem>
    bool assemble_element(TElem& elem, DimFV1Geometry<dim>& geo, grid_type& grid,TGridFunction& uNew,const TGridFunction& uOld,aaGrad& aaGradient, aaVol& aaVolume );
 
    bool assign_dirichlet(TGridFunction&);
    bool limit_grad(TGridFunction& uOld, aaGrad& aaGradient);
 
  private:
        number m_dt;
    number m_time;
      number m_gamma;
        number m_delta;
      bool m_divFree;
      size_t m_nrOfSteps;
    number m_maxCFL;
    bool m_print;
    size_t m_timestep_nr;
    size_t m_limiter;
    SubsetGroup m_neumann_sg;
    SubsetGroup m_dirichlet_sg;
    SubsetGroup m_inactive_sg;
    int m_inside_elements_si;
    int m_outside_elements_si;
    int m_onls_elements_si;
    int m_inside_nodes_si;
    int m_outside_nodes_si;
    int m_onls_nodes_si;
 
    SmartPtr<CplUserData<MathVector<dim>, dim> > m_imVelocity;
    SmartPtr<CplUserData<number,dim> > m_imSource;
    SmartPtr<CplUserData<number,dim> > m_imDirichlet;
};
 
} // end namespace LevelSet
} // end namespace ug
 
// include implementation
#include "level_set_impl.h"
 
#endif /* __H__UG__PLUGINS__LEVEL_SET__LEVEL_SET_H__ */