fract_eval_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 evaluation of values of functions at fractures.
 */
 
namespace ug {
namespace d3f {
 
template <typename TGridFunc>
FractGFEval<TGridFunc>::FractGFEval
(
  SmartPtr<grid_func_type> spGridFunc, 
  const char * cmp_name, 
  SmartPtr<deg_layer_mngr_type> spDegLayerMngr, 
  const position_type & point_0, 
  const position_type & point_1 
)
: m_spGridFunc (spGridFunc), m_spDegLayerMngr (spDegLayerMngr),
  m_point_0 (point_0), m_point_1 (point_1)
{
  init_and_check (cmp_name);
}
 
#ifdef UG_FOR_LUA
template <typename TGridFunc>
FractGFEval<TGridFunc>::FractGFEval
(
  SmartPtr<grid_func_type> spGridFunc, 
  const char * cmp_name, 
  SmartPtr<deg_layer_mngr_type> spDegLayerMngr, 
  const std::vector<number> & point_0, 
  const std::vector<number> & point_1 
)
: m_spGridFunc (spGridFunc), m_spDegLayerMngr (spDegLayerMngr)
{
  if (point_0.size () != dim || point_1.size () != dim)
    UG_THROW ("FractGFEval: Dimension mismatch of the coordinates");
  for (size_t i = 0; i < dim; i++)
  {
    m_point_0 [i] = point_0 [i];
    m_point_1 [i] = point_1 [i];
  }
 
  init_and_check (cmp_name);
}
#endif
 
template <typename TGridFunc>
void FractGFEval<TGridFunc>::init_and_check
(
  const char * cmp_name 
)
{
//  check the data
  if (m_spGridFunc.invalid ())
    UG_THROW ("FractGFEval: No grid function specified.");
  if (m_spDegLayerMngr.invalid ())
    UG_THROW ("FractGFEval: No fracture manager specified.");
 
//  get the component
  m_fct = m_spGridFunc->fct_id_by_name (cmp_name);
  if (m_fct >= m_spGridFunc->num_fct ())
    UG_THROW ("FractGFEval: Function space does not contain a function with name '" << cmp_name << "'.");
    
//  local the finite element id
  m_lfeID = m_spGridFunc->local_finite_element_id (m_fct);
}
 
template <typename TGridFunc>
bool FractGFEval<TGridFunc>::is_intersected
(
  element_type * elem, 
  Grid & grid, 
  position_accessor_type & aaPos, 
  position_type & intersection, 
  size_t & num_fract_co, 
  side_type * & inner_side, 
  size_t inner_side_corners [], 
  size_t ass_co [] 
)
{
  side_type * outer_side;
  size_t inner_side_idx, outer_side_idx, outer_side_corners [maxLayerSideCorners];
  m_spDegLayerMngr->get_layer_sides (elem, num_fract_co,
    inner_side, inner_side_idx, inner_side_corners,
    outer_side, outer_side_idx, outer_side_corners, ass_co);
    
  MathVector<dim> dir;
  number lambda_min, lambda_max;
  VecSubtract (dir, m_point_1, m_point_0);
  if (RayElementIntersection (lambda_min, lambda_max, m_point_0, dir, inner_side, grid, aaPos))
  {
    if (std::fabs (lambda_max - lambda_min) > SMALL)
      UG_THROW ("FractGFEval: The line lyes directly in the fracture!");
    if (lambda_min < 0 || lambda_max > 1)
      return false; // not inside the segment
    intersection = m_point_0;
    VecScaleAppend (intersection, lambda_min, dir);
    return true;
  }
  
  return false;
}
 
template <typename TGridFunc>
void FractGFEval<TGridFunc>::get_local_coord
(
  side_type * side, 
  size_t num_co, 
  position_accessor_type & aaPos, 
  position_type & global, 
  MathVector<dim-1> & local 
)
{
  position_type corners [maxLayerSideCorners];
  for (size_t i = 0; i < num_co; i++)
    corners [i] = aaPos [side->vertex (i)];
  
  DimReferenceMapping<dim-1, dim> & rmap
    = ReferenceMappingProvider::get<dim-1, dim> (side->reference_object_id (), corners);
  rmap.global_to_local (local, global);
}
 
template <typename TGridFunc>
void FractGFEval<TGridFunc>::get_co_values
(
  element_type * elem, 
  side_type * side, 
  size_t num_fract_co, 
  size_t inner_side_corners [], 
  size_t ass_co [], 
  MathVector<dim-1> & local, 
  number & inner_val, 
  number & outer_val 
)
{
  const ReferenceObjectID roid = side->reference_object_id ();
  const LocalShapeFunctionSet<dim-1>& r_trial_space =
      LocalFiniteElementProvider::get<dim-1> (roid, m_lfeID);
  std::vector<number> shape (maxLayerSideCorners);
  r_trial_space.shapes (shape, local);
  if (shape.size () != num_fract_co)
    UG_THROW ("FractGFEval: Illegal number of shapes.");
  
  inner_val = outer_val = 0;
  std::vector<DoFIndex> inner_ind (1), outer_ind (1);
  for (size_t side_co = 0; side_co < num_fract_co; side_co++)
  {
    size_t inner_co = inner_side_corners [side_co];
    size_t outer_co = ass_co [inner_co];
    
    m_spGridFunc->dof_indices (elem->vertex (inner_co), m_fct, inner_ind);
    m_spGridFunc->dof_indices (elem->vertex (outer_co), m_fct, outer_ind);
    if (inner_ind.size () != 1 || outer_ind.size () != 1)
      UG_THROW ("FractGFEval: A non-scalar function specified.");
    inner_val += DoFRef (*m_spGridFunc, inner_ind[0]) * shape [side_co];
    outer_val += DoFRef (*m_spGridFunc, outer_ind[0]) * shape [side_co];
  }
}
 
template <typename TGridFunc>
void FractGFEval<TGridFunc>::append_values
(
  const side_type * side, 
  const position_type & intersection, 
  number inner_val, 
  number outer_val 
)
{
  for (size_t i = 0; i < m_values.size (); i++)
  {
    t_fract_pnt_data & data = m_values [i];
    if (data.face == side)
    {
      if (data.closed)
        UG_THROW ("FractGFEval: More than 2 fracture sides!");
      data.side_val[1] = outer_val;
      data.closed = true;
      return;
    }
  }
  m_values.push_back (t_fract_pnt_data (side, intersection, inner_val, outer_val));
}
 
template <typename TGridFunc>
void FractGFEval<TGridFunc>::check_values ()
{
  for (size_t i = 0; i < m_values.size (); i++)
    if (! m_values[i].closed)
      UG_THROW ("FractGFEval: Less than 2 fracture sides!");
}
 
template <typename TGridFunc>
void FractGFEval<TGridFunc>::get_all_values ()
{
  typedef typename grid_func_type::template dim_traits<dim>::const_iterator t_elem_iter;
  
  Grid & grid = (Grid &) * (m_spGridFunc->domain()->grid ());
  position_accessor_type & aaPos = m_spGridFunc->domain()->position_accessor ();
  
  m_values.clear ();
  for (size_t i_si = 0; i_si < m_spDegLayerMngr->num_subsets (); i_si++)
  {
    int si = m_spDegLayerMngr->subset (i_si);
    t_elem_iter iter_end = m_spGridFunc->template end<element_type> (si);
    
    for (t_elem_iter iter = m_spGridFunc->template begin<element_type> (si); iter != iter_end; ++iter)
    {
      element_type * elem = *iter;
      side_type * side;
      position_type intersection;
      size_t num_fract_co, inner_side_corners [maxLayerSideCorners], ass_co [2 * maxLayerSideCorners];
      
      if (! is_intersected (elem, grid, aaPos, intersection, num_fract_co, side, inner_side_corners, ass_co))
        continue;
      
      MathVector<dim-1> local;
      number inner_val, outer_val;
      get_local_coord (side, num_fract_co, aaPos, intersection, local);
      get_co_values (elem, side, num_fract_co, inner_side_corners, ass_co, local, inner_val, outer_val);
      
      append_values (side, intersection, inner_val, outer_val);
    }
  }
  
  check_values ();
}
 
template <typename TGridFunc>
void FractGFEval<TGridFunc>::write
(
  std::ostream & out, 
  const char * inner_sep, 
  const char * outer_sep 
) const
{
  for (size_t k = 0; k < num_intersections (); k++)
  {
    out << outer_sep;
    position_type pos = position (k);
    for (size_t i = 0; i < dim; i++)
      out << pos[i] << inner_sep;
    out << fract_value (k)
      << inner_sep << side_value_1 (k)
      << inner_sep << side_value_2 (k);
  }
}
 
} // namespace d3f
} // end namespace ug
 
/* End of File */