domain_indirect_seq.cc

Go to the documentation of this file.

#include "rheolef/domain_indirect.h"
#include "rheolef/geo.h"
#include "rheolef/rheostream.h"
 
namespace rheolef {
 
// ----------------------------------------------------------------------------
// @brief build a domain from a list of element indexes
// ----------------------------------------------------------------------------
template<class M>
void 
domain_indirect_base_rep<M>::build_from_list (
    const std::string&            name,
    size_type                     map_dim,
    const communicator&           comm,
    const std::vector<size_type>& ie_list)
{
  _name    = name;
  _map_dim = map_dim;
  distributor dom_ownership (distributor::decide, comm, ie_list.size());
  disarray<geo_element_indirect,M>::resize (dom_ownership);
  typename disarray<geo_element_indirect,M>::iterator oige_iter = disarray<geo_element_indirect,M>::begin();
  for (std::vector<size_type>::const_iterator iter = ie_list.begin(),
                                              last = ie_list.end();
        iter != last; iter++, oige_iter++) {
    (*oige_iter).set (1, *iter);
  }
}
// ----------------------------------------------------------------------------
// @brief build the union of two domains
// ----------------------------------------------------------------------------
// c := a union b with c=*this
template<class M>
void
domain_indirect_base_rep<M>::build_union (
  const domain_indirect_base_rep<M>& a,
  const domain_indirect_base_rep<M>& b)
{
  check_macro (a._map_dim == b._map_dim,
        "union: domains "<<a._name<<" and " << b._name << " have incompatible dimensions");
  _name = a._name + "+" + b._name;
  _map_dim = a._map_dim;
  // use a map to perform the union efficiently
  std::map<size_type,geo_element_indirect> c_map;
  for (typename disarray<geo_element_indirect,M>::const_iterator
        iter = a.begin(), last = a.end(); iter != last; ++iter) {
    size_type ie = (*iter).index();
    c_map [ie] = *iter;
  }
  for (typename disarray<geo_element_indirect,M>::const_iterator
        iter = b.begin(), last = b.end(); iter != last; ++iter) {
    size_type ie = (*iter).index();
    c_map [ie] = *iter; // not inserted if already present
  }
  distributor c_ownership (distributor::decide, a.comm(), c_map.size());
  disarray<geo_element_indirect,M>::resize (c_ownership);
  typename disarray<geo_element_indirect,M>::iterator oige_iter = disarray<geo_element_indirect,M>::begin();
  for (typename std::map<size_type,geo_element_indirect>::const_iterator
        iter = c_map.begin(), last = c_map.end(); iter != last; ++iter, ++oige_iter) {
    *oige_iter = (*iter).second;
  }
}
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
odiststream&
domain_indirect_rep<sequential>::put (odiststream& ops)  const {
  using namespace std;
  ops << "domain" << endl
      << base::_name << endl
      << "2 " << base::_map_dim << " " << size() << endl;
  disarray<geo_element_indirect,sequential>::put_values (ops);
  return ops;
}
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
void
build_set_that_contains_S (
        const geo_element&            S,
        const std::vector<index_set>& ball,
        index_set&                    contains_S)
{
    typedef geo_element::size_type size_type;
    contains_S.clear();
    if (S.size() == 0) return;
    contains_S = ball[S[0]];
    for (size_type i = 1; i < S.size(); i++) {
        contains_S.inplace_intersection (ball[S[i]]);
    }
}
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
template<class U>
idiststream&
domain_indirect_rep<sequential>::get (
    idiststream&                      ips,
    const geo_rep<U,sequential>&      omega,
    std::vector<index_set>           *ball)
{ 
  std::istream& is = ips.is();
  if (!scatch(is,"\ndomain")) {
    is.setstate (std::ios::badbit);
    return ips;
  }
  size_type version, n_side;
  is >> base::_name >> version >> base::_map_dim >> n_side;
  check_macro (version <= 2, "unexpected domain format version="<<version);
  if (version == 2 || base::_map_dim == 0) {
    disarray<geo_element_indirect,sequential>::resize (n_side);
    disarray<geo_element_indirect,sequential>::get_values (ips);
    // check that elements are sorted by increasing variant
    {
      size_type prev_variant = 0;
      for (size_type ioige = 0, noige = size(); ioige < noige; ioige++) {
        size_type ige = oige (ioige).index();
        const geo_element& S = omega.get_geo_element (base::_map_dim, ige);
        check_macro (prev_variant <= S.variant(), "elements should be sorted by increasing variant order");
        prev_variant = S.variant();
      }
    }
    return ips;
  }
  check_macro (version == 1, "unexpected domain format version="<<version);
  check_macro (base::_map_dim <= omega.dimension(), "unexpected domain dimension="
        <<base::_map_dim<<" > geometry dimension = " << omega.dimension());
  //
  // old version=1 format: searh for index of sides
  // => this computation is slow, so format 2 is preferred 
  //    when using large 3d meshes with large boundary domains
  //    or with 3d regions defined as domains
  //
  geo_element_auto<> elem_init;
  typedef disarray<geo_element_auto<>,sequential> 
          disarray_t;
  disarray_t d_tmp (n_side, elem_init);
  d_tmp.get_values (ips);
  build_from_data (omega, d_tmp, ball);
  return ips;
}
template<class U>
void
domain_indirect_rep<sequential>::build_from_data (
    const geo_rep<U,sequential>&      omega,
    disarray<geo_element_auto<>,sequential>&
                                      d_tmp,
    std::vector<index_set>*           ball)
{
  typedef disarray<geo_element_auto<>,sequential> 
          disarray_t;
 
  disarray<geo_element_indirect,sequential>::resize (d_tmp.size());
  if (d_tmp.size() == 0) return;
  base::_map_dim = d_tmp[0].dimension(); // get map_dim from d_tmp disarray elements
 
  if (ball [base::_map_dim].size() == 0) {
    ball [base::_map_dim].resize (omega.n_vertex());
    for (size_type ige = 0, nge = omega.geo_element_ownership(base::_map_dim).size(); ige < nge; ige++) {
        const geo_element& E = omega.get_geo_element (base::_map_dim,ige);
        for (size_type iloc = 0; iloc < E.size(); iloc++) {
            ball [base::_map_dim] [E[iloc]] += ige;
        }
    }
  }
  disarray_t::const_iterator q = d_tmp.begin();
  size_type prev_variant = 0;
  for (iterator_ioige p = ioige_begin(), last = ioige_end(); p != last; ++p, ++q) {
    const geo_element& S = *q;
    check_macro (prev_variant <= S.variant(), "elements should be sorted by increasing variant order");
    prev_variant = S.variant();
    index_set contains_S;
    build_set_that_contains_S (S, ball[base::_map_dim], contains_S);
    check_macro (contains_S.size() >= 1, "domain element not in mesh: S.dis_ie=" << S.dis_ie());
    check_macro (contains_S.size() <= 1, "problem with domain element: S.dis_ie=" << S.dis_ie());
    // now, choose the only one mesh (sub-)element matching domain element S
    size_type i_side = *(contains_S.begin());
    const geo_element& S_orig = omega.get_geo_element (base::_map_dim, i_side);
    geo_element::orientation_type orient;
    geo_element::shift_type       shift;
    check_macro (S_orig.get_orientation_and_shift (S, orient, shift),
                        "problem with domain element: S.dis_ie=" << S.dis_ie());
    (*p).set (orient, i_side, shift);
  }
}
template <class U>
domain_indirect_basic<sequential>::domain_indirect_basic (
        disarray<geo_element_auto<>,sequential>&
                                        d_tmp,
        const geo_basic<U, sequential>& omega,
        std::vector<index_set>*         ball)
 : smart_pointer<domain_indirect_rep<sequential> > (new_macro(domain_indirect_rep<sequential>))
{
  check_macro (omega.variant() == geo_abstract_base_rep<U>::geo,
        "build a domain for a geo_domain: not yet");
  const geo_rep<U,sequential>& omega_data = dynamic_cast<const geo_rep<U,sequential>&>(omega.data());
  data().build_from_data (omega_data, d_tmp, ball);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template
idiststream&
domain_indirect_rep<sequential>::get (
    idiststream&                      ips,
    const geo_rep<Float,sequential>&  omega,
    std::vector<index_set>           *ball);
 
template
void
domain_indirect_rep<sequential>::build_from_data (
    const geo_rep<Float,sequential>&  omega,
    disarray<geo_element_auto<>,sequential>&
                                      d_tmp,
    std::vector<index_set>*           ball);
 
template
domain_indirect_basic<sequential>::domain_indirect_basic (
        disarray<geo_element_auto<>,sequential>&
                                        d_tmp,
        const geo_basic<Float, sequential>& omega,
        std::vector<index_set>*         ball);
 
#define _RHEOLEF_instanciation(M)                                       \
template class domain_indirect_base_rep<M>;                             \
 
_RHEOLEF_instanciation(sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(distributed)
#endif // _RHEOLEF_HAVE_MPI
 
} // namespace rheolef