smile_rescale_impl.hpp

Go to the documentation of this file.

 
#include "lib_grid/algorithms/attachment_util.h"
#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
#include "lib_disc/spatial_disc/disc_util/geom_provider.h"
#ifdef UG_PARALLEL
#include "lib_grid/parallelization/util/attachment_operations.hpp"
#endif
 
namespace ug {
namespace smile {
 
template <typename TDomain, typename TAlgebra>
template <typename TElem>
void SmartRescaler<TDomain, TAlgebra>::get_scaling_factors_for_
(
  int si, 
  gf_type & gf, 
  aa_factor_type & aa_factor 
)
{
  typedef typename gf_type::template traits<TElem>::const_iterator iterator;
  typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
  typedef FV1Geometry<TElem, WDim> FVGeom_type;
  
  const size_t n_co = ref_elem_type::numCorners;
  const size_t n_scv = FVGeom_type::numSCV; // actually, we assume n_scv == n_co
  
  LocalIndices loc_ind;
  LocalVector loc_vec;
  
//  get position accessor
  domain_type & dom = * gf.domain ();
  position_accessor_type & aaPos = dom.position_accessor ();
  
  FVGeom_type & geo = GeomProvider<FVGeom_type>::get();
  
  position_type e_corners [n_co];
  Vertex * e_vrt [n_co];
  
// loop the elements in the subset
  iterator end_e_iter = gf.template end<TElem>(si);
  for (iterator e_iter = gf.template begin<TElem>(si); e_iter != end_e_iter; ++e_iter)
  {
    TElem * e = *e_iter;
    
  //  get vertices and extract corner coordinates
    for (size_t i = 0; i < n_co; i++)
      e_corners[i] = aaPos [e_vrt[i] = e->vertex (i)];
  
  //  compute the FV geometry
    try
    {
      geo.update (e, e_corners, dom.subset_handler().get ());
    }
    UG_CATCH_THROW("SmartRescaler: Cannot update the Finite Volume Geometry for an element.");
  
  //  get the local vector
    gf.indices (e, loc_ind);
    loc_vec.resize (loc_ind);
    GetLocalVector (loc_vec, gf);
    
  //  get the values of the factor
    number factor_values [n_scv];
    (* m_spFactor) (factor_values, e_corners, 0, si, e,
      geo.scv_global_ips (), geo.scv_local_ips (), n_scv, &loc_vec);
 
  //  evaluate the factor in at the corners
    for (size_t i = 0; i < n_scv; i++)
    {
      const typename FVGeom_type::SCV & scv = geo.scv (i);
      aa_factor [e_vrt [scv.node_id()]] += factor_values[i] * scv.volume ();
    }
  }
}
 
template <typename TDomain, typename TAlgebra>
void SmartRescaler<TDomain, TAlgebra>::get_scaling_factors
(
  gf_type & gf, 
  a_factor_type & a_factor 
)
{
  domain_type & dom = * gf.domain ();
  grid_type & g = * dom.grid ();
  
  if (m_spFactor.invalid ())
    UG_THROW ("SmartRescaler: No factor specified!");
  
  m_spFactor->set_function_pattern (gf.function_pattern ());
  
  aa_factor_type aa_factor (g, a_factor);
  SetAttachmentValues (aa_factor, g.template begin<Vertex> (), g.template end<Vertex> (), 0);
  
//  The full-dim. grid element types for this dimension:
  typedef typename domain_traits<WDim>::DimElemList ElemList;
  subset_handler_type & ssh = * dom.subset_handler();
  for (int si = 0; si < ssh.num_subsets (); si++)
    boost::mpl::for_each<ElemList> (GetScalingFactors (this, si, gf, aa_factor));
  
# ifdef UG_PARALLEL
  AttachmentAllReduce<Vertex> (g, a_factor, PCL_RO_SUM);
# endif
}
 
template <typename TDomain, typename TAlgebra>
void SmartRescaler<TDomain, TAlgebra>::multiply
(
  SmartPtr<gf_type> gf, 
  const char * cmp 
)
{
  typedef typename gf_type::template traits<Vertex>::const_iterator iterator;
  
//  get the components
  size_t fct = gf->fct_id_by_name (cmp);
  
//  create the attachment
  domain_type & dom = * gf->domain ();
  grid_type & g = * dom.grid ();
  a_factor_type a_factor;
  g.attach_to_vertices (a_factor);
  
//  compute the scaling factors
  get_scaling_factors (*gf, a_factor);
  aa_factor_type aa_factor (g, a_factor);
  
  iterator end_v_iter = gf->template end<Vertex> ();
  for (iterator v_iter = gf->template begin<Vertex> (); v_iter != end_v_iter; ++v_iter)
  {
    Vertex * vrt = *v_iter;
    number factor = aa_factor [vrt];
    
    std::vector<DoFIndex> ind (1);
    gf->inner_dof_indices (vrt, fct, ind);
    DoFRef(*gf, ind[0]) *= factor; // this '*=' should be the only difference to 'divide'
  }
  
  g.detach_from_vertices (a_factor);
}
  
template <typename TDomain, typename TAlgebra>
void SmartRescaler<TDomain, TAlgebra>::divide
(
  SmartPtr<gf_type> gf, 
  const char * cmp 
)
{
  typedef typename gf_type::template traits<Vertex>::const_iterator iterator;
  
//  get the components
  size_t fct = gf->fct_id_by_name (cmp);
  
//  create the attachment
  domain_type & dom = * gf->domain ();
  grid_type & g = * dom.grid ();
  a_factor_type a_factor;
  g.attach_to_vertices (a_factor);
  
//  compute the scaling factors
  get_scaling_factors (*gf, a_factor);
  aa_factor_type aa_factor (g, a_factor);
  
  iterator end_v_iter = gf->template end<Vertex> ();
  for (iterator v_iter = gf->template begin<Vertex> (); v_iter != end_v_iter; ++v_iter)
  {
    Vertex * vrt = *v_iter;
    number factor = aa_factor [vrt];
    
    std::vector<DoFIndex> ind (1);
    gf->inner_dof_indices (vrt, fct, ind);
    DoFRef(*gf, ind[0]) /= factor; // this '/=' should be the only difference to 'multiply'
  }
  
  g.detach_from_vertices (a_factor);
}
 
} // end namespace smile
} // end namespace ug
 
/* End of File */