level_set_pos_impl.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: Implementation
 */
#include <limits>
 
#include "common/util/file_util.h"
#include "lib_disc/domain_traits.h"
 
#ifdef UG_PARALLEL
#include "pcl/pcl_process_communicator.h"
#include "lib_grid/parallelization/util/attachment_operations.hpp"
#endif
 
namespace ug{
namespace d3f{
 
/*----- Class template 'LSPositionZ': -----*/
 
template<typename TGridFunction>
int LSPositionZ<TGridFunction>::lsf_sign
(
  size_t noc, 
  number lsf [] 
)
{
  int pos = 0, neg = 0;
  for (size_t co = 0; co < noc; co++)
  {
    if (lsf [co] >= m_level_set_zero - lsf_threshold ()
     && lsf [co] <= m_level_set_zero + lsf_threshold ())
      return 0; // we consider such elements as intersected; this is important for computation of the gradients
    if (lsf [co] < m_level_set_zero) neg = 1;
    else pos = 1;
  }
  return pos - neg;
}
 
template <typename TGridFunction>
void LSPositionZ<TGridFunction>::fill_cut ()
{
  typedef typename TGridFunction::template dim_traits<dim>::const_iterator t_elem_iter;
 
//  reset the array
  m_vpCut.clear ();
  
//  if no LSF then nothing to do
  if (! m_spLSF.valid ()) return;
  
  LocalIndices locInd;
  LocalVector locLSF;
  
//  loop the vertices
  t_elem_iter iterEnd = m_spLSF->template end<elem_type> ();
  for (t_elem_iter iter = m_spLSF->template begin<elem_type> (); iter != iterEnd; ++iter)
  {
    elem_type * elem = *iter;
    size_t noc = elem->num_vertices ();
    number lsf [maxNumCorners];
    
  //  local values of the LSF
    m_spLSF->indices (elem, locInd);
    locLSF.resize (locInd);
    GetLocalVector (locLSF, *m_spLSF);
    for (size_t i = 0; i < noc; i++) lsf[i] = locLSF (m_lsf_fct, i);
    
  //  check the sign of the lsf
    if (lsf_sign (noc, lsf) == 0)
      m_vpCut.push_back (elem);
  }
}
 
template <typename TGridFunction>
bool LSPositionZ<TGridFunction>::get_intersected
(
  const MathVector<dim-1> & xy 
)
{
  MathVector<dim> origin, dir;
  
//  prepare the full-dimensional vectors
  VecSet (dir, 0); dir[dim-1] = 1;
  for (size_t i = 0; i < dim-1; i++) origin[i] = xy[i];
  origin[dim-1] = 0;
  
//  get the intersections
  m_intersections.clear ();
  if (! RayElementIntersections (m_intersectionRecords, m_tree, origin, dir))
    return false;
  for (size_t i = 0; i < m_intersectionRecords.size (); i++)
  {
    intersect_rec_type & r = m_intersectionRecords[i];
    m_intersections.push_back (elem_intersect_data (r.elem,
      PointOnRay(origin, dir, r.smin), PointOnRay(origin, dir, r.smax)));
  }
  return true;
}
 
template <typename TGridFunction>
void LSPositionZ<TGridFunction>::get_z
(
  std::vector<number> & z, 
  number small_z 
)
{
  LocalIndices locInd;
  LocalVector locLSF;
  std::vector<number> vShape;
    
  for (size_t i = 0; i < m_intersections.size (); i++)
  {
    elem_intersect_data & e_data = m_intersections[i];
    elem_type * elem = e_data.elem;
    
  //  get corners of element
    std::vector<MathVector<dim> > vCornerCoords;
    CollectCornerCoordinates (vCornerCoords, *elem, *m_spLSF->domain ());
 
  //  reference object id, the reference mapping and the trial space
    const ReferenceObjectID roid = elem->reference_object_id ();
    DimReferenceMapping<dim, dim> & map
      = ReferenceMappingProvider::get<dim, dim> (roid, vCornerCoords);
    const LocalShapeFunctionSet<dim>& rTrialSpace
      = LocalFiniteElementProvider::get<dim> (roid, m_spLSF->local_finite_element_id (m_lsf_fct));
 
  //  get local position of DoF
    MathVector<dim> locPos [2];
    VecSet (locPos[0], 0.5); VecSet (locPos[1], 0.5);
    map.global_to_local(locPos, e_data.pnt, 2, m_traNI, m_traTol);
 
  //  local values of the LSF
    m_spLSF->indices (elem, locInd);
    locLSF.resize (locInd);
    GetLocalVector (locLSF, *m_spLSF);
 
  //  process the points
    number lsf_value [2]; // these are the level-set values with the interface at 0 (independently on m_level_set_zero value)
    for (size_t j = 0; j < 2; j++)
    {
    //  evaluate at shapes at the point
      rTrialSpace.shapes (vShape, locPos[j]);
 
    //  compute value of the LSF at the point
      lsf_value[j] = - m_level_set_zero; // we shift the values here to 0 to the the interpolation coeffs.
      for(size_t sh = 0; sh < vShape.size(); sh++)
        lsf_value[j] += locLSF (m_lsf_fct, sh) * vShape[sh];
      
    //  "correct" the value if approx. 0
      if (fabs (lsf_value[j]) < 1e-10)
        lsf_value[j] = 0;
    }
    
  //  check if the element contains the intersection of the level set with the line
    if (lsf_value[0] * lsf_value[1] > 0) // stricktly ">" not to loose the intersections at faces
      continue;
    
  //  find the intersection (assume the linearity of the LSF along the line)
    number elem_z;
    if (lsf_value[0] == lsf_value[1]) // only possible if both are 0 (i.e. the original diff. < 2e-10 in abs. val.)
      elem_z = e_data.pnt[1][dim-1]; // in this case e_data.pnt[0] == e_data.pnt[1]
    else
      elem_z = (lsf_value[0] * e_data.pnt[1][dim-1] - lsf_value[1] * e_data.pnt[0][dim-1])
            / (lsf_value[0] - lsf_value[1]);
    
  //  check if we have already this value
    size_t k;
    for (k = 0; k < z.size (); k++)
      if (std::fabs (z[k] - elem_z) < small_z) break;
    if (k >= z.size ()) // i.e. if none found
      z.push_back (elem_z);
  }
}
 
/*----- Auxiliary class template 'LSPosZData_ElemEvalPnt': -----*/
 
template <int WDim, typename TAAPos, typename TElem>
class LSPosZData_ElemEvalPnt
{
  static const int dim = WDim;
  typedef TAAPos pos_acc_t;
  typedef TElem elem_t;
  typedef typename reference_element_traits<elem_t>::reference_element_type ref_elem_t;
  
public:
  
  static inline void get
  (
    elem_t * elem,
    pos_acc_t & aa_pos,
    MathVector<dim> & eval_pnt
  )
  {
  //  get the 'center' of the element
    eval_pnt = aa_pos [elem->vertex (0)];
    for (size_t i = 1; i < (size_t) ref_elem_t::numCorners; i++)
      eval_pnt += aa_pos [elem->vertex (i)];
    eval_pnt /= ref_elem_t::numCorners;
  }
};
 
template <int WDim, typename TAAPos>
class LSPosZData_ElemEvalPnt <WDim, TAAPos, RegularVertex>
{
  static const int dim = WDim;
  typedef TAAPos pos_acc_t;
  typedef RegularVertex elem_t;
  
public:
  
  static inline void get
  (
    elem_t * elem,
    pos_acc_t & aa_pos,
    MathVector<dim> & eval_pnt
  )
  {
  //  return the coordinates of the vertex itself
    eval_pnt = aa_pos [elem];
  }
};
 
/*----- Class template 'LSPosZData': -----*/
 
template <typename TGridFunction>
LSPosZData<TGridFunction>::LSPosZData
(
  SmartPtr<measurer_type>& sp_ls_pos_z, 
  const char * subsetNames 
)
: m_sp_ls_pos_z (sp_ls_pos_z), m_use_ss_z_as_default (true), m_elem_dim (-1)
{
  if (! sp_ls_pos_z.valid ())
    UG_THROW ("LSPosZData: The measurer not specified.");
    
  SmartPtr<gf_type> spLSF (m_sp_ls_pos_z->lsf ());
  m_SsGrp.set_subset_handler (spLSF->domain()->subset_handler ());
  m_SsGrp.add (TokenizeString (subsetNames));
}
 
template <typename TGridFunction>
LSPosZData<TGridFunction>::~LSPosZData ()
{
  detach_for (m_elem_dim);
}
 
template <typename TGridFunction>
template <typename TElem>
void LSPosZData<TGridFunction>::get_xy_for_all_
(
  MathVector<wDim-1> * xy_pos, 
  number * z_pos 
)
{
  typedef typename gf_type::template traits<TElem>::const_iterator elem_iter_t;
  typedef typename geometry_traits<TElem>::grid_base_object base_elem_t;
  typedef typename Grid::AttachmentAccessor<base_elem_t, num_a_type> num_aa_t;
  
  SmartPtr<gf_type> spLSF (m_sp_ls_pos_z->lsf ());
  pos_acc_type aa_pos = spLSF->domain()->position_accessor ();
  num_aa_t aa_elem_num (* spLSF->domain()->grid (), m_a_elem_num);
  
  for (size_t i = 0; i < m_SsGrp.size (); i++)
  {
    int si = m_SsGrp[i];
    for (elem_iter_t iter = spLSF->template begin<TElem> (si);
            iter != spLSF->template end<TElem> (si); ++iter)
    {
      TElem * elem = * iter;
      
    //  get the 'center' of the element
      MathVector<wDim> xyz;
      LSPosZData_ElemEvalPnt<wDim,pos_acc_type,TElem>::get (elem, aa_pos, xyz);
    
    //  get the position of the level set at the center
      const size_t index = aa_elem_num [(base_elem_t *) elem];
      MathVector<wDim-1> & xy = xy_pos [index];
      for (size_t k = 0; k < wDim - 1; k++) xy[k] = xyz[k];
      if (z_pos != NULL) z_pos [index] = xyz [wDim-1];
    }
  }
}
 
template <typename TGridFunction>
template <int elemDim>
void LSPosZData<TGridFunction>::reinit_for_ ()
{
  typedef typename grid_dim_traits<elemDim>::grid_base_object elem_t;
  typedef typename gf_type::template traits<elem_t>::const_iterator elem_iter_t;
  typedef typename Grid::AttachmentAccessor<elem_t, num_a_type> num_aa_t;
  typedef typename domain_traits<elemDim>::DimElemList ElemList;
  
  SmartPtr<gf_type> spLSF (m_sp_ls_pos_z->lsf ());
  bool default_used = false;
  
  if (m_elem_dim != elemDim)
  {
    detach_for (m_elem_dim);
    spLSF->domain()->grid()->template attach_to_dv<elem_t> (m_a_elem_num, (size_t) -1);
    m_elem_dim = elemDim;
  }
  
//  initialize the attachment
  num_aa_t aa_elem_num (* spLSF->domain()->grid (), m_a_elem_num);
  size_t elem_num = 0;
  for (size_t i = 0; i < m_SsGrp.size (); i++)
  {
    int si = m_SsGrp[i];
    for (elem_iter_t iter = spLSF->template begin<elem_t> (si);
            iter != spLSF->template end<elem_t> (si); ++iter)
      aa_elem_num [*iter] = elem_num++;
  }
  
# ifndef UG_PARALLEL
  {
  //  initialize the values
    m_elem_z.resize (elem_num);
    std::vector<MathVector<wDim-1> > eval_pos (elem_num);
    std::vector<number> z_pos; number * z_pos_ptr = NULL;
    if (m_use_ss_z_as_default)
    {
      z_pos.resize (elem_num);
      z_pos_ptr = & (z_pos[0]);
    }
    boost::mpl::for_each<ElemList> (GetXYForElements (this, & (eval_pos[0]), z_pos_ptr));
    for (size_t j = 0; j < eval_pos.size (); j++)
      if (! m_sp_ls_pos_z->get_height_at (eval_pos[j], m_elem_z[j]))
      {
        m_elem_z [j] = m_use_ss_z_as_default? z_pos[j] : m_default_height;
        default_used = true;
      }
  }
  
# else // i.e. ifdef UG_PARALLEL
  
  {
    pcl::ProcessCommunicator proc_comm;
    
  //  get the points where the values must be computed
    std::vector<MathVector<wDim-1> > proc_pos (elem_num);
    std::vector<number> z_pos; number * z_pos_ptr = NULL;
    if (m_use_ss_z_as_default)
    {
      z_pos.resize (elem_num);
      z_pos_ptr = & (z_pos[0]);
    }
    boost::mpl::for_each<ElemList> (GetXYForElements (this, & (proc_pos[0]), z_pos_ptr));
    
    std::vector<MathVector<wDim-1> > all_pos;
    std::vector<int> offset;
    proc_comm.allgatherv (all_pos, proc_pos, NULL, &offset);
    
  //  compute the values at all the positions
    std::vector<number> tmp_elem_z (all_pos.size ());
    for (size_t j = 0; j < tmp_elem_z.size (); j++)
      if (! m_sp_ls_pos_z->get_height_at (all_pos[j], tmp_elem_z [j]))
        tmp_elem_z [j] = - std::numeric_limits<number>::max (); // the position is not covered by the free surface in this process
    
  //  compose the common values
    std::vector<number> all_elem_z (tmp_elem_z.size ());
    proc_comm.allreduce (tmp_elem_z, all_elem_z, PCL_RO_MAX);
    
  //  extract the values from the common array
    m_elem_z.resize (elem_num);
    int my_offset = offset [proc_comm.get_local_proc_id ()];
    if (! m_use_ss_z_as_default)
    {
      for (size_t j = 0; j < elem_num; j++)
      {
        number z_val = all_elem_z [my_offset + j];
        if (z_val == - std::numeric_limits<number>::max ()) // i.e. is not set in any process
        {
          m_elem_z [j] = m_default_height;
          default_used = true;
        }
        else
          m_elem_z [j] = z_val;
      }
    }
    else
    {
      for (size_t i = 0; i < m_SsGrp.size (); i++)
      {
        int si = m_SsGrp[i];
        for (elem_iter_t iter = spLSF->template begin<elem_t> (si);
                iter != spLSF->template end<elem_t> (si); ++iter)
        {
          elem_t * vrt = * iter;
          size_t j = aa_elem_num [vrt];
          number z_val = all_elem_z [my_offset + j];
          if (z_val == - std::numeric_limits<number>::max ()) // i.e. is not set in any process
          {
            m_elem_z [j] = z_pos[j];
            default_used = true;
          }
          else
            m_elem_z [j] = z_val;
        }
      }
    }
  }
  
# endif // UG_PARALLEL
  
  if (default_used)
  {
    UG_LOG ("LSPosZData::reinit: Some parts of the level set is not under the specified subsets.");
    if (m_use_ss_z_as_default)
    {
      UG_LOG (" Position of the subsets is used instead.\n");
    }
    else
    {
      UG_LOG (" Default value " << m_default_height << " used.\n");
    }
  }
}
 
template <typename TGridFunction>
template <int elemDim>
void LSPosZData<TGridFunction>::detach_for_ ()
{
  typedef typename grid_dim_traits<elemDim>::grid_base_object elem_type;
  m_sp_ls_pos_z->lsf()->domain()->grid()->template detach_from<elem_type> (m_a_elem_num);
}
 
template <typename TGridFunction>
void LSPosZData<TGridFunction>::detach_for
(
  int dim 
)
{
  if (dim == 0) detach_for_<0> ();
  else if (dim == 1) detach_for_<1> ();
#if defined (UG_DIM_2) || defined (UG_DIM_3)
  else if (dim == 2) detach_for_<2> ();
#endif
#if defined (UG_DIM_3)
  else if (dim == 3) detach_for_<3> ();
#endif
}
 
template <typename TGridFunction>
void LSPosZData<TGridFunction>::reinit_for
(
  int dim 
)
{
  if (dim < 0) reinit (); // restart with the correct dimensionality
  else if (dim == 0) reinit_for_<0> ();
  else if (dim == 1) reinit_for_<1> ();
#if defined (UG_DIM_2) || defined (UG_DIM_3)
  else if (dim == 2) reinit_for_<2> ();
#endif
#if defined (UG_DIM_3)
  else if (dim == 3) reinit_for_<3> ();
#endif
  else
    UG_THROW ("LSPosZData::reinit_for: Illegal dimensionality " << dim << " specified.");
}
 
template <typename TGridFunction>
void LSPosZData<TGridFunction>::reinit ()
{
//  get the highest dimension of the subset elements
  int dim = m_SsGrp.get_highest_subset_dimension ();
  
# ifdef UG_PARALLEL
  {
    pcl::ProcessCommunicator proc_comm;
    dim = proc_comm.allreduce (dim, PCL_RO_MAX);
  }
# endif
 
//  reinit for the computed dimensionality
  reinit_for (dim);
}
 
template <typename TGridFunction>
template <int refDim>
void LSPosZData<TGridFunction>::eval_and_deriv
(
  number vValue[], 
  const MathVector<wDim> vGlobIP[],
  number time,
  int si,
  GridObject* grid_obj,
  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
)
{
  typedef typename gf_type::template dim_traits<refDim>::grid_base_object elem_type;
  
//  May be used only on the specified subsets!
  if (! m_SsGrp.contains (si))
    UG_THROW ("LSPosZData: Attempt to get the value for a subset that has not been specified.");
  
//  Interpolate the values
 
  const ReferenceObjectID roid = grid_obj->reference_object_id ();
  const elem_type * elem = dynamic_cast<elem_type *> (grid_obj);
 
  const LocalShapeFunctionSet<refDim> & rTrialSpace =
      LocalFiniteElementProvider::get<refDim> (roid, LFEID (LFEID::LAGRANGE, refDim, 1));
 
  if (m_elem_dim == refDim)
  {
    typedef typename grid_dim_traits<refDim>::grid_base_object base_elem_t;
    typedef typename Grid::AttachmentAccessor<base_elem_t, num_a_type> num_aa_t;
    num_aa_t aa_elem_num (* m_sp_ls_pos_z->lsf()->domain()->grid (), m_a_elem_num);
    size_t elem_no = aa_elem_num [elem];
    if (elem_no == (size_t) -1)
      UG_THROW ("LSPosZData: The object is not properly initialized.");
    number z_val = m_elem_z [elem_no];
    for (size_t ip = 0; ip < nip; ip++)
      vValue [ip] = z_val;
    if (bDeriv)
      UG_THROW ("LSPosZData: No derivative for " << m_elem_dim << "d > 0d.");
  }
  else if (m_elem_dim == 0)
  {
    Grid::VertexAttachmentAccessor<num_a_type> aa_vrt_num (* m_sp_ls_pos_z->lsf()->domain()->grid (), m_a_elem_num);
    std::vector<number> shape;
 
    for (size_t ip = 0; ip < nip; ip++)
    {
      //  evaluate at shapes at ip
      rTrialSpace.shapes (shape, vLocIP [ip]);
      UG_ASSERT (shape.size () == elem->num_vertices (), "LSPosZData: Hanging nodes are not supported!");
 
      //  compute solution at integration point
      vValue [ip] = 0;
      for(size_t sh = 0; sh < shape.size (); sh++)
      {
        const Vertex * vrt = elem->vertex (sh);
        size_t vrt_no = aa_vrt_num [vrt];
        if (vrt_no == (size_t) -1)
          UG_THROW ("LSPosZData: The object is not properly initialized.");
        vValue [ip] += m_elem_z [vrt_no] * shape [sh];
      }
    }
 
    if (bDeriv)
    {
      for (size_t ip = 0; ip < nip; ip++)
      {
        rTrialSpace.shapes(shape, vLocIP [ip]);
        for (size_t sh = 0; sh < shape.size(); sh++) vvvDeriv [ip] [0] [sh] = shape [sh];
      }
    }
  }
  else
    UG_THROW ("LSPosZData: The object is reinitialized for a wrong dimensionality " << m_elem_dim << ".");
}
 
/*----- Class template 'LSPosZCellToVrtData': -----*/
 
template <typename TGridFunction>
LSPosZCellToVrtData<TGridFunction>::LSPosZCellToVrtData
(
  SmartPtr<measurer_type>& sp_ls_pos_z, 
  const char * subsetNames 
)
: m_sp_ls_pos_z (sp_ls_pos_z), m_use_ss_z_as_default (true), m_elem_dim (-1)
{
  if (! sp_ls_pos_z.valid ())
    UG_THROW ("LSPosZCellToVrtData: The measurer not specified.");
    
  SmartPtr<gf_type> spLSF (m_sp_ls_pos_z->lsf ());
  m_SsGrp.set_subset_handler (spLSF->domain()->subset_handler ());
  m_SsGrp.add (TokenizeString (subsetNames));
  
  grid_type & grid = * spLSF->domain()->grid ();
  grid.template attach_to_dv<Vertex> (m_a_vrt_div, 0);
  grid.template attach_to_dv<Vertex> (m_a_height, 0);
}
 
template <typename TGridFunction>
LSPosZCellToVrtData<TGridFunction>::~LSPosZCellToVrtData ()
{
  grid_type & grid = * m_sp_ls_pos_z->lsf()->domain()->grid ();
  
  grid.template detach_from<Vertex> (m_a_vrt_div);
  grid.template detach_from<Vertex> (m_a_height);
}
 
template <typename TGridFunction>
template <typename TElem>
void LSPosZCellToVrtData<TGridFunction>::get_xy_for_all_
(
  num_a_type & a_elem_num, 
  MathVector<wDim-1> * xy_pos, 
  number * z_pos 
)
{
  typedef typename gf_type::template traits<TElem>::const_iterator elem_iter_t;
  typedef typename geometry_traits<TElem>::grid_base_object base_elem_t;
  typedef typename Grid::AttachmentAccessor<base_elem_t, num_a_type> num_aa_t;
  
  SmartPtr<gf_type> spLSF (m_sp_ls_pos_z->lsf ());
  pos_acc_type aa_pos = spLSF->domain()->position_accessor ();
  num_aa_t aa_elem_num (* spLSF->domain()->grid (), a_elem_num);
  
  for (size_t i = 0; i < m_SsGrp.size (); i++)
  {
    int si = m_SsGrp[i];
    for (elem_iter_t iter = spLSF->template begin<TElem> (si);
            iter != spLSF->template end<TElem> (si); ++iter)
    {
      TElem * elem = * iter;
      
    //  get the 'center' of the element
      MathVector<wDim> xyz;
      LSPosZData_ElemEvalPnt<wDim,pos_acc_type,TElem>::get (elem, aa_pos, xyz);
    
    //  get the position of the level set at the center
      const size_t index = aa_elem_num [(base_elem_t *) elem];
      MathVector<wDim-1> & xy = xy_pos [index];
      for (size_t k = 0; k < wDim - 1; k++) xy[k] = xyz[k];
      if (z_pos != NULL) z_pos [index] = xyz [wDim-1];
    }
  }
}
 
template <typename TGridFunction>
template <int elemDim>
void LSPosZCellToVrtData<TGridFunction>::reinit_for_ ()
{
  typedef typename grid_dim_traits<elemDim>::grid_base_object elem_t;
  typedef typename gf_type::template traits<elem_t>::const_iterator elem_iter_t;
  typedef typename gf_type::template traits<Vertex>::const_iterator vert_iter_t;
  typedef typename Grid::AttachmentAccessor<elem_t, num_a_type> num_aa_t;
  typedef typename domain_traits<elemDim>::DimElemList ElemList;
  
  SmartPtr<gf_type> spLSF (m_sp_ls_pos_z->lsf ());
  bool default_used = false;
  
  m_elem_dim = elemDim;
  
  grid_type & grid = * spLSF->domain()->grid ();
  
//  initialize the vertex attachments
  
  cnt_aa_t aa_vrt_div (grid, m_a_vrt_div);
  SetAttachmentValues (aa_vrt_div, spLSF->template begin<Vertex> (), spLSF->template end<Vertex> (), 0);
  
  h_aa_t aa_height (grid, m_a_height);
  SetAttachmentValues (aa_height, spLSF->template begin<Vertex> (), spLSF->template end<Vertex> (), 0);
  
//  initialize the element attachment
  num_a_type a_elem_num; // attachment keeping a numbering of the elements (vertices etc.) in the subsets
  grid.template attach_to_dv<elem_t> (a_elem_num, (size_t) -1);
  num_aa_t aa_elem_num (grid, a_elem_num);
  size_t elem_num = 0;
  for (size_t i = 0; i < m_SsGrp.size (); i++)
  {
    int si = m_SsGrp[i];
    for (elem_iter_t iter = spLSF->template begin<elem_t> (si);
            iter != spLSF->template end<elem_t> (si); ++iter)
      aa_elem_num [*iter] = elem_num++;
  }
  
//  array for the depths
  std::vector<number> elem_z;
  
# ifndef UG_PARALLEL
  {
  //  initialize the values
    elem_z.resize (elem_num);
    std::vector<MathVector<wDim-1> > eval_pos (elem_num);
    std::vector<number> z_pos; number * z_pos_ptr = NULL;
    if (m_use_ss_z_as_default)
    {
      z_pos.resize (elem_num);
      z_pos_ptr = & (z_pos[0]);
    }
    boost::mpl::for_each<ElemList> (GetXYForElements (this, a_elem_num, & (eval_pos[0]), z_pos_ptr));
    for (size_t j = 0; j < eval_pos.size (); j++)
      if (! m_sp_ls_pos_z->get_height_at (eval_pos[j], elem_z[j]))
      {
        elem_z [j] = m_use_ss_z_as_default? z_pos[j] : m_default_height;
        default_used = true;
      }
  }
  
# else // i.e. ifdef UG_PARALLEL
  
  {
    pcl::ProcessCommunicator proc_comm;
    
  //  get the points where the values must be computed
    std::vector<MathVector<wDim-1> > proc_pos (elem_num);
    std::vector<number> z_pos; number * z_pos_ptr = NULL;
    if (m_use_ss_z_as_default)
    {
      z_pos.resize (elem_num);
      z_pos_ptr = & (z_pos[0]);
    }
    boost::mpl::for_each<ElemList> (GetXYForElements (this, a_elem_num, & (proc_pos[0]), z_pos_ptr));
    
    std::vector<MathVector<wDim-1> > all_pos;
    std::vector<int> offset;
    proc_comm.allgatherv (all_pos, proc_pos, NULL, &offset);
    
  //  compute the values at all the positions
    std::vector<number> tmp_elem_z (all_pos.size ());
    for (size_t j = 0; j < tmp_elem_z.size (); j++)
      if (! m_sp_ls_pos_z->get_height_at (all_pos[j], tmp_elem_z [j]))
        tmp_elem_z [j] = - std::numeric_limits<number>::max (); // the position is not covered by the free surface in this process
    
  //  compose the common values
    std::vector<number> all_elem_z (tmp_elem_z.size ());
    proc_comm.allreduce (tmp_elem_z, all_elem_z, PCL_RO_MAX);
    
  //  extract the values from the common array
    elem_z.resize (elem_num);
    int my_offset = offset [proc_comm.get_local_proc_id ()];
    if (! m_use_ss_z_as_default)
    {
      for (size_t j = 0; j < elem_num; j++)
      {
        number z_val = all_elem_z [my_offset + j];
        if (z_val == - std::numeric_limits<number>::max ()) // i.e. is not set in any process
        {
          elem_z [j] = m_default_height;
          default_used = true;
        }
        else
          elem_z [j] = z_val;
      }
    }
    else
    {
      for (size_t i = 0; i < m_SsGrp.size (); i++)
      {
        int si = m_SsGrp[i];
        for (elem_iter_t iter = spLSF->template begin<elem_t> (si);
                iter != spLSF->template end<elem_t> (si); ++iter)
        {
          elem_t * vrt = * iter;
          size_t j = aa_elem_num [vrt];
          number z_val = all_elem_z [my_offset + j];
          if (z_val == - std::numeric_limits<number>::max ()) // i.e. is not set in any process
          {
            elem_z [j] = z_pos[j];
            default_used = true;
          }
          else
            elem_z [j] = z_val;
        }
      }
    }
  }
  
# endif // UG_PARALLEL
  
//  average the values in the vertices
  for (size_t i = 0; i < m_SsGrp.size (); i++)
  {
    int si = m_SsGrp[i];
    for (elem_iter_t iter = spLSF->template begin<elem_t> (si);
            iter != spLSF->template end<elem_t> (si); ++iter)
    {
      elem_t * elem = *iter;
      for (size_t j = 0; j < elem->size (); j++)
      {
        aa_height [(* elem) [j]] += elem_z [aa_elem_num [elem]];
        aa_vrt_div [(* elem) [j]] ++;
      }
    }
  }
  
# ifdef UG_PARALLEL
  AttachmentAllReduce<Vertex> (grid, m_a_vrt_div, PCL_RO_SUM);
  AttachmentAllReduce<Vertex> (grid, m_a_height, PCL_RO_SUM);
# endif
  
  for (vert_iter_t v_iter = spLSF->template begin<Vertex> ();
      v_iter != spLSF->template end<Vertex> (); ++v_iter)
  {
    Vertex * v = *v_iter;
    if (aa_vrt_div [v] != 0) // otherwise we did not visit that vertex
      aa_height [v] /= aa_vrt_div [v];
  }
  
//  detach the element attachment
  grid.template detach_from<elem_t> (a_elem_num);
  
//  warn about the default values if used
  if (default_used)
  {
    UG_LOG ("LSPosZData::reinit: Some parts of the level set is not under the specified subsets.");
    if (m_use_ss_z_as_default)
    {
      UG_LOG (" Position of the subsets is used instead.\n");
    }
    else
    {
      UG_LOG (" Default value " << m_default_height << " used.\n");
    }
  }
}
 
template <typename TGridFunction>
void LSPosZCellToVrtData<TGridFunction>::reinit_for
(
  int dim 
)
{
  if (dim < 0) reinit (); // restart with the correct dimensionality
  else if (dim == 0)
  {
    UG_THROW ("LSPosZCellToVrtData::reinit_for is not implemented for vertices (0d). Use LSPosZData instead.");
  }
  else if (dim == 1) reinit_for_<1> ();
#if defined (UG_DIM_2) || defined (UG_DIM_3)
  else if (dim == 2) reinit_for_<2> ();
#endif
#if defined (UG_DIM_3)
  else if (dim == 3) reinit_for_<3> ();
#endif
  else
    UG_THROW ("LSPosZCellToVrtData::reinit_for: Illegal dimensionality " << dim << " specified.");
}
 
template <typename TGridFunction>
void LSPosZCellToVrtData<TGridFunction>::reinit ()
{
//  get the highest dimension of the subset elements
  int dim = m_SsGrp.get_highest_subset_dimension ();
  
# ifdef UG_PARALLEL
  {
    pcl::ProcessCommunicator proc_comm;
    dim = proc_comm.allreduce (dim, PCL_RO_MAX);
  }
# endif
 
//  reinit for the computed dimensionality
  reinit_for (dim);
}
 
template <typename TGridFunction>
template <int refDim>
void LSPosZCellToVrtData<TGridFunction>::eval_and_deriv
(
  number vValue[], 
  const MathVector<wDim> vGlobIP[],
  number time,
  int si,
  GridObject* grid_obj,
  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
)
{
  typedef typename gf_type::template dim_traits<refDim>::grid_base_object elem_type;
  
//  May be used only on the specified subsets!
  if (! m_SsGrp.contains (si))
    UG_THROW ("LSPosZCellToVrtData: Attempt to get the value for a subset that has not been specified.");
  
//  Interpolate the values
 
  const ReferenceObjectID roid = grid_obj->reference_object_id ();
  const elem_type * elem = dynamic_cast<elem_type *> (grid_obj);
 
  const LocalShapeFunctionSet<refDim> & rTrialSpace =
      LocalFiniteElementProvider::get<refDim> (roid, LFEID (LFEID::LAGRANGE, refDim, 1));
 
  cnt_aa_t aa_vrt_div (* m_sp_ls_pos_z->lsf()->domain()->grid (), m_a_vrt_div);
  h_aa_t aa_vrt_h (* m_sp_ls_pos_z->lsf()->domain()->grid (), m_a_height);
  std::vector<number> shape;
 
  for (size_t ip = 0; ip < nip; ip++)
  {
    //  evaluate at shapes at ip
    rTrialSpace.shapes (shape, vLocIP [ip]);
    UG_ASSERT (shape.size () == elem->num_vertices (), "LSPosZCellToVrtData: Hanging nodes are not supported!");
 
    //  compute solution at integration point
    vValue [ip] = 0;
    for(size_t sh = 0; sh < shape.size (); sh++)
    {
      const Vertex * vrt = elem->vertex (sh);
      if (aa_vrt_div [vrt] == 0)
        UG_THROW ("LSPosZCellToVrtData: Attempt to get the position of the level set in a subset where the object is not initialized");
      number vrt_h = aa_vrt_h [vrt];
      vValue [ip] += vrt_h * shape [sh];
    }
  }
 
  if (bDeriv)
  {
    for (size_t ip = 0; ip < nip; ip++)
    {
      rTrialSpace.shapes(shape, vLocIP [ip]);
      for (size_t sh = 0; sh < shape.size(); sh++) vvvDeriv [ip] [0] [sh] = shape [sh];
    }
  }
}
 
/*----- Class template 'LSPositionAtPoints': -----*/
 
template <typename TLSPosZ>
void LSPositionAtPoints<TLSPosZ>::load_points_from
(
  std::istream & in, 
  bool z_given 
)
{
  std::string input_line;
  std::istringstream line_stream;
  
  line_stream.exceptions (std::istream::failbit | std::istream::badbit);
  
  while (! in.eof ())
  {
    MathVector<dim-1> xy;
    number z_spec = 0;
    
    if (in.fail ())
      UG_THROW ("LSPositionAtPoints: Could not load the points from the file!");
    std::getline (in, input_line);
    if (input_line.length () == 0)
      continue;
    try
    {
      line_stream.str (input_line);
      line_stream.clear ();
            for (size_t i = 0; i < dim-1; i++)
                line_stream >> xy[i];
      if (z_given)
        line_stream >> z_spec;
    }
    catch (std::istream::failure & e)
    {
      UG_THROW ("LSPositionAtPoints: Failed to parse line '" << input_line << "' in the input file.");
    };
    
    add_pnt (xy, z_spec, z_given);
  }
}
 
template <typename TLSPosZ>
void LSPositionAtPoints<TLSPosZ>::load_points_from
(
  const char * file_name, 
  bool z_given 
)
{
  std::string full_file_name = FindFileInStandardPaths (file_name);
  UG_COND_THROW (full_file_name.empty (), "LSPositionAtPoints: Coulnd't locate file '" << file_name << "'.");
  
  std::ifstream input (full_file_name.c_str (), std::ifstream::in);
  UG_COND_THROW (input.fail (), "LSPositionAtPoints: Cannot open data file '" << full_file_name << "' for input!");
  load_points_from (input, z_given);
}
 
template <typename TLSPosZ>
void LSPositionAtPoints<TLSPosZ>::compute ()
{
  if (m_sp_ls_measurer.invalid ())
    UG_THROW ("LSPositionAtPoints: The measurer is not specified!");
  
# ifndef UG_PARALLEL
 
  for (size_t i = 0; i < m_v_points.size (); i++)
  {
    t_pnt_data & pnt = m_v_points[i];
    pnt.valid = m_sp_ls_measurer->get_height_at (pnt.xy, pnt.z);
  }
  
# else // i.e. ifdef UG_PARALLEL
  
  {
    pcl::ProcessCommunicator proc_comm;
    std::vector<number> loc_z (m_v_points.size ()), red_z (m_v_points.size ());
    
    for (size_t i = 0; i < m_v_points.size (); i++)
      if (! m_sp_ls_measurer->get_height_at (m_v_points[i].xy, loc_z[i]))
        loc_z[i] = - std::numeric_limits<number>::max (); // the position is not covered by the free surface in this process
    
    proc_comm.allreduce (loc_z, red_z, PCL_RO_MAX);
    for (size_t i = 0; i < m_v_points.size (); i++)
    {
      t_pnt_data & pnt = m_v_points[i];
      number z_val = red_z[i];
      if (z_val != - std::numeric_limits<number>::max ()) // i.e. is not set in any process
      {
        pnt.z = z_val;
        pnt.valid = true;
      }
      else
        pnt.valid = false;
    }
  }
  
# endif // UG_PARALLEL
}
 
template <typename TLSPosZ>
void LSPositionAtPoints<TLSPosZ>::print_to
(
  std::ostream & out, 
  bool print_difference 
)
{
  if (m_sp_ls_measurer.invalid ())
    UG_THROW ("LSPositionAtPoints: The measurer is not specified!");
  
  for (size_t i = 0; i < m_v_points.size (); i++)
  {
    const t_pnt_data & pnt = m_v_points[i];
    for (size_t j = 0; j < dim - 1; j++)
      out << pnt.xy[j] << '\t';
    if (pnt.valid)
      out << pnt.z;
    else
      out << "nAn";
    if (pnt.z_given)
    {
      out << '\t' << pnt.z_spec;
      if (print_difference)
      {
        if (pnt.valid)
          out << '\t' << pnt.z - pnt.z_spec;
        else
          out << "\tnAn";
      }
    }
    out << std::endl;
  }
}
 
template <typename TLSPosZ>
void LSPositionAtPoints<TLSPosZ>::print_to
(
  const char * file_name, 
  bool print_difference 
)
{
# ifdef UG_PARALLEL
  pcl::ProcessCommunicator proc_comm;
  if (pcl::ProcRank () != proc_comm.get_proc_id (0))
    return; // only the master process should print
# endif
 
  std::ofstream output (file_name, std::ofstream::out);
  UG_COND_THROW (output.fail (), "LSPositionAtPoints: Cannot open data file '" << file_name << "' for output!");
  print_to (output, print_difference);
}
 
template <typename TLSPosZ>
void LSPositionAtPoints<TLSPosZ>::append_to_table
(
  const char * file_name, 
  size_t step, 
  number time 
)
{
# ifdef UG_PARALLEL
  pcl::ProcessCommunicator proc_comm;
  if (pcl::ProcRank () != proc_comm.get_proc_id (0))
    return; // only the master process should print
# endif
 
  std::ofstream out (file_name, std::ofstream::app);
  UG_COND_THROW (out.fail (), "LSPositionAtPoints: Cannot open data file '" << file_name << "' for output!");
  
  out.unsetf (std::ios::fixed | std::ios::scientific);
  out.precision (16);
  
  for (size_t i = 0; i < m_v_points.size (); i++)
  {
    out << step << '\t' << time << '\t';
    const t_pnt_data & pnt = m_v_points[i];
    if (pnt.valid)
      out << pnt.z;
    else
      out << "nAn";
    out << std::endl;
  }
}
 
template <typename TLSPosZ>
void LSPositionAtPoints<TLSPosZ>::print_table_header
(
  const char * file_name 
)
{
# ifdef UG_PARALLEL
  pcl::ProcessCommunicator proc_comm;
  if (pcl::ProcRank () != proc_comm.get_proc_id (0))
    return; // only the master process should print
# endif
 
  std::ofstream out (file_name, std::ofstream::out);
  UG_COND_THROW (out.fail (), "LSPositionAtPoints: Cannot open data file '" << file_name << "' for output!");
  
  out.unsetf (std::ios::fixed | std::ios::scientific);
  out.precision (8);
  
  for (size_t k = 0; k < dim - 1; k++)
  {
    out << "# ----\t----";
    for (size_t i = 0; i < m_v_points.size (); i++)
    {
      const t_pnt_data & pnt = m_v_points[i];
      out << '\t' << pnt.xy[k];
    }
    out << std::endl;
  }
  out << "# step\ttime";
  for (size_t i = 0; i < m_v_points.size (); i++)
  {
    const t_pnt_data & pnt = m_v_points[i];
    if (pnt.z_given)
      out << '\t' << pnt.z_spec;
    else
      out << '\t' << "----";
  }
  out << std::endl;
}
 
} // end namespace d3f
} // end namespace ug
 
/* End of File */