val_surf_proj_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
 */
 
 
/*
 * Projection of values of a given grid function from a given z-plane to a given subset.
 */
#include <limits>
 
#include "lib_disc/local_finite_element/local_finite_element_provider.h"
#include "lib_disc/reference_element/reference_mapping_provider.h"
#ifdef UG_PARALLEL
#include "pcl/pcl_process_communicator.h"
#endif
 
namespace ug{
namespace d3f{
 
template <typename TGridFunction>
ValSurfProjData<TGridFunction>::ValSurfProjData
(
  SmartPtr<gf_type> spLSF, 
  const char * cmp, 
  number z_coord, 
  const char * subsetNames 
)
: m_spLSF (spLSF), m_z_coord (z_coord), m_glob_userdata (spLSF, cmp), m_default_value (0)
{
  m_SsGrp.set_subset_handler (spLSF->domain()->subset_handler ());
  m_SsGrp.add (TokenizeString (subsetNames));
  spLSF->domain()->grid()->attach_to_vertices_dv (m_a_vrt_num, (size_t) -1);
}
 
template <typename TGridFunction>
void ValSurfProjData<TGridFunction>::reinit ()
{
  typedef typename gf_type::template traits<Vertex>::const_iterator vrt_iter_t;
  
//  initialize the attachment
  num_aa_type aa_vrt_num (* m_spLSF->domain()->grid (), m_a_vrt_num);
  size_t vrt_num = 0;
  for (size_t i = 0; i < m_SsGrp.size (); i++)
  {
    int si = m_SsGrp[i];
    for (vrt_iter_t iter = m_spLSF->template begin<Vertex> (si);
            iter != m_spLSF->template end<Vertex> (si); ++iter)
      aa_vrt_num [*iter] = vrt_num++;
  }
  
# ifndef UG_PARALLEL
  
//  initialize the values
  m_vrt_val.resize (vrt_num);
  const pos_acc_type aa_pos = m_spLSF->domain()->position_accessor ();
  for (size_t i = 0; i < m_SsGrp.size (); i++)
  {
    int si = m_SsGrp[i];
    for (vrt_iter_t iter = m_spLSF->template begin<Vertex> (si);
            iter != m_spLSF->template end<Vertex> (si); ++iter)
    {
      Vertex * vrt = * iter;
      const MathVector<dim> & pos = aa_pos [vrt];
      MathVector<dim> the_pos;
      for (size_t k = 0; k < dim - 1; k++) the_pos[k] = pos[k];
      the_pos[dim-1] = m_z_coord;
      if (! m_glob_userdata.evaluate (m_vrt_val [aa_vrt_num [vrt]], the_pos))
        m_vrt_val [aa_vrt_num [vrt]] = m_default_value;
    }
  }
  
# else // i.e. ifdef UG_PARALLEL
  
  {
    pcl::ProcessCommunicator proc_comm;
    
  //  get the points where the values must be computed
    std::vector<MathVector<dim> > proc_pos (vrt_num);
    pos_acc_type aa_pos = m_spLSF->domain()->position_accessor ();
    for (size_t i = 0; i < m_SsGrp.size (); i++)
    {
      int si = m_SsGrp[i];
      for (vrt_iter_t iter = m_spLSF->template begin<Vertex> (si);
              iter != m_spLSF->template end<Vertex> (si); ++iter)
      {
        Vertex * vrt = * iter;
        const MathVector<dim> & pos = aa_pos [vrt];
        MathVector<dim> & the_pos = proc_pos [aa_vrt_num [vrt]];
        for (size_t k = 0; k < dim - 1; k++) the_pos[k] = pos[k];
        the_pos[dim-1] = m_z_coord;
      }
    }
    std::vector<MathVector<dim> > 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_vrt_val (all_pos.size ());
    for (size_t j = 0; j < tmp_vrt_val.size (); j++)
      if (! m_glob_userdata.evaluate (tmp_vrt_val [j], all_pos[j]))
        tmp_vrt_val [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_vrt_val (tmp_vrt_val.size ());
    proc_comm.allreduce (tmp_vrt_val, all_vrt_val, PCL_RO_MAX);
    
  //  extract the values from the common array
    m_vrt_val.resize (vrt_num);
    int my_offset = offset [proc_comm.get_local_proc_id ()];
    for (size_t j = 0; j < vrt_num; j++)
    {
      number val = all_vrt_val [my_offset + j];
      if (val == - std::numeric_limits<number>::max ()) // i.e. is not set in any process
        m_vrt_val [j] = m_default_value;
      else
        m_vrt_val [j] = all_vrt_val [my_offset + j];
    }
  }
  
# endif // UG_PARALLEL
}
 
template <typename TGridFunction>
template <int refDim>
void ValSurfProjData<TGridFunction>::eval_and_deriv
(
  number vValue[], 
  const MathVector<dim> vGlobIP[],
  number time,
  int si,
  GridObject* grid_obj,
  const MathVector<dim> vCornerCoords[],
  const MathVector<refDim> vLocIP[],
  const size_t nip,
  LocalVector* u,
  bool bDeriv,
  int s,
  std::vector<std::vector<number> > vvvDeriv[],
  const MathMatrix<refDim, dim>* 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 ("ValSurfProjData: 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));
 
  num_aa_type aa_vrt_num (* m_spLSF->domain()->grid (), m_a_vrt_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 (), "ValSurfProjData: 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 ("ValSurfProjData: The object is not properly initialized.");
      vValue [ip] += m_vrt_val [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];
    }
  }
}
 
} // end namespace d3f
} // end namespace ug
 
/* End of File */