rivers_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
 */
 
 
/*
 * Implementation of the surface rivers represented as sets of 1d segments in 3d.
 */
 
#ifdef UG_PARALLEL
#include "pcl/pcl_process_communicator.h"
#endif
 
namespace ug {
namespace d3f {
 
/*----- Class 'VertexLowDimRecharge': -----*/
 
template <typename TGridFunction>
void VertexLowDimRecharge<TGridFunction>::set_relative_to
(
  const char* subsetNames, 
  number top_tolerance, 
  int top_grid_level 
)
{
  domain_type& domain = * m_lsPosZ.lsf()->domain ();
  MultiGrid& mg = * domain.grid ();
  MGSubsetHandler& sh = * domain.subset_handler ();
  
//  parse the subset names
  SubsetGroup ssGrp (domain.subset_handler ());
  ssGrp.add (TokenizeString (subsetNames));
  
//  create the tree
  m_top_tolerance = top_tolerance;
  m_spTopTracerTree = SmartPtr<top_tracer_tree_t> (new top_tracer_tree_t (mg, domain.position_attachment ()));
  
//  fill the tree
  std::vector<g_surf_elem_t*> top_faces;
  for (size_t i = 0; i < ssGrp.size (); i++)
  {
    int si = ssGrp [i];
    if (top_grid_level >= 0) // if the grid level specified
      for (g_surf_iter_t it = sh.begin<g_surf_elem_t> (si, top_grid_level);
                  it != sh.end<g_surf_elem_t> (si, top_grid_level); ++it)
        top_faces.push_back (*it);
    else
      for (int lvl = 0; lvl < (int) sh.num_levels(); lvl++)
        for (g_surf_iter_t it = sh.begin<g_surf_elem_t> (si, lvl);
                    it != sh.end<g_surf_elem_t> (si, lvl); ++it)
        {
          g_surf_elem_t * t = *it;
          if (! mg.has_children (t))
            top_faces.push_back (t);
        }
  }
  m_spTopTracerTree->create_tree (top_faces.begin (), top_faces.end ());
}
 
template <typename TGridFunction>
size_t VertexLowDimRecharge<TGridFunction>::get_top_z
(
  const MathVector<dim> & over 
)
{
  m_topTracePoints.clear ();
  
  MathVector<dim> up_dir;
  up_dir = 0;
  up_dir [dim - 1] = 1;
  RayElementIntersections (m_topIntersectionRecords, * (m_spTopTracerTree.get ()), over, up_dir, m_top_tolerance);
  
  for_each_in_vec (top_intersection_record_t& r, m_topIntersectionRecords)
  {
    m_topTracePoints.push_back (PointOnRay (over, up_dir, r.smin));
  }
  end_for;
  
  return m_topTracePoints.size ();
}
 
template <typename TGridFunction>
void VertexLowDimRecharge<TGridFunction>::simple_compute
(
  SmartPtr<grid_func_type> spRecharge 
)
{
  typedef typename domain_type::position_accessor_type position_accessor_type;
  typedef typename grid_func_type::template traits<Vertex>::const_iterator const_vrt_iter_t;
 
//  Coordinates of the points
  position_accessor_type & aaPos = spRecharge->domain()->position_accessor ();
 
//  Loop the vertices
  std::vector<DoFIndex> ind (1);
  const_vrt_iter_t vrtIterEnd = spRecharge->template end<Vertex> ();
  for (const_vrt_iter_t iter = spRecharge->template begin<Vertex> (); iter != vrtIterEnd; ++iter)
  {
    Vertex * vrt = * iter;
    spRecharge->inner_dof_indices (vrt, 0, ind);
    DoFRef (*spRecharge, ind[0]) = simple_smoothed_recharge_at (aaPos[vrt]);
  }
}
 
template <typename TGridFunction>
void VertexLowDimRecharge<TGridFunction>::compute
(
  SmartPtr<grid_func_type> spRecharge 
)
{
  typedef typename domain_type::position_accessor_type position_accessor_type;
  typedef typename grid_func_type::template traits<Vertex>::const_iterator const_vrt_iter_t;
  
//  Coordinates of the points
  position_accessor_type & aaPos = spRecharge->domain()->position_accessor ();
  
//  Reinitialize the depth measurer and compute the depths of the free surface
  compute_depths ();
  
//  Loop the vertices
  std::vector<DoFIndex> ind (1);
  const_vrt_iter_t vrtIterEnd = spRecharge->template end<Vertex> ();
  for (const_vrt_iter_t iter = spRecharge->template begin<Vertex> (); iter != vrtIterEnd; ++iter)
  {
    Vertex * vrt = * iter;
    spRecharge->inner_dof_indices (vrt, 0, ind);
    DoFRef (*spRecharge, ind[0]) = smoothed_recharge_at (aaPos[vrt]);
  }
}
 
/*----- Class 'VertexWellRecharge': -----*/
 
template <typename TGridFunction>
void VertexWellRecharge<TGridFunction>::compute_depths ()
{
//  Prepare the object for a new level-set function
  this->prepare_ls_height ();
  
//  Loop the wells
  for (size_t well_i = 0; well_i < m_wells.size (); well_i++)
  {
    well_recharge_data_t& well = m_wells [well_i];
    
    well.well_depth = well.well_depth_spec;
    
  //  get the depth of the free surface at the center
    well.fs_depth = this->get_ls_height_at (well.x);
    
  //  take into the account the actual top
    if (this->ls_depth_is_relative ())
    {
    // get the height of the "top of the domain" over this segment
      number top_depth;
      if (this->min_top_depth (well.x, top_depth) == 0)
      {
        well.well_depth = std::numeric_limits<number>::max (); // the well is not covered by the top
        continue; // skip the computation of the river depth for this segment
      }
    // correct the depth of the river for this segment
      well.well_depth += top_depth;
    }
    
  //  restrict the depth of the river
    if (this->restrict_depth () && well.fs_depth != - std::numeric_limits<number>::max ()
      && well.well_depth < well.fs_depth - this->max_depth_diff ())
      well.well_depth = well.fs_depth - this->max_depth_diff ();
  }
  
#ifdef UG_PARALLEL
//  reduce the depths
  pcl::ProcessCommunicator procComm;
  std::vector<number> fs_depth_vec (m_wells.size ());
  std::vector<number> well_depth_vec (m_wells.size ());
  for (size_t well_i = 0; well_i < m_wells.size (); well_i++)
  {
    well_recharge_data_t& well = m_wells [well_i];
    fs_depth_vec [well_i] = well.fs_depth;
    well_depth_vec [well_i] = well.well_depth;
  }
  std::vector<number> red_fs_depth_vec;
  std::vector<number> red_well_depth_vec;
  procComm.allreduce (fs_depth_vec, red_fs_depth_vec, PCL_RO_MAX);
  procComm.allreduce (well_depth_vec, red_well_depth_vec, PCL_RO_MIN);
  if (red_fs_depth_vec.size () != fs_depth_vec.size () || red_well_depth_vec.size () != well_depth_vec.size ())
    UG_THROW ("VertexWellRecharge: Failed to reduce the depth vectors.");
  for (size_t well_i = 0; well_i < m_wells.size (); well_i++)
  {
    well_recharge_data_t& well = m_wells [well_i];
    well.fs_depth = red_fs_depth_vec [well_i];
    well.well_depth = red_well_depth_vec [well_i];
  }
#endif // UG_PARALLEL
}
 
template <typename TGridFunction>
number VertexWellRecharge<TGridFunction>::smoothed_recharge_at
(
  MathVector<dim> & x 
) const
{
//  get the low-dimensional vector for x
  MathVector<dim-1> xc;
  for (size_t i = 0; i < dim-1; i++) xc[i] = x[i];
  
//  loop the wells
  number sum = 0;
  for (size_t well_i = 0; well_i < m_wells.size (); well_i++)
  {
    const well_recharge_data_t& well = m_wells [well_i];
    number dist = VecDistance (xc, well.x);
    
    if (dist > this->smooth_len ())
      continue; // the well is too far away
      
  //  get the depth of the free surface at the well
    number fs_depth = well.fs_depth;
    if (fs_depth == - std::numeric_limits<number>::max ()) // this value should be set explicitly
      continue; // the segment is not under the free surface; no contribution
    
  //  get the depth of the well
    number well_depth = well.well_depth;
    if (well_depth == std::numeric_limits<number>::max ()) // this value should be set explicitly
      continue; // the segment is not under the top of the domain; no contribution
    
  //  add the contribution
    number well_q = this->object_recharge (well_depth, fs_depth, this->m_width)
                * cone_mollifier<this_type::mol_dim> (dist, this->smooth_len());
    if (this->no_inflow () && well_q >= 0)
      continue;
    sum += well_q;
  }
  
  return sum;
}
 
template <typename TGridFunction>
number VertexWellRecharge<TGridFunction>::simple_smoothed_recharge_at
(
    MathVector<dim> & x 
) const
{
  return 0.0;
}
 
/*----- Class 'VertexRiverRecharge': -----*/
 
template <typename TGridFunction>
void VertexRiverRecharge<TGridFunction>::set_river_width
(
  number riverWidth 
)
{
  typedef Grid::traits<Edge>::const_iterator EdgeConstIterator;
  
  Grid & riverGrid = m_spRiverNetwork->grid ();
  
  EdgeConstIterator iterEnd = riverGrid.template end<Edge> ();
  for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
    m_aaRiverWidth[*iter] = riverWidth;
}
 
template <typename TGridFunction>
void VertexRiverRecharge<TGridFunction>::set_river_width
(
  number riverWidth, 
  int ssi 
)
{
  typedef Grid::traits<Edge>::const_iterator EdgeConstIterator;
  
  const GridSubsetHandler& SH = m_spRiverNetwork->subset_handler ();
  
//  Loop the segments
  EdgeConstIterator iterEnd = SH.template end<Edge> (ssi);
  for (EdgeConstIterator iter = SH.template begin<Edge> (ssi); iter != iterEnd; ++iter)
    m_aaRiverWidth[*iter] = riverWidth;
}
 
template <typename TGridFunction>
void VertexRiverRecharge<TGridFunction>::set_river_width
(
  number riverWidth, 
  const char * ss_name 
)
{
  const GridSubsetHandler& SH = m_spRiverNetwork->subset_handler ();
  std::vector<std::string> names;
  
  TokenizeTrimString (ss_name, names);
  
  for (size_t i = 0; i < names.size (); i++)
  {
    const int ssi = SH.get_subset_index (names[i].c_str());
    if (ssi < 0)
      UG_THROW ("VertexRiverRecharge::set_river_width: no subset '" << names[i] << "' found in the river grid!");
    set_river_width (riverWidth, ssi);
  }
}
 
template <typename TGridFunction>
void VertexRiverRecharge<TGridFunction>::compute_fs_depths_at_nodes ()
{
  typedef SegmentNetwork::position_accessor_t river_vrt_pos_acc_t;
  Grid & riverGrid = m_spRiverNetwork->grid ();
  river_vrt_pos_acc_t river_vrt_pos_acc = river_network()->position_accessor ();
  
  typedef Grid::traits<Edge>::const_iterator EdgeConstIterator;
  typedef Grid::traits<Vertex>::const_iterator VertConstIterator;
  
  typedef Grid::VertexAttachmentAccessor<depth_attachment_type> nodal_depth_accessor_type;
  depth_attachment_type aNodalFSDepth;
  nodal_depth_accessor_type aaNodalFSDepth;
  riverGrid.attach_to_vertices (aNodalFSDepth);
  aaNodalFSDepth.access (riverGrid, aNodalFSDepth);
  
//  Prepare the object for a new level-set function
  this->prepare_ls_height ();
  
//  Loop the river nodes
  VertConstIterator vertIterEnd = riverGrid.template end<Vertex> ();
  for (VertConstIterator vertIter = riverGrid.template begin<Vertex> (); vertIter != vertIterEnd; ++vertIter)
  {
    Vertex * vert = *vertIter;
    aaNodalFSDepth[vert] = this->get_ls_height_at (river_vrt_pos_acc [vert]);
  }
  
#ifdef UG_PARALLEL
//  reduce the depths
  pcl::ProcessCommunicator procComm;
  std::vector<number> fs_depth_vec (riverGrid.num_vertices ());
  size_t vert_i = 0;
  for (VertConstIterator vertIter = riverGrid.template begin<Vertex> (); vertIter != vertIterEnd; ++vertIter)
    fs_depth_vec [vert_i++] = aaNodalFSDepth [*vertIter];
  std::vector<number> red_fs_depth_vec;
  procComm.allreduce (fs_depth_vec, red_fs_depth_vec, PCL_RO_MAX);
  if (red_fs_depth_vec.size () != fs_depth_vec.size ())
    UG_THROW ("VertexRiverRecharge: Failed to reduce the free surface depth vector.");
  vert_i = 0;
  for (VertConstIterator vertIter = riverGrid.template begin<Vertex> (); vertIter != vertIterEnd; ++vertIter)
    aaNodalFSDepth [*vertIter] = red_fs_depth_vec [vert_i++];
#endif // UG_PARALLEL
  
//  Loop the segments
  EdgeConstIterator iterEnd = riverGrid.template end<Edge> ();
  for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
  {
    Edge * seg = *iter;
    
  //  get the depths of the free surface at the ends of the segment
    number fs_depth_1 = aaNodalFSDepth [seg->vertex (0)];
    number fs_depth_2 = aaNodalFSDepth [seg->vertex (1)];
  
  //  check the validity
    if (fs_depth_1 != - std::numeric_limits<number>::max () && fs_depth_2 != - std::numeric_limits<number>::max ())
      m_aaFSDepth[seg] = (fs_depth_1 + fs_depth_2) / 2;
    else // otherwise we skip the segment
      m_aaFSDepth[seg] = - std::numeric_limits<number>::max ();
  }
  
  riverGrid.detach_from_vertices (aNodalFSDepth);
}
 
template <typename TGridFunction>
void VertexRiverRecharge<TGridFunction>::compute_fs_depths_at_segments ()
{
  typedef SegmentNetwork::position_accessor_t river_vrt_pos_acc_t;
  typedef Grid::traits<Edge>::const_iterator EdgeConstIterator;
  
  Grid & riverGrid = m_spRiverNetwork->grid ();
  river_vrt_pos_acc_t river_vrt_pos_acc = river_network()->position_accessor ();
  
//  Prepare the object for a new level-set function
  this->prepare_ls_height ();
  
//  Loop the segments
  EdgeConstIterator iterEnd = riverGrid.template end<Edge> ();
  for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
  {
    Edge * seg = *iter;
    
  //  get the center of the segment
    MathVector<dim-1> seg_center_c;
    seg_center_c = river_vrt_pos_acc [seg->vertex (0)];
    seg_center_c += river_vrt_pos_acc [seg->vertex (1)];
    seg_center_c /= 2;
    
  //  get the depth of the free surface at the center
    number fs_depth = this->get_ls_height_at (seg_center_c);
    m_aaFSDepth[seg] = fs_depth;
  }
  
#ifdef UG_PARALLEL
//  reduce the depths
  pcl::ProcessCommunicator procComm;
  std::vector<number> fs_depth_vec (riverGrid.num_edges ());
  size_t seg_i = 0;
  for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
    fs_depth_vec [seg_i++] = m_aaFSDepth [*iter];
  std::vector<number> red_fs_depth_vec;
  procComm.allreduce (fs_depth_vec, red_fs_depth_vec, PCL_RO_MAX);
  if (red_fs_depth_vec.size () != fs_depth_vec.size ())
    UG_THROW ("VertexRiverRecharge: Failed to reduce the free surface depth vector.");
  seg_i = 0;
  for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
    m_aaFSDepth [*iter] = red_fs_depth_vec [seg_i++];
#endif // UG_PARALLEL
}
 
template <typename TGridFunction>
void VertexRiverRecharge<TGridFunction>::compute_depths ()
{
  typedef SegmentNetwork::position_accessor_t river_vrt_pos_acc_t;
  typedef SegmentNetwork::gauge_accessor_t river_depth_acc_t;
  typedef Grid::traits<Edge>::const_iterator EdgeConstIterator;
  
  Grid & riverGrid = m_spRiverNetwork->grid ();
  river_vrt_pos_acc_t river_vrt_pos_acc = river_network()->position_accessor ();
  river_depth_acc_t river_vrt_depth_acc = river_network()->gauge_accessor ();
  
//  Compute the depth of the free surface
  if (m_bNodalFSDepth)
    this->compute_fs_depths_at_nodes ();
  else
    this->compute_fs_depths_at_segments ();
  
//  Loop the segments
  EdgeConstIterator iterEnd = riverGrid.template end<Edge> ();
  for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
  {
    Edge * seg = *iter;
    
  //  get the data at the ends of the segment
  
    Vertex * vrt_0 = seg->vertex (0);
    number depth_0 = river_vrt_depth_acc[vrt_0];
    
    Vertex * vrt_1 = seg->vertex (1);
    number depth_1 = river_vrt_depth_acc[vrt_1];
    
    number seg_depth = (depth_0 + depth_1) / 2;
    
  //  take into the account the actual top
    if (this->ls_depth_is_relative ())
    {
    //  get the data at the center
      MathVector<dim-1> pos_0 = river_vrt_pos_acc[vrt_0];
      MathVector<dim-1> pos_1 = river_vrt_pos_acc[vrt_1];
      MathVector<dim-1> seg_center_c;
      seg_center_c = pos_0; seg_center_c += pos_1; seg_center_c /= 2;
      
    // get the height of the "top of the domain" over this segment
      number top_depth;
      if (this->min_top_depth (seg_center_c, top_depth) == 0)
      {
        m_aaRiverDepth[seg] = std::numeric_limits<number>::max (); // the segment is not covered by the top
        continue; // skip the computation of the river depth for this segment
      }
    // correct the depth of the river for this segment
      seg_depth += top_depth;
    }
    
  //  restrict the depth of the river
    number fs_depth = m_aaFSDepth[seg];
    if (this->restrict_depth () && fs_depth != - std::numeric_limits<number>::max ()
      && seg_depth < fs_depth - this->max_depth_diff ())
      seg_depth = fs_depth - this->max_depth_diff ();
  
  //  set the river depth for the segment
    m_aaRiverDepth[seg] = seg_depth;
  }
  
#ifdef UG_PARALLEL
//  reduce the depths
  pcl::ProcessCommunicator procComm;
  std::vector<number> river_depth_vec (riverGrid.num_edges ());
  size_t seg_i = 0;
  for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
    river_depth_vec [seg_i++] = m_aaRiverDepth [*iter];
  std::vector<number> red_river_depth_vec;
  procComm.allreduce (river_depth_vec, red_river_depth_vec, PCL_RO_MIN);
  if (red_river_depth_vec.size () != river_depth_vec.size ())
    UG_THROW ("VertexRiverRecharge: Failed to reduce the river depth vector.");
  seg_i = 0;
  for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
    m_aaRiverDepth [*iter] = red_river_depth_vec [seg_i++];
#endif // UG_PARALLEL
 
//  check the depths
  if (m_bCheckTopCovering)
    for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
    {
      Edge * seg = *iter;
      if (m_aaRiverDepth [seg] == std::numeric_limits<number>::max ())
        UG_THROW ("VertexRiverRecharge: River segment ["
          << river_vrt_pos_acc[seg->vertex(0)] << ", " << river_vrt_pos_acc[seg->vertex(1)]
            << "] is not (completely) covered by the top subset.");
    }
  if (m_bCheckFSCovering)
    for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
    {
      Edge * seg = *iter;
      if (m_aaRiverDepth [seg] != std::numeric_limits<number>::max ()
        && m_aaFSDepth [seg] == - std::numeric_limits<number>::max ())
        UG_THROW ("VertexRiverRecharge: River segment ["
          << river_vrt_pos_acc[seg->vertex(0)] << ", " << river_vrt_pos_acc[seg->vertex(1)]
            << "] is not (completely) covered by the free surface.");
    }
    
//  print the total recharge
  if (m_bPrintTotalRecharge)
    do_print_total_recharge ();
}
 
template <typename TGridFunction>
void VertexRiverRecharge<TGridFunction>::do_print_total_recharge ()
{
  typedef SegmentNetwork::position_accessor_t river_vrt_pos_acc_t;
  typedef SegmentNetwork::gauge_accessor_t river_depth_acc_t;
  typedef Grid::traits<Edge>::const_iterator EdgeConstIterator;
 
  Grid & riverGrid = m_spRiverNetwork->grid ();
  river_vrt_pos_acc_t river_vrt_pos_acc = river_network()->position_accessor ();
  number total_recharge = 0;
  
  EdgeConstIterator iterEnd = riverGrid.template end<Edge> ();
  for (EdgeConstIterator iter = riverGrid.template begin<Edge> (); iter != iterEnd; ++iter)
  {
    Edge * seg = *iter;
  
  //  get the length of the segment
    MathVector<dim-1> pos_0 = river_vrt_pos_acc[seg->vertex (0)];
    MathVector<dim-1> pos_1 = river_vrt_pos_acc[seg->vertex (1)];
    number seg_length = VecDistance (pos_0, pos_1);
  
  //  get the depth of the free surface at the center
    number fs_depth = fs_depth_for (seg);
    if (fs_depth == - std::numeric_limits<number>::max ()) // this value should be set explicitly
      continue; // the segment is not under the free surface; no contribution
  
  //  get the depth of the river
    number seg_depth = river_depth_for (seg);
    if (seg_depth == std::numeric_limits<number>::max ()) // this value should be set explicitly
      continue; // the segment is not under the top of the domain; no contribution
    
    number seg_recharge = this->object_recharge (seg_depth, fs_depth, river_width_for (seg));
    if (this->no_inflow () && seg_recharge >= 0)
      continue;
    
    total_recharge += seg_recharge * seg_length;
  }
  UG_LOG ("----> Total contribution by the rivers: " << total_recharge << '\n');
  
  if (! m_sTotalRechargeFile.empty ())
  {
#   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 (m_sTotalRechargeFile.c_str (), std::ofstream::app);
    UG_COND_THROW (output.fail (), "VertexRiverRecharge: Cannot open data file '" << m_sTotalRechargeFile << "' for output!");
    output << m_sOutputPrefix << ' ' << total_recharge << '\n';
  }
}
 
template <typename TGridFunction>
number VertexRiverRecharge<TGridFunction>::smoothed_recharge_at
(
  MathVector<dim> & x 
) const
{
  typedef Grid::traits<Edge>::const_iterator edge_iter_t;
  typedef SegmentNetwork::position_accessor_t river_vrt_pos_acc_t;
  
  Grid & riverGrid = river_network()->grid ();
  river_vrt_pos_acc_t river_vrt_pos_acc = river_network()->position_accessor ();
  
//  get the low-dimensional vector for x
  MathVector<dim-1> xc;
  for (size_t i = 0; i < dim-1; i++) xc[i] = x[i];
  
//  loop the segments of the river grid
  number sum = 0;
  for (edge_iter_t s_iter = riverGrid.begin<Edge>(); s_iter != riverGrid.end<Edge>(); ++s_iter)
  {
    Edge * seg = *s_iter;
    
  //  get the the center and the distance to the center
  
    Vertex * vrt_0 = seg->vertex (0);
    MathVector<dim-1> pos_0 = river_vrt_pos_acc[vrt_0];
    
    Vertex * vrt_1 = seg->vertex (1);
    MathVector<dim-1> pos_1 = river_vrt_pos_acc[vrt_1];
    
    MathVector<dim-1> seg_center_c;
    seg_center_c = pos_0; seg_center_c += pos_1; seg_center_c /= 2;
    
    number dist = VecDistance (xc, seg_center_c);
    if (dist >= this->smooth_len ())
      continue; // this segment is too far, we skip it! (Note that this imposes some assumptions to the mollifier.)
    
    number seg_length = VecDistance (pos_0, pos_1);
    
  //  get the depth of the free surface at the center
    number fs_depth = fs_depth_for (seg);
    if (fs_depth == - std::numeric_limits<number>::max ()) // this value should be set explicitly
      continue; // the segment is not under the free surface; no contribution
    
  //  get the depth of the river
    number seg_depth = river_depth_for (seg);
    if (seg_depth == std::numeric_limits<number>::max ()) // this value should be set explicitly
      continue; // the segment is not under the top of the domain; no contribution
    
  //  add the contribution
    number river_q = this->object_recharge (seg_depth, fs_depth, river_width_for (seg))
              * cone_mollifier<this_type::mol_dim> (dist, this->smooth_len()) * seg_length;
    if (this->no_inflow () && river_q >= 0)
      continue;
    sum += river_q;
  }
  
  return sum;
}
 
template <typename TGridFunction>
void VertexRiverRecharge<TGridFunction>::reset_river_depth ()
{
  typedef SegmentNetwork::position_accessor_t river_vrt_pos_acc_t;
  typedef SegmentNetwork::gauge_accessor_t river_depth_acc_t;
  Grid & riverGrid = m_spRiverNetwork->grid ();
  river_vrt_pos_acc_t river_vrt_pos_acc = river_network()->position_accessor ();
  river_depth_acc_t river_vrt_depth_acc = river_network()->gauge_accessor ();
  
  typedef Grid::traits<Edge>::const_iterator EdgeConstIterator;
  typedef Grid::traits<Vertex>::const_iterator VertConstIterator;
  
  typedef Grid::VertexAttachmentAccessor<depth_attachment_type> nodal_depth_accessor_type;
  depth_attachment_type aNodalFSDepth;
  nodal_depth_accessor_type aaNodalFSDepth;
  riverGrid.attach_to_vertices (aNodalFSDepth);
  aaNodalFSDepth.access (riverGrid, aNodalFSDepth);
  
//  Prepare the object for a new level-set function
  this->prepare_ls_height ();
  
//  Loop the river nodes
  VertConstIterator vertIterEnd = riverGrid.template end<Vertex> ();
  for (VertConstIterator vertIter = riverGrid.template begin<Vertex> (); vertIter != vertIterEnd; ++vertIter)
  {
    Vertex * vert = *vertIter;
    aaNodalFSDepth[vert] = this->get_ls_height_at (river_vrt_pos_acc [vert]);
  }
  
#ifdef UG_PARALLEL
//  reduce the depths
  pcl::ProcessCommunicator procComm;
  std::vector<number> fs_depth_vec (riverGrid.num_vertices ());
  size_t vert_i = 0;
  for (VertConstIterator vertIter = riverGrid.template begin<Vertex> (); vertIter != vertIterEnd; ++vertIter)
    fs_depth_vec [vert_i++] = aaNodalFSDepth [*vertIter];
  std::vector<number> red_fs_depth_vec;
  procComm.allreduce (fs_depth_vec, red_fs_depth_vec, PCL_RO_MAX);
  if (red_fs_depth_vec.size () != fs_depth_vec.size ())
    UG_THROW ("VertexRiverRecharge: Failed to reduce the free surface depth vector.");
  vert_i = 0;
  for (VertConstIterator vertIter = riverGrid.template begin<Vertex> (); vertIter != vertIterEnd; ++vertIter)
    aaNodalFSDepth [*vertIter] = red_fs_depth_vec [vert_i++];
#endif // UG_PARALLEL
 
//  Write the depth back to the rivers
  for (VertConstIterator vertIter = riverGrid.template begin<Vertex> (); vertIter != vertIterEnd; ++vertIter)
  {
    number depth = aaNodalFSDepth [*vertIter];
    if (depth != - std::numeric_limits<number>::max ()) // if the vertex is not covered by the FS, do not reset the orig. depth
      river_vrt_depth_acc[*vertIter] = depth;
  }
  
  riverGrid.detach_from_vertices (aNodalFSDepth);
}
 
template <typename TGridFunction>
void VertexRiverRecharge<TGridFunction>::save_river_depth
(
  const char * file_name 
)
{
  typedef SegmentNetwork::position_accessor_t river_vrt_pos_acc_t;
  typedef SegmentNetwork::gauge_accessor_t river_depth_acc_t;
  Grid & riverGrid = m_spRiverNetwork->grid ();
  river_vrt_pos_acc_t river_vrt_pos_acc = river_network()->position_accessor ();
  river_depth_acc_t river_vrt_depth_acc = river_network()->gauge_accessor ();
  
  typedef Grid::traits<Edge>::const_iterator EdgeConstIterator;
  typedef Grid::traits<Vertex>::const_iterator VertConstIterator;
  
  typedef Grid::VertexAttachmentAccessor<depth_attachment_type> nodal_depth_accessor_type;
  depth_attachment_type aNodalFSDepth;
  nodal_depth_accessor_type aaNodalFSDepth;
  riverGrid.attach_to_vertices (aNodalFSDepth);
  aaNodalFSDepth.access (riverGrid, aNodalFSDepth);
  
//  Prepare the object for a new level-set function
  this->prepare_ls_height ();
  
//  Loop the river nodes
  VertConstIterator vertIterEnd = riverGrid.template end<Vertex> ();
  for (VertConstIterator vertIter = riverGrid.template begin<Vertex> (); vertIter != vertIterEnd; ++vertIter)
  {
    Vertex * vert = *vertIter;
    aaNodalFSDepth[vert] = this->get_ls_height_at (river_vrt_pos_acc [vert]);
  }
  
#ifdef UG_PARALLEL
//  reduce the depths
  pcl::ProcessCommunicator procComm;
  std::vector<number> fs_depth_vec (riverGrid.num_vertices ());
  size_t vert_i = 0;
  for (VertConstIterator vertIter = riverGrid.template begin<Vertex> (); vertIter != vertIterEnd; ++vertIter)
    fs_depth_vec [vert_i++] = aaNodalFSDepth [*vertIter];
  std::vector<number> red_fs_depth_vec;
  procComm.allreduce (fs_depth_vec, red_fs_depth_vec, PCL_RO_MAX);
  if (red_fs_depth_vec.size () != fs_depth_vec.size ())
    UG_THROW ("VertexRiverRecharge: Failed to reduce the free surface depth vector.");
  vert_i = 0;
  for (VertConstIterator vertIter = riverGrid.template begin<Vertex> (); vertIter != vertIterEnd; ++vertIter)
    aaNodalFSDepth [*vertIter] = red_fs_depth_vec [vert_i++];
#endif // UG_PARALLEL
 
//  Save the file
  m_spRiverNetwork->save_to_file_with_ (aaNodalFSDepth, file_name);
  
  riverGrid.detach_from_vertices (aNodalFSDepth);
}
 
template <typename TGridFunction>
number VertexRiverRecharge<TGridFunction>::simple_smoothed_recharge_at
(
  MathVector<dim> & x 
) const
{
  typedef SegmentNetwork::position_accessor_t river_vrt_pos_acc_t;
  typedef SegmentNetwork::gauge_accessor_t river_depth_acc_t;
  typedef Grid::traits<Edge>::const_iterator edge_iter_t;
 
  Grid & riverGrid = river_network()->grid ();
  river_vrt_pos_acc_t river_vrt_pos_acc = river_network()->position_accessor();
  river_depth_acc_t river_vrt_depth_acc = river_network()->gauge_accessor();
 
//  get the low-dimensional vector for x
  MathVector<dim-1> xc;
  for (size_t i = 0; i < dim-1; i++) xc[i] = x[i];
 
 
  //  loop the segments of the river grid
    number sum = 0;
    for (edge_iter_t s_iter = riverGrid.begin<Edge>(); s_iter != riverGrid.end<Edge>(); ++s_iter)
    {
      Edge * seg = *s_iter;
 
    //  get the the center and the distance to the center
 
      Vertex * vrt_0 = seg->vertex (0);
      MathVector<dim-1> pos_0 = river_vrt_pos_acc[vrt_0];
      number depth_0 = river_vrt_depth_acc[vrt_0];
 
      Vertex * vrt_1 = seg->vertex (1);
      MathVector<dim-1> pos_1 = river_vrt_pos_acc[vrt_1];
      number depth_1 = river_vrt_depth_acc[vrt_1];
 
      MathVector<dim-1> seg_center_c;
      seg_center_c = pos_0; seg_center_c += pos_1; seg_center_c *= 0.5;
 
      number dist = VecDistance (xc, seg_center_c);
      if (dist >= this->smooth_len ())
        continue; // this segment is too far, we skip it! (Note that this imposes some assumptions to the mollifier.)
 
      number seg_length = VecDistance (pos_0, pos_1);
      number seg_depth = (depth_0 + depth_1)*0.5;
 
    //  add the contribution
      number river_q = seg_depth * cone_mollifier<this_type::mol_dim> (dist, this->smooth_len()) * seg_length * 0.5;
      //  this->object_recharge (seg_depth, fs_depth)
      if (this->no_inflow () && river_q >= 0) continue;
 
      sum += river_q;
    }
 
  return sum;
}
 
} // namespace d3f
} // end namespace ug
 
/* End of File */