geo_build_from_list.cc

Go to the documentation of this file.

// build geo from element & node lists
#include "rheolef/geo.h"
#include "rheolef/geo_domain.h"
#include "rheolef/dis_macros.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
#include "rheolef/index_set.h"
 
namespace rheolef {
 
// --------------------------------------------------------------------------
// base class, common to seq & dist derived classes
// --------------------------------------------------------------------------
template <class T, class M>
void
geo_base_rep<T,M>::build_from_list (
  const geo_basic<T,M>&                               lambda,
  const disarray<point_basic<T>,M>&                   node_list,
  const std::array<disarray<geo_element_auto<heap_allocator<size_type> >,M>,
                     reference_element::max_variant>& ge_list)
{
  typedef geo_element_auto<heap_allocator<size_type> > geo_element_tmp_type;
  // ------------------
  // 0) header data
  // ------------------
  _name          = "level_set_from_" + lambda.name(); // TODO: send name by fct argument
  _version       = 4;
  _dimension     = lambda.dimension();
  _piola_basis   = lambda.get_piola_basis();
  _sys_coord     = lambda.coordinate_system();
  _have_connectivity = true;
  
  // compute map_dimension from data:
  _gs._map_dimension = 0;
  for (size_type variant = 0; variant < reference_element::max_variant; variant++) {
    if (ge_list[variant].dis_size() != 0) {
      _gs._map_dimension = reference_element::dimension(variant);
    }
  }
  // ------------------
  // 1) nodes
  // ------------------
  communicator comm = lambda.sizes().node_ownership.comm();
  _node = node_list;
  _gs.node_ownership = _node.ownership();
  compute_bbox();
  // vertices=geo_element[p] differs from nodes when order > 1:
  check_macro (order() == 1, "order > 1: not yet supported"); // TODO: with high order level set
  _gs.ownership_by_variant [reference_element::p] = _gs.node_ownership;
  // ------------------
  // 2) count elements
  // ------------------
  size_type ne = 0;
  for (size_type variant = reference_element::first_variant_by_dimension(_gs._map_dimension);
                 variant < reference_element:: last_variant_by_dimension(_gs._map_dimension); variant++) {
    size_type nge = ge_list [variant].size();
    _gs.ownership_by_variant [variant] = distributor (distributor::decide, comm, nge);
    geo_element::parameter_type param (variant, order());
    _geo_element [variant].resize (_gs.ownership_by_variant [variant], param);
    ne += nge;
  }
  _gs.ownership_by_dimension [_gs._map_dimension] = distributor (distributor::decide, comm, ne);
  //
  // 1.4) create 0d vertex-elements
  //
  // set ios_dis_iv index as fisrt field of the idx_vertex pair:
  // # of node that are vertices:
  {
    check_macro (order() == 1, "order > 1: not yet supported");
    _gs.ownership_by_dimension [0] = _gs.node_ownership;
    _gs.ownership_by_variant   [reference_element::p] = _gs.ownership_by_dimension [0];
    geo_element::parameter_type param (reference_element::p, order());
    _geo_element [reference_element::p].resize (_gs.ownership_by_dimension [0], param);
    size_type first_dis_iv = _gs.ownership_by_dimension [0].first_index();
    for (size_type iv = 0, nv = _gs.ownership_by_dimension [0].size(); iv < nv; iv++) {
      geo_element& P = _geo_element [reference_element::p] [iv];
      size_type dis_iv = first_dis_iv + iv;
      P [0] = dis_iv;
      P.set_dis_ie     (dis_iv);
      P.set_ios_dis_ie (dis_iv);
    }
  }
  // ------------------
  // 3) copy elements
  // ------------------
  size_type first_dis_ie = _gs.ownership_by_dimension [_gs._map_dimension].first_index();
  size_type dis_ie = first_dis_ie;
  size_type first_dis_v = 0;
  for (size_type variant = reference_element::first_variant_by_dimension(_gs._map_dimension);
                 variant < reference_element:: last_variant_by_dimension(_gs._map_dimension); variant++) {
    size_type first_dis_igev = _gs.ownership_by_variant [variant].first_index();
    size_type dis_igev = first_dis_igev;
    hack_array<geo_element>::iterator ge_iter = _geo_element [variant].begin();
    for (typename disarray<geo_element_tmp_type,M>::const_iterator
          iter = ge_list[variant].begin(),
          last = ge_list[variant].end();
          iter != last; iter++, ge_iter++, dis_ie++, dis_igev++) {
      *ge_iter = *iter;
      size_type ios_dis_ie = first_dis_v + dis_igev;
      (*ge_iter).set_dis_ie     (dis_ie);
      (*ge_iter).set_ios_dis_ie (ios_dis_ie);
    }
    first_dis_v += _gs.ownership_by_variant[variant].dis_size();
  }
  if (_gs._map_dimension >= 2) {
    _have_connectivity = false; // TODO
    trace_macro ("connectivity: not yet computed");
  }
}
// --------------------------------------------------------------------------
// seq
// --------------------------------------------------------------------------
template <class T>
geo_rep<T,sequential>::geo_rep (
  const geo_basic<T,sequential>&                      lambda,
  const disarray<point_basic<T>,sequential>&          node_list,
  const std::array<disarray<geo_element_auto<heap_allocator<size_type> >,sequential>,
                     reference_element::max_variant>& ge_list)
 : geo_base_rep<T,sequential>()
{
  base::build_from_list (lambda, node_list, ge_list);
}
// --------------------------------------------------------------------------
// dist
// --------------------------------------------------------------------------
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
geo_rep<T,distributed>::geo_rep (
  const geo_basic<T,distributed>&                     lambda,
  const disarray<point_basic<T>,distributed>&         node_list,
  const std::array<disarray<geo_element_auto<heap_allocator<size_type> >,distributed>,
                     reference_element::max_variant>& ge_list)
 : geo_base_rep<T,distributed>(),
   _inod2ios_dis_inod(),
   _ios_inod2dis_inod(),
   _ios_ige2dis_ige(),
   _ios_gs(),
   _igev2ios_dis_igev(),
   _ios_igev2dis_igev()
{
  base::build_from_list (lambda, node_list, ge_list);
  // ------------------------------------------------------------------------
  // 1) set external entities, at partition boundaries
  // ------------------------------------------------------------------------
  build_external_entities ();
  // ------------------------------------------------------------------------
  // 2) set ios numbering
  // TODO: ios: not yet nproc independant
  // ------------------------------------------------------------------------
  communicator comm = lambda.sizes().node_ownership.comm();
  if (comm.size() > 1) { dis_trace_macro ("ios: not yet fully computed"); }
  _ios_gs = base::_gs;
  //
  // 2.1) node ios (TODO)
  //
  _inod2ios_dis_inod.resize (base::_node.ownership());
  size_type first_dis_inod = base::_node.ownership().first_index();
  for (size_type inod = 0, nnod = _inod2ios_dis_inod.size(); inod < nnod; inod++) {
    _inod2ios_dis_inod [inod] = first_dis_inod + inod;
  } 
  _ios_inod2dis_inod = _inod2ios_dis_inod;
  //
  // 2.2) 0d-elts = vertices (TODO)
  //
  _igev2ios_dis_igev[reference_element::p].resize (base::_gs.ownership_by_variant [reference_element::p]);
  size_type first_dis_iv = base::_gs.ownership_by_variant [reference_element::p].first_index();
  for (size_type iv = 0, nv = _igev2ios_dis_igev[reference_element::p].size(); iv < nv; iv++) {
    _igev2ios_dis_igev [reference_element::p][iv] = first_dis_iv + iv;
  } 
  _ios_igev2dis_igev[reference_element::p] = _igev2ios_dis_igev[reference_element::p];
  //
  // 2.3) element ios (TODO)
  //
  _ios_ige2dis_ige[base::_gs._map_dimension].resize  (base::_gs.ownership_by_dimension [base::_gs._map_dimension]);
  size_type first_dis_ige = base::_gs.ownership_by_dimension [base::_gs._map_dimension].first_index();
  size_type ige = 0;
  for (size_type variant = reference_element::first_variant_by_dimension(base::_gs._map_dimension);
                 variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension); variant++) {
    _igev2ios_dis_igev[variant].resize (base::_gs.ownership_by_variant [variant]);
    size_type first_dis_igev = base::_gs.ownership_by_variant [variant].first_index();
    for (size_type igev = 0, ngev = base::_geo_element[variant].size(); igev < ngev; igev++, ige++) {
      const geo_element& K = base::_geo_element [variant][igev];
      _ios_ige2dis_ige [base::_gs._map_dimension][ige] = first_dis_ige + ige;
      _igev2ios_dis_igev [variant][igev] = first_dis_igev + igev;
    }
    _ios_igev2dis_igev[variant] = _igev2ios_dis_igev [variant];
  }
}
#endif // _RHEOLEF_HAVE_MPI
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M)                                             \
template                                                                        \
void                                                                            \
geo_base_rep<T,M>::build_from_list (                                            \
  const geo_basic<T,M>&                               lambda,                   \
  const disarray<point_basic<T>,M>&                   node_list,                \
  const std::array<disarray<geo_element_auto<heap_allocator<size_type> >,M>,    \
                     reference_element::max_variant>& ge_list);                 \
template                                                                        \
geo_rep<T,M>::geo_rep (                                                         \
  const geo_basic<T,M>&                               lambda,                   \
  const disarray<point_basic<T>,M>&                      node_list,             \
  const std::array<disarray<geo_element_auto<heap_allocator<size_type> >,M>,    \
                     reference_element::max_variant>& ge_list);                 \
 
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
 
} // namespace rheolef