fract_lexorder_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 ordering of vertices (dofs) for domains with degenerated
 * fractures.
 */
 
#include "lib_disc/ordering_strategies/algorithms/lexorder.h"
#include "lib_disc/function_spaces/dof_position_util.h"
 
namespace ug {
namespace d3f {
 
template <typename TDomain>
void CorrectFractPositions
(
  ConstSmartPtr<TDomain> domain, 
  SmartPtr<DoFDistribution> dd, 
  std::vector<typename std::pair<MathVector<TDomain::dim>, size_t> >& vPositions, 
  DegeneratedLayerManager<TDomain::dim>* dLayerManager, 
  number step 
)
{
  static const int dim = TDomain::dim;
  typedef typename DegeneratedLayerManager<dim>::element_type element_type;
  typedef typename DegeneratedLayerManager<dim>::side_type side_type;
  typedef typename DoFDistribution::traits<element_type>::const_iterator t_elem_iter;
  static const size_t maxLayerSideCorners = DegeneratedLayerManager<dim>::maxLayerSideCorners;
  
//  do we have fractures?
  if (dLayerManager == NULL) return;
  
//  algebra indices vector
  std::vector<size_t> ind;
  
//  loop over fracture elements
  for (size_t ss_grp_i = 0; ss_grp_i < dLayerManager->subset_grp().size (); ss_grp_i++)
  {
    size_t si = dLayerManager->subset_grp() [ss_grp_i];
    t_elem_iter e_end = dd->template end<element_type> (si);
    for (t_elem_iter e_iter = dd->template begin<element_type> (si);
      e_iter != e_end; ++e_iter)
    {
      element_type * elem = *e_iter;
      size_t n_co, inner_side_idx, outer_side_idx;
      size_t inner_side_corners [maxLayerSideCorners];
      size_t outer_side_corners [maxLayerSideCorners];
      side_type * inner_side, * outer_side;
      size_t ass_co [maxLayerSideCorners];
      
    //  get the fracture sides
      dLayerManager->get_layer_sides (elem, n_co,
        inner_side, inner_side_idx, inner_side_corners,
        outer_side, outer_side_idx, outer_side_corners, ass_co);
      
    //  get the vertices of the outer side (ordered as in the ref. elem. to get the OUTER normal)
      const ReferenceElement & ref_elem = ReferenceElementProvider::get (elem->reference_object_id ());
      ReferenceObjectID sideRoid = ref_elem.roid (dim - 1, outer_side_idx);
      size_t n_side_co = ref_elem.num (dim - 1, outer_side_idx, 0);
      Vertex * side_vert [maxLayerSideCorners];
      for (size_t co = 0; co < n_side_co; co++)
        side_vert [co] = elem->vertex (ref_elem.id (dim - 1, outer_side_idx, 0, co));
      
    //  get the outer normal of the fracture
      const typename TDomain::position_accessor_type & aaPos = domain->position_accessor();
      MathVector<dim> vSideCoords [maxLayerSideCorners];
      MathVector<dim> outer_normal;
      number outer_normal_norm;
      for (size_t co = 0; co < n_side_co; ++co)
        vSideCoords [co] = aaPos [side_vert [co]];
      ElementNormal<dim> (sideRoid, outer_normal, vSideCoords);
      if ((outer_normal_norm = VecLength (outer_normal)) < 1e-32)
        UG_THROW ("Cannot get the unit normal to a fracture.");
      
    //  Compute the reference distance: The minimum distance between the corners of the side
      UG_ASSERT (n_side_co >= 2, "To few corners of a side: 1d?");
      number min_edge, t;
      min_edge = VecDistanceSq (vSideCoords[0], vSideCoords[1]);
      for (size_t i = 0; i < n_side_co; i++)
        for (size_t j = i + 1; j < n_side_co; j++)
          if (min_edge > (t = VecDistanceSq (vSideCoords[i], vSideCoords[j])))
            min_edge = t;
      if ((min_edge = std::sqrt (min_edge)) < 1e-32)
        UG_THROW ("Too short edges in the element!");
      
    //  correct the corresponding positions
      outer_normal *= step * min_edge / outer_normal_norm;
      for (size_t co = 0; co < n_side_co; ++co)
      {
        dd->inner_algebra_indices (side_vert [co], ind);
        for(size_t i = 0; i < ind.size (); i++)
          vPositions [ind[i]].first += outer_normal;
      }
    }
  }
}
 
template <typename TDomain>
void FractOrderLexForDofDist
(
  SmartPtr<DoFDistribution> dd, 
  ConstSmartPtr<TDomain> domain, 
  DegeneratedLayerManager<TDomain::dim>* dLayerManager, 
  number step 
)
{
  static const int dim = TDomain::dim;
  
//  We work with Lagrange P1 only:
  for (size_t fct = 0; fct < dd->num_fct(); fct++)
    if (dd->local_finite_element_id (fct) != LFEID (LFEID::LAGRANGE, dim, 1))
      UG_THROW ("FractOrderLexForDofDist: Lex. order for fract. domains implemented for LagrangeP1 only!");
  
//  get positions
  std::vector<typename std::pair<MathVector<TDomain::dim>, size_t> > vPositions;
  ExtractPositions<TDomain> (domain, dd, vPositions);
  
//  correct the positions at the fractures
  CorrectFractPositions (domain, dd, vPositions, dLayerManager, step);
 
//  get mapping: old -> new index
  std::vector<size_t> vNewIndex (dd->num_indices());
  ComputeLexicographicOrder<TDomain::dim> (vNewIndex, vPositions);
 
//  reorder indices
  dd->permute_indices (vNewIndex);
}
 
template <typename TDomain>
void FractOrderLex
(
  ApproximationSpace<TDomain>& approxSpace, 
  DegeneratedLayerManager<TDomain::dim>* dLayerManager, 
  number step 
)
{
  std::vector<SmartPtr<DoFDistribution> > vDD = approxSpace.dof_distributions ();
 
  for(size_t i = 0; i < vDD.size (); i++)
    FractOrderLexForDofDist<TDomain> (vDD[i], approxSpace.domain (), dLayerManager, step);
}
 
} // namespace d3f
} // end namespace ug
 
/* End of File */