levset_simple_extrapol.h

Go to the documentation of this file.

/*
 * Copyright (c) 2020:  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.
 */
 
/*
 * Extrapolation over the interface based on the simple linear extrapolation along the ray.
 */
#ifndef __H__UG__PLUGINS__LEVSET_SIMPLE_EXTRAPOL__
#define __H__UG__PLUGINS__LEVSET_SIMPLE_EXTRAPOL__
 
#include <vector>
 
// ug4 headers
#include "lib_disc/local_finite_element/local_finite_element_provider.h"
#include "lib_disc/spatial_disc/dom_disc_embb.h"
 
namespace ug {
namespace LevelSet {
 
 
template <typename TDomain, typename TAlgebra>
class LevSetGFsimpleExtrapolation
{
  typedef LevSetGFsimpleExtrapolation<TDomain, TAlgebra> this_type;
  
public:
 
  static constexpr number small_lsf_value = 1e-8;
  
 
  static bool lsf_inside
  (
    number lsf_val 
  )
  {
    return lsf_val < - small_lsf_value;
  }
  
  typedef TDomain domain_type;
  
  typedef typename TDomain::grid_type grid_type;
  
  typedef TAlgebra algebra_type;
  
  typedef CPUAlgebra lsf_algebra_type;
  
  typedef typename algebra_type::vector_type vector_type;
  
  typedef typename algebra_type::matrix_type matrix_type;
  
  typedef GridFunction<domain_type, lsf_algebra_type> ls_grid_func_type;
  
  typedef StdInjection<domain_type, lsf_algebra_type> projection_type;
  
  typedef ApproximationSpace<domain_type> approx_space_type;
  
  typedef typename domain_type::position_accessor_type position_accessor_type;
  
  static const int dim = domain_type::dim;
  
private:
 
  enum t_fs_bc
  {
    DIRICHLET_BC = 0,
    DIRICHLET_PLAIN_BC,
    NEUMANN_0_BC,
    MAX_FS_BC
  };
  
  static const size_t max_num_corners = grid_dim_traits<dim>::MaxNumVerticesOfElem;
  
public:
 
  LevSetGFsimpleExtrapolation ()
  : m_spLSF (SPNULL)
  {}
  
  void set_LSF
  (
    SmartPtr<ls_grid_func_type> spLSF 
  )
  {
    if (spLSF.valid () && (spLSF->local_finite_element_id (0) != LFEID (LFEID::LAGRANGE, dim, 1)
                || spLSF->dd()->num_fct () != 1))
      UG_THROW ("LevSetGFExtrapolation: The Level-Set Function must be scalar piecewise linear.");
    m_spLSF = spLSF;
    prepare_grid_levels ();
  }
  
  void project_LSF ();
  
  void check_dd
  (
    ConstSmartPtr<DoFDistribution> dd 
  ) const
  {
    if (m_spLSF.invalid ()) return; // there is no LSF, so no matter!
    
    if (m_vICData.size () != dd->num_fct ())
      UG_THROW ("LevSetGFExtrapolation: Mismatch of number of the functions!");
    for (size_t fct = 0; fct < dd->num_fct (); fct++)
      if (dd->local_finite_element_id (fct) != LFEID (LFEID::LAGRANGE, dim, 1))
        UG_THROW ("LevSetGFExtrapolation: The Level-Set Function works currently only with piecewise linear functions.");
  }
  
  void prepare_interface_bc
  (
    SmartPtr<approx_space_type> spApproxSpace 
  )
  {
    m_vICData.resize (spApproxSpace->num_fct ());
  }
  
  void set_Dirichlet_for
  (
    SmartPtr<approx_space_type> spApproxSpace, 
    const char* fct_name, 
    number value 
  )
  {
    prepare_interface_bc (spApproxSpace);
    size_t fct = spApproxSpace->fct_id_by_name (fct_name);
    m_vICData[fct].bc_kind = DIRICHLET_BC;
    m_vICData[fct].func = SPNULL; // we do not use the function here
    m_vICData[fct].value = value;
  }
  
  void set_Dirichlet_for
  (
    SmartPtr<approx_space_type> spApproxSpace, 
    const char* fct_name, 
    SmartPtr<CplUserData<number, dim> > func 
  )
  {
    prepare_interface_bc (spApproxSpace);
    size_t fct = spApproxSpace->fct_id_by_name (fct_name);
    m_vICData[fct].bc_kind = DIRICHLET_BC;
    m_vICData[fct].func = func;
    m_vICData[fct].value = 0; // a dummy value
  }
  
  void set_plain_Dirichlet_for
  (
    SmartPtr<approx_space_type> spApproxSpace, 
    const char* fct_name, 
    number value 
  )
  {
    prepare_interface_bc (spApproxSpace);
    size_t fct = spApproxSpace->fct_id_by_name (fct_name);
    m_vICData[fct].bc_kind = DIRICHLET_PLAIN_BC;
    m_vICData[fct].func = SPNULL; // we do not use the function here
    m_vICData[fct].value = value;
  }
  
  void set_plain_Dirichlet_for
  (
    SmartPtr<approx_space_type> spApproxSpace, 
    const char* fct_name, 
    SmartPtr<CplUserData<number, dim> > func 
  )
  {
    prepare_interface_bc (spApproxSpace);
    size_t fct = spApproxSpace->fct_id_by_name (fct_name);
    m_vICData[fct].bc_kind = DIRICHLET_PLAIN_BC;
    m_vICData[fct].func = func;
    m_vICData[fct].value = 0; // a dummy value
  }
  
  void set_Neumann0_for
  (
    SmartPtr<approx_space_type> spApproxSpace, 
    const char* fct_name 
  )
  {
    prepare_interface_bc (spApproxSpace);
    size_t fct = spApproxSpace->fct_id_by_name (fct_name);
    m_vICData[fct].bc_kind = NEUMANN_0_BC;
    m_vICData[fct].func = SPNULL; // a dummy value
    m_vICData[fct].value = 0; // a dummy value
  }
  
  void exclude_subsets
  (
    SmartPtr<approx_space_type> spApproxSpace, 
    const char* subset_names 
  )
  {
    std::string ss_names (subset_names);
    std::vector<std::string> v_ss_names;
    TokenizeTrimString (ss_names, v_ss_names);
    m_excl_ssg.set_subset_handler (spApproxSpace->subset_handler ());
    m_excl_ssg.add (v_ss_names);
  }
  
private:
 
  number boundary_value
  (
    size_t fct 
  ) const;
 
  bool corner_excluded
  (
    size_t co 
  ) const
  {
    return m_co_excluded [co];
  }
 
//----  Extrapolation ----//
public:
  
  int check_elem_lsf
  (
    size_t n_co, 
    GridObject * pElem, 
    int si, 
    int g_level, 
    bool use_hanging, 
    const MathVector<dim> vCornerCoords [], 
    number time 
  );
  
  bool corner_inside
  (
    size_t co 
  ) const
  {return lsf_inside (m_locLSF.value (0, co));}
  
  number lsf_at
  (
    size_t co 
  ) const
  {return m_locLSF.value (0, co);}
  
  template <typename TElem>
  void extrapolate_sol_by_lsf
  (
    LocalVector& locU, 
    size_t base_co 
  ) const;
  
  void extrapolate_by_lsf
  (
    size_t num_co, 
    size_t base_co, 
    number * u, 
    size_t fct 
  ) const;
  
  void extrapolate_by_lsf
  (
    size_t num_co, 
    number * u, 
    size_t fct 
  ) const;
 
  template <typename TElem>
  void clear_outer_vectors
  (
    LocalVector& locD, 
    size_t base_co 
  );
  
  template <typename TElem>
  void eliminate_extrapolated
  (
    LocalMatrix& locM, 
    size_t base_co 
  );
 
  template <typename TElem>
  void eliminate_extrapolated
  (
    LocalMatrix& locM, 
    LocalVector& locB, 
    size_t base_co 
  );
 
  template <typename TElem>
  void eliminate_extrapolated
  (
    LocalMatrix& locM, 
    LocalVector& locB 
  );
 
//----  Resetting data in the outer parts ----//
public:
 
  void clear_outer_values
  (
    vector_type & d, 
    const DoFDistribution * dd 
  ) const;
  
  void set_outer_values
  (
    vector_type & u, 
    const DoFDistribution * dd, 
    number time 
  );
  
  void set_outer_matrices
  (
    matrix_type & A, 
    const DoFDistribution * dd 
  ) const;
  
private:
 
  template <typename TElem>
  void sum_up_near_if_outer_values_for_
  (
    vector_type & u, 
    size_t fct, 
    const DoFDistribution * dd, 
    number time, 
    Grid::VertexAttachmentAccessor<ANumber> & aaBC, 
    Grid::VertexAttachmentAccessor<AUInt> & aaNumElem 
  );
  
  struct SumUpNearIfOuterValues
  {
    this_type * m_pThis;
    vector_type & m_u;
    size_t m_fct;
    const DoFDistribution * m_dd;
    number m_time;
    Grid::VertexAttachmentAccessor<ANumber> & m_aaBC;
    Grid::VertexAttachmentAccessor<AUInt> & m_aaNumElem;
    
    SumUpNearIfOuterValues
    (
      this_type * pThis, 
      vector_type & u, 
      size_t fct, 
      const DoFDistribution * dd, 
      number time, 
      Grid::VertexAttachmentAccessor<ANumber> & aaBC, 
      Grid::VertexAttachmentAccessor<AUInt> & aaNumElem 
    ) : m_pThis (pThis), m_u (u), m_fct (fct), m_dd (dd), m_time (time), m_aaBC (aaBC), m_aaNumElem (aaNumElem) {};
    
    template <typename TElem> void operator () (TElem &)
    {
      m_pThis->template sum_up_near_if_outer_values_for_<TElem> (m_u, m_fct, m_dd, m_time, m_aaBC, m_aaNumElem);
    }
  };
 
  void prepare_grid_levels ();
  
private:
 
  struct GLData
  {
    SmartPtr<ls_grid_func_type> lsf_on_gl; 
    SmartPtr<projection_type> inject; 
  };
  
  SmartPtr<ls_grid_func_type> m_spLSF; 
  std::vector<GLData> m_vGLData; 
  
  struct ICData
  {
    t_fs_bc bc_kind; 
    SmartPtr<CplUserData<number, dim> > func; 
    number value; 
    
    ICData () : bc_kind (DIRICHLET_BC), func (SPNULL), value (0) {} // dummy values
  };
  
  std::vector<ICData> m_vICData; 
  
  SubsetGroup m_excl_ssg; 
  
//  temporaty data (valid for the initialization with a particular element)
 
  LocalIndices m_indLSF; 
  LocalVector m_locLSF; 
 
  GridObject * m_pElem; 
  int m_si; 
  MathVector<dim> m_elemCenter; 
  const MathVector<dim> * m_vCornerCoords; 
  bool m_co_excluded [max_num_corners]; 
  
  number m_time; 
  
}; // class LevSetGFExtrapolation
 
 
template <typename TDomain, typename TAlgebra>
class LSGFsimpleDomainDiscretization
: public LSGFDomainDiscretization<TDomain, TAlgebra, LevSetGFsimpleExtrapolation<TDomain, TAlgebra> >
{
  typedef LSGFDomainDiscretization<TDomain, TAlgebra, LevSetGFsimpleExtrapolation<TDomain, TAlgebra> > base_type;
  
  typedef ApproximationSpace<TDomain> approx_space_type;
  
public:
 
  LSGFsimpleDomainDiscretization (SmartPtr<approx_space_type> pApproxSpace)
  :   base_type (pApproxSpace)
  {}
};
 
} // namespace LevelSet
} // end namespace ug
 
#include "levset_simple_extrapol_impl.h"
 
#endif // __H__UG__PLUGINS__LEVSET_SIMPLE_EXTRAPOL__
 
/* End of File */