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
put
void put(idiststream &in, odiststream &out, bool do_proj, bool do_lumped_mass, bool def_fill_opt, size_type extract_id, const Float &scale_value, const std::pair< Float, Float > &u_range, render_type render)
Definition: branch.cc:500
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
rheolef::disarray::size_type
rep::base::size_type size_type
Definition: disarray.h:463
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::distributor::decide
static const size_type decide
Definition: distributor.h:76
rheolef::domain_indirect_base_rep
Definition: domain_indirect.h:71
rheolef::domain_indirect_base_rep::build_from_list
void build_from_list(const std::string &name, size_type map_dim, const communicator &comm, const std::vector< size_type > &ie_list)
Definition: domain_indirect_seq.cc:32
rheolef::domain_indirect_base_rep::build_union
void build_union(const domain_indirect_base_rep< M > &a, const domain_indirect_base_rep< M > &b)
Definition: domain_indirect_seq.cc:55
rheolef::domain_indirect_base_rep::size_type
geo_element_indirect::size_type size_type
Definition: domain_indirect.h:77
rheolef::domain_indirect_basic
the finite element boundary domain
Definition: domain_indirect.h:335
rheolef::domain_indirect_rep< sequential >::iterator_ioige
base::iterator_ioige iterator_ioige
Definition: domain_indirect.h:158
rheolef::domain_indirect_rep
Definition: domain_indirect.h:148
rheolef::geo_abstract_base_rep
abstract base interface class
Definition: geo.h:248
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
rheolef::geo_element_auto
Definition: geo_element.h:376
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_element::orientation_type
geo_element_indirect::orientation_type orientation_type
Definition: geo_element.h:134
rheolef::geo_element::get_orientation_and_shift
bool get_orientation_and_shift(const geo_element &S, orientation_type &orient, shift_type &shift) const
return orientation and shift between *this element and S
Definition: geo_element.cc:114
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element::size_type
reference_element::size_type size_type
Definition: geo_element.h:125
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::geo_element::shift_type
geo_element_indirect::shift_type shift_type
Definition: geo_element.h:135
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
rheolef::geo_rep
sequential mesh representation
Definition: geo.h:778
rheolef::index_set
Definition: index_set_header.icc:58
rheolef::index_set::inplace_intersection
void inplace_intersection(const index_set &b)
Definition: index_set_body.icc:129
rheolef::index_set::clear
void clear()
Definition: index_set_header.icc:152
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
rheolef::smart_pointer
see the smart_pointer page for the full documentation
Definition: smart_pointer.h:351
rheolef::std
Definition: vec_expr_v2.h:402
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::distributed
distributed
Definition: asr.cc:228
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
mkgeo_contraction.name
string name
Definition: mkgeo_contraction.sh:133
rheolef-config.version
version
Definition: rheolef-config.in:126
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::scatch
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
scatch: see the rheostream page for the full documentation
Definition: scatch.icc:52
rheolef::build_set_that_contains_S
void build_set_that_contains_S(const geo_element &S, const std::vector< index_set > &ball, index_set &contains_S)
builds a set of elements that all contain S.
Definition: domain_indirect_seq.cc:99
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
a
Definition: diffusion_isotropic.h:25
p
Definition: sphere.icc:25
rheolef::sequential
Definition: distributed.h:28