level_set_pos.h

Go to the documentation of this file.

/* 
 * SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
 * SPDX-License-Identifier: EUPL-1.2
 * SPDX-FileContributor: Dmitry Logashenko
 * SPDX-FileContributor: Goethe Universität Frankfurt
 * SPDX-FileType: SOURCE
 * 
 * This file is part of d3f++.
 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
 */
 
 
/*
 * Tools for determining the position of the level set
 */
#ifndef __H__UG__PLUGINS__LEVEL_SET_POSITION_H__
#define __H__UG__PLUGINS__LEVEL_SET_POSITION_H__
 
#include <vector>
#include <iostream>
#include <limits>
#include <cmath>
 
#include "common/common.h"
 
#include "lib_grid/tools/subset_group.h"
#include "lib_grid/algorithms/space_partitioning/lg_ntree.h"
#include "lib_disc/common/function_group.h"
#include "lib_disc/common/groups_util.h"
#include "lib_disc/domain_util.h"
#include "lib_disc/local_finite_element/local_finite_element_provider.h"
#include "lib_disc/reference_element/reference_mapping_provider.h"
#include "lib_disc/reference_element/element_list_traits.h"
#include "lib_disc/function_spaces/grid_function_user_data.h"
 
namespace ug {
namespace d3f {
 
template <typename TGridFunction>
class LSPositionZ
{
public:
 
  typedef TGridFunction grid_func_type;
  
  typedef typename TGridFunction::domain_type domain_type;
  
  typedef typename TGridFunction::algebra_type algebra_type;
 
  typedef typename TGridFunction::element_type elem_type;
  
  static const int dim = TGridFunction::dim;
  
  static const size_t maxNumCorners = (size_t)
    element_list_traits<typename domain_traits<dim>::DimElemList>::maxCorners;
  
private:
  
  static number lsf_threshold () {return 1e-8;}
  
  typedef lg_ntree<dim-1, dim, elem_type> tree_type;
 
    typedef RayElemIntersectionRecord<elem_type *> intersect_rec_type;
    
  struct elem_intersect_data
  {
    elem_type * elem; 
    MathVector<dim> pnt[2]; 
    
  //  constructor
    elem_intersect_data (elem_type * e, MathVector<dim> p_min, MathVector<dim> p_max)
    : elem (e)
    {
      pnt[0] = p_min; pnt[1] = p_max;
    }
  };
    
public:
 
  LSPositionZ
  (
    SmartPtr<grid_func_type> spLSF, 
    size_t lsf_fct, 
    number ls_value = 0 
  )
  : m_spLSF (spLSF), m_lsf_fct (lsf_fct),  m_level_set_zero (ls_value),
    m_tree (* (spLSF->domain()->grid ()), spLSF->domain()->position_attachment ()),
    m_traTol (1e-10), m_traNI (1000)
  {
    m_tree.enable_warnings(false);
  }
  
  LSPositionZ
  (
    SmartPtr<grid_func_type> spLSF, 
    const char * lsf_cmp_name, 
    number ls_value = 0 
  )
  : LSPositionZ (spLSF, spLSF->fct_id_by_name (lsf_cmp_name), ls_value)
  {}
  
  LSPositionZ
  (
    SmartPtr<grid_func_type> spLSF 
  )
  : LSPositionZ (spLSF, (size_t) 0, 0)
  {}
  
  SmartPtr<grid_func_type> lsf () {return m_spLSF;}
  
  size_t lsf_fct () {return m_lsf_fct;}
  
  number level_set_zero () {return m_level_set_zero;}
  
  void reinit ()
  {
  //  find the intersected elements
    fill_cut ();
  
  //  create the tree
    m_tree.create_tree (m_vpCut.begin (), m_vpCut.end ());
  };
  
  void get_height
  (
    const MathVector<dim-1> & xy, 
    std::vector<number> & z, 
    number small_z = 1e-12 
  )
  {
    z.clear ();
    if (! get_intersected (xy)) return;
    get_z (z, small_z);
  }
  
  bool get_height_at
  (
    const MathVector<dim-1> & xy, 
    number& height 
  )
  {
    std::vector<number> z;
    get_height (xy, z);
    if (z.size () == 0) return false;
    number max_z = z[0];
    for (size_t i = 1; i < z.size (); i++)
      if (z[i] > max_z)
        max_z = z[i];
    height = max_z;
    return true;
  }
  
  number height_at_dv
  (
    const MathVector<dim-1> & xy, 
    number default_height 
  )
  {
    number v;
    if (get_height_at (xy, v))
      return v;
    return default_height;
  }
 
#ifdef UG_FOR_LUA
  number height_at
  (
    std::vector<number> xy 
  )
  {
    if (xy.size () != dim - 1)
      UG_THROW ("LSPositionZ::height_at: Only " << dim - 1 << "coordinates may be specified!");
    
    MathVector<dim-1> xy_vec;
    number v;
    for (size_t i = 0; i < dim - 1; i++) xy_vec [i] = xy [i];
    if (! get_height_at (xy_vec, v))
      return std::numeric_limits<number>::quiet_NaN ();
    return v;
  }
#endif
  
  void set_transform_tol
  (
    number tol, 
    size_t nIter 
  )
  {
    m_traTol = tol; m_traNI = nIter;
  }
  
private:
 
  inline int lsf_sign (size_t noc, number lsf []);
  
  void fill_cut ();
  
  bool get_intersected
  (
    const MathVector<dim-1> & xy 
  );
  
  void get_z
  (
    std::vector<number> & z, 
    number small_z = 1e-12 
  );
  
private:
 
  SmartPtr<grid_func_type> m_spLSF;
  
  size_t m_lsf_fct;
  
  number m_level_set_zero;
 
  tree_type m_tree;
  
  std::vector<elem_type *> m_vpCut;
  
  std::vector<intersect_rec_type> m_intersectionRecords;
  
  std::vector<elem_intersect_data> m_intersections;
  
  number m_traTol;
  size_t m_traNI;
};
 
template <typename TGridFunction>
class LSPosZData
: public StdDependentUserData<LSPosZData<TGridFunction>, number, TGridFunction::dim>
{
public:
  typedef TGridFunction gf_type;
  
  static const int wDim = gf_type::dim;
  
  typedef LSPosZData<gf_type> this_type;
 
  typedef StdDependentUserData<this_type, number, wDim> base_type;
  
  typedef typename gf_type::domain_type domain_type;
  
  typedef typename domain_type::position_accessor_type pos_acc_type;
  
private:
  typedef LSPositionZ<TGridFunction> measurer_type;
  
  typedef Attachment<size_t> num_a_type;
  
public:
 
  LSPosZData
  (
    SmartPtr<measurer_type>& sp_ls_pos_z, 
    const char * subsets 
  );
 
  ~LSPosZData ();
  
  void set_default_height
  (
    number defaultHeight 
  )
  {
    m_default_height = defaultHeight;
    m_use_ss_z_as_default = false;
  }
  
  void reinit ();
 
  void reinit_for
  (
    int dim 
  );
 
  virtual bool continuous () const {return m_elem_dim == 0;}
 
  virtual bool requires_grid_fct () const {return true;}
 
  template <int refDim>
  void eval_and_deriv
  (
    number vValue[], 
    const MathVector<wDim> vGlobIP[],
    number time,
    int si,
    GridObject* elem,
    const MathVector<wDim> vCornerCoords[],
    const MathVector<refDim> vLocIP[],
    const size_t nip,
    LocalVector* u,
    bool bDeriv,
    int s,
    std::vector<std::vector<number> > vvvDeriv[],
    const MathMatrix<refDim, wDim>* vJT = NULL
  );
  
private:
 
  template <typename TElem>
  void get_xy_for_all_
  (
    MathVector<wDim-1> * xy_pos, 
    number * z_pos 
  );
  
  struct GetXYForElements
  {
    LSPosZData * m_pThis;
    MathVector<wDim-1> * m_xy_pos;
    number * m_z_pos;
    GetXYForElements (LSPosZData * pThis, MathVector<wDim-1> * xy_pos, number * z_pos)
    : m_pThis (pThis), m_xy_pos (xy_pos), m_z_pos (z_pos) {};
    template <typename TElem> void operator() (TElem &)
    {
      m_pThis->template get_xy_for_all_<TElem> (m_xy_pos, m_z_pos);
    }
  };
  
  template <int elemDim>
  void reinit_for_ ();
  
  template <int elemDim>
  void detach_for_ ();
  
  void detach_for
  (
    int dim 
  );
  
private:
 
  SubsetGroup m_SsGrp; 
 
  SmartPtr<measurer_type> m_sp_ls_pos_z; 
  
  bool m_use_ss_z_as_default; 
  number m_default_height; 
  
  int m_elem_dim; 
  num_a_type m_a_elem_num; 
  std::vector<number> m_elem_z; 
};
 
template <typename TGridFunction>
class LSPosZCellToVrtData
: public StdDependentUserData<LSPosZCellToVrtData<TGridFunction>, number, TGridFunction::dim>
{
public:
  typedef TGridFunction gf_type;
  
  static const int wDim = gf_type::dim;
  
  typedef LSPosZCellToVrtData<gf_type> this_type;
 
  typedef StdDependentUserData<this_type, number, wDim> base_type;
  
  typedef typename gf_type::domain_type domain_type;
  
  typedef typename domain_type::grid_type grid_type;
 
  typedef typename domain_type::position_accessor_type pos_acc_type;
  
private:
  typedef LSPositionZ<TGridFunction> measurer_type;
  
  typedef Attachment<size_t> num_a_type;
  
  typedef Attachment<size_t> cnt_a_type;
  
  typedef typename Grid::AttachmentAccessor<Vertex, cnt_a_type> cnt_aa_t;
  
  typedef Attachment<number> h_a_type;
 
  typedef typename Grid::AttachmentAccessor<Vertex, h_a_type> h_aa_t;
  
public:
 
  LSPosZCellToVrtData
  (
    SmartPtr<measurer_type>& sp_ls_pos_z, 
    const char * subsets 
  );
 
  ~LSPosZCellToVrtData ();
  
  void set_default_height
  (
    number defaultHeight 
  )
  {
    m_default_height = defaultHeight;
    m_use_ss_z_as_default = false;
  }
  
  void reinit ();
 
  void reinit_for
  (
    int dim 
  );
 
  virtual bool continuous () const {return true;}
 
  virtual bool requires_grid_fct () const {return true;}
 
  template <int refDim>
  void eval_and_deriv
  (
    number vValue[], 
    const MathVector<wDim> vGlobIP[],
    number time,
    int si,
    GridObject* elem,
    const MathVector<wDim> vCornerCoords[],
    const MathVector<refDim> vLocIP[],
    const size_t nip,
    LocalVector* u,
    bool bDeriv,
    int s,
    std::vector<std::vector<number> > vvvDeriv[],
    const MathMatrix<refDim, wDim>* vJT = NULL
  );
  
private:
 
  template <typename TElem>
  void get_xy_for_all_
  (
    num_a_type & a_elem_num, 
    MathVector<wDim-1> * xy_pos, 
    number * z_pos 
  );
  
  struct GetXYForElements
  {
    LSPosZCellToVrtData * m_pThis;
    num_a_type & m_a_elem_num;
    MathVector<wDim-1> * m_xy_pos;
    number * m_z_pos;
    GetXYForElements (LSPosZCellToVrtData * pThis, num_a_type & a_elem_num, MathVector<wDim-1> * xy_pos, number * z_pos)
    : m_pThis (pThis), m_a_elem_num (a_elem_num), m_xy_pos (xy_pos), m_z_pos (z_pos) {};
    template <typename TElem> void operator() (TElem &)
    {
      m_pThis->template get_xy_for_all_<TElem> (m_a_elem_num, m_xy_pos, m_z_pos);
    }
  };
  
  template <int elemDim>
  void reinit_for_ ();
  
private:
 
  SubsetGroup m_SsGrp; 
 
  SmartPtr<measurer_type> m_sp_ls_pos_z; 
  
  bool m_use_ss_z_as_default; 
  number m_default_height; 
  
  int m_elem_dim; 
  cnt_a_type m_a_vrt_div; 
  h_a_type m_a_height; 
};
 
template <typename TLSPosZ>
class LSPositionAtPoints
{
public:
 
  typedef TLSPosZ t_ls_measurer;
  
  static const int dim = t_ls_measurer::dim;
  
public:
 
  LSPositionAtPoints
  (
    SmartPtr<t_ls_measurer> sp_ls_measurer 
  )
  : m_sp_ls_measurer (sp_ls_measurer)
  {}
  
  void clear_pnt () {m_v_points.resize (0);}
  
  void add_pnt
  (
    MathVector<dim-1>& xy, 
    number z_spec = 0, 
    bool z_given = false 
  )
  {
    m_v_points.push_back (t_pnt_data (xy, z_spec, z_given));
  }
 
#ifdef UG_FOR_LUA
  void add_pnt
  (
    std::vector<number> xyz 
  )
  {
    bool z_given = false;
    number z_spec = 0;
    
    if (xyz.size () == dim)
    {
      z_spec = xyz [dim-1];
      z_given = true;
    }
    else if (xyz.size () != dim - 1)
      UG_THROW ("LSPositionAtPoints::add_pnt: Only " << dim - 1 << "coordinates may be specified!");
    
    MathVector<dim-1> xy_vec;
    for (size_t i = 0; i < dim - 1; i++) xy_vec [i] = xyz [i];
    add_pnt (xy_vec, z_spec, z_given);
  }
#endif
  
  void load_points_from
  (
    std::istream & in, 
    bool z_given = false 
  );
 
  void load_points_from
  (
    const char * file_name, 
    bool z_given = false 
  );
 
  void compute ();
 
  size_t n_points () {return m_v_points.size ();}
  
  const MathVector<dim-1>& coord (size_t i) {return m_v_points[i].xy;}
  
  number z_at_pnt (size_t i) {return m_v_points[i].z;}
  
  number spec_z_at_pnt (size_t i) {return m_v_points[i].z_spec;}
  
  bool z_is_spec (size_t i) {return m_v_points[i].z_given;}
  
  void print_to
  (
    std::ostream & out, 
    bool print_difference = false 
  );
  
  void print_to
  (
    const char * file_name, 
    bool print_difference = false 
  );
  
  void print
  (
    bool print_difference = false 
  )
  {
    print_to (ug::GetLogAssistant().logger(), print_difference);
  }
  
  void append_to_table
  (
    const char * file_name, 
    size_t step, 
    number time 
  );
 
  void print_table_header
  (
    const char * file_name 
  );
  
private:
 
  struct t_pnt_data
  {
    MathVector<dim-1> xy; 
    number z; 
    number z_spec; 
    bool z_given; 
    
    bool valid; 
    
    t_pnt_data (MathVector<dim-1> & the_xy, number the_z_spec, bool if_z_given)
    : xy (the_xy), z (0), z_spec (the_z_spec), z_given (if_z_given), valid (false) {};
  };
 
  SmartPtr<t_ls_measurer> m_sp_ls_measurer; 
  
  std::vector<t_pnt_data> m_v_points; 
};
 
} // namespace d3f
} // end namespace ug
 
// include implementation
#include "level_set_pos_impl.h"
 
#endif // __H__UG__PLUGINS__LEVEL_SET_POSITION_H__
 
/* End of File */