geo_neighbour.cc
Go to the documentation of this file.
// Banded level set routines
//
// Author: Pierre Saramito
//
#include "rheolef/geo.h"
namespace rheolef {
// ----------------------------------------------------------------------------
// 1) utilities for neighbour initialization
// ----------------------------------------------------------------------------
// Implementation note: as member function in geo_base_rep<T,M> could be simpler
// but function member of template class cannot be specialized
// and the sequential case do not define scatter_map_type and such
// so use a specialization of a non-member function
template<class T, class M>
static
void
add_ball_externals (
const geo_base_rep<T,M>& omega,
const disarray<index_set,M>& ball)
{
}
// explicit M=sequential specialization (otherwise: linker problems)
template<class T>
static
void
add_ball_externals (
const geo_base_rep<T,sequential>& omega,
const disarray<index_set,sequential>& ball)
{
}
#ifdef _RHEOLEF_HAVE_MPI
template<class T>
static
void
add_ball_externals (
const geo_base_rep<T,distributed>& omega,
const disarray<index_set,distributed>& ball)
{
typedef typename geo_base_rep<T,distributed>::size_type size_type;
index_set ext_iv_set;
// faudrait attraper _geo_element[0].get_dis_map_entries() : suit les vertices au lieu nodes
const hack_array<geo_element_hack,distributed>& vertex = omega._geo_element[0];
typedef typename hack_array<geo_element_hack,distributed>::scatter_map_type map_t;
const map_t& map_vertex = vertex.get_dis_map_entries();
size_type dis_nv = omega.dis_size (0);
for (typename map_t::const_iterator
iter = map_vertex.begin(),
last = map_vertex.end(); iter != last; ++iter) {
size_type dis_iv = (*iter).first;
check_macro (dis_iv < dis_nv, "dis_iv="<<dis_iv<<" out of range [0:"<<dis_nv<<"[");
ext_iv_set += dis_iv;
}
ball.set_dis_indexes (ext_iv_set);
}
#endif // _RHEOLEF_HAVE_MPI
// ----------------------------------------------------------------------------
// 2) (lazy) initialization of neighbour informations
// ----------------------------------------------------------------------------
template<class T, class M>
void
geo_base_rep<T,M>::init_neighbour() const
{
if (is_broken()) return; // skip broken domain e.g. "sides" or "internal_sides": no neighbours
size_type map_dim = map_dimension();
if (map_dim <= 0) return; // a set of points have no neighbours
size_type sid_dim = map_dim-1;
// 0.0) clean-up neighbour information (could come from a copy of elements, from domain, e.g. band)
for (const_iterator iter_S = begin(sid_dim), last_S = end(sid_dim); iter_S != last_S; ++iter_S) {
const geo_element& S = *iter_S;
S.set_master (0, std::numeric_limits<size_type>::max());
S.set_master (1, std::numeric_limits<size_type>::max());
}
// 0.1) manage all external sides S that appears on some local elements K
distributor is_ownership = _gs.ownership_by_dimension[map_dim-1];
std::array<index_set, reference_element::max_variant> ext_isv_set; // external sides, by variants
index_set ext_is_set; // external sides, all variants merged
for (const_iterator iter_K = begin(map_dim), last_K = end(map_dim); iter_K != last_K; ++iter_K) {
const geo_element& K = *iter_K;
for (size_type loc_is = 0, loc_ns = K.n_subgeo(sid_dim); loc_is < loc_ns; ++loc_is) {
size_type dis_is = (map_dim == 1) ? K[loc_is] : ((map_dim == 2) ? K.edge(loc_is) : K.face(loc_is));
if (is_ownership.is_owned (dis_is)) continue;
// have an external side S: convert its dis_is to dis_isv and get its variant
size_type variant;
size_type dis_isv = _gs.dis_ige2dis_igev_by_dimension (sid_dim, dis_is, variant);
ext_isv_set [variant] += dis_isv;
ext_is_set += dis_is;
}
}
// for external side S that appears in a local element K may be available
// => enlarge the external set of _geo_element for sid_dim
// but ext_is_set may be classified by variants and all dis_is may be
// converted to dis_isv
for (size_type variant = reference_element::first_variant_by_dimension(sid_dim);
variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
_geo_element [variant].append_dis_indexes (ext_isv_set [variant]);
}
// 1) build ball(xi) := { K; xi is a vertex of K }
index_set emptyset;
distributor vertex_ownership = sizes().ownership_by_dimension[0];
disarray<index_set,M> ball (vertex_ownership, emptyset);
for (const_iterator iter = begin(map_dim), last = end(map_dim); iter != last; ++iter) {
const geo_element& K = *iter;
index_set dis_ie_set;
dis_ie_set += K.dis_ie();
for (size_type loc_inod = 0, loc_nnod = K.size(); loc_inod < loc_nnod; loc_inod++) {
size_type dis_inod = K [loc_inod];
size_type dis_iv = dis_inod2dis_iv (dis_inod);
ball.dis_entry (dis_iv) += dis_ie_set; // union with {dis_ie}
}
}
ball.dis_entry_assembly();
// 2) set ball available at partition boundaries
add_ball_externals (*this, ball);
// 3) set master elements on all sides
distributor ie_ownership = _gs.ownership_by_dimension[map_dim];
distributor inod_ownership = _gs.node_ownership;
size_type first_dis_ie = ie_ownership.first_index();
std::array<index_set, reference_element::max_variant> ext_igev_set; // external elts, by variants
index_set ext_ie_set; // external elements, all variants merged
index_set ext_inod_set;
disarray<size_type, M> side_marked (is_ownership, 0); // boolean: 0,1
side_marked.set_dis_indexes (ext_is_set);
for (const_iterator iter_K = begin(map_dim), last_K = end(map_dim); iter_K != last_K; ++iter_K) {
const geo_element& K = *iter_K;
for (size_type loc_is = 0, loc_ns = K.n_subgeo(sid_dim); loc_is < loc_ns; ++loc_is) {
size_type dis_is = (sid_dim == 0) ? K[loc_is] : ((sid_dim == 1) ? K.edge(loc_is) : K.face(loc_is));
if (side_marked.dis_at (dis_is)) continue;
side_marked.dis_entry (dis_is) = 1;
const geo_element& S = dis_get_geo_element (sid_dim, dis_is);
size_type dis_iv0 = dis_inod2dis_iv (S[0]);
index_set ball_intersect = ball.dis_at (dis_iv0);
for (size_type loc_iv = 1, loc_nv = S.size(); loc_iv < loc_nv; loc_iv++) {
size_type dis_iv = dis_inod2dis_iv (S[loc_iv]);
ball_intersect.inplace_intersection (ball.dis_at (dis_iv));
}
switch (ball_intersect.size()) {
case 1: { // one master K: S is on the boundary
index_set::const_iterator iter = ball_intersect.begin();
size_type dis_ie = *iter;
S.set_master (0, dis_ie);
S.set_master (1, std::numeric_limits<size_type>::max());
break;
}
case 2: { // two masters K1, K2: S is an internal side
index_set::const_iterator iter = ball_intersect.begin();
size_type dis_ie1 = *iter++;
size_type dis_ie2 = *iter;
// one of K1 and K2 is available in distributed environment
// the other is not a priori available
if (! ie_ownership.is_owned (dis_ie1)) std::swap (dis_ie1, dis_ie2);
if (! ie_ownership.is_owned (dis_ie2)) {
// convert dis_ie2 to dis_igev2 and get its variant2:
size_type variant2;
size_type dis_igev2 = _gs.dis_ige2dis_igev_by_dimension (map_dim, dis_ie2, variant2);
ext_igev_set [variant2] += dis_igev2;
ext_ie_set += dis_ie2;
}
check_macro (ie_ownership.is_owned (dis_ie1), "unexpected external K1");
// let K1 indicated by the normal on S given by S orient
// let K2, the other one
size_type ie1 = dis_ie1 - first_dis_ie;
const geo_element& K1 = get_geo_element (map_dim, ie1);
geo_element::orientation_type orient = K1.get_side_orientation (S);
if (orient < 0) std::swap (dis_ie1, dis_ie2);
S.set_master (0, dis_ie1);
S.set_master (1, dis_ie2);
break;
}
default: {
error_macro ("neighbour: side #"<< ball_intersect.size()<<" connectivity problem");
}
}
}
}
// side values have been changed:
// propagate these modifs for sides at the partition boundaries
for (size_type variant = reference_element::first_variant_by_dimension(sid_dim);
variant < reference_element:: last_variant_by_dimension(sid_dim); variant++) {
_geo_element [variant].update_dis_entries();
}
// for all S in the partition boundary, both K1 & K2 may be available
// => enlarge the external set of _geo_element for map_dim
// but ext_ie_set may be classified by variants and all dis_ie may be
// converted to dis_igev
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
_geo_element [variant].append_dis_indexes (ext_igev_set [variant]);
}
// manage also external nodes of external elements K:
for (index_set::const_iterator iter = ext_ie_set.begin(), last = ext_ie_set.end(); iter != last; ++iter) {
size_type dis_ie = *iter;
const geo_element& K = dis_get_geo_element (map_dim, dis_ie);
for (size_type loc_inod = 0, loc_nnod = K.n_node(); loc_inod < loc_nnod; loc_inod++) {
size_type dis_inod = K[loc_inod];
if (! inod_ownership.is_owned (dis_inod)) {
ext_inod_set += dis_inod;
}
}
}
_node.append_dis_indexes (ext_inod_set);
}
// ----------------------------------------------------------------------------
// 3) compute the neighbour K2 of an element K1 across a side S
// ----------------------------------------------------------------------------
// returns the dis_ie2 global index associated to K2,
// the neighbour element of K1, with index ie1,
// across the side S of local index loc_isid in K1
// when S is on the boundary, returns size_t(-1)
template<class T, class M>
typename geo_base_rep<T,M>::size_type
geo_base_rep<T,M>::neighbour (size_type ie1, size_type loc_isid) const
{
check_macro (_have_neighbour, "neighbour_guard() may be called globally before local neighbour()");
// 1) get the i-th side S in K1
size_type map_dim = map_dimension();
check_macro (ie1 < size(map_dim), "neighbour: element index ie=" << ie1
<< " is out of range [0:" << size(map_dim) << "[");
const geo_element& K1 = get_geo_element (map_dim, ie1);
size_type dis_isid;
switch (map_dim) {
case 3: dis_isid = K1.face (loc_isid); break;
case 2: dis_isid = K1.edge (loc_isid); break;
case 1: {
size_type dis_inod = K1[loc_isid];
size_type dis_iv = dis_inod2dis_iv (dis_inod);
dis_isid = dis_iv;
break;
}
case 0:
default: return std::numeric_limits<size_type>::max();
}
size_type sid_dim = map_dim-1;
const geo_element& S = dis_get_geo_element (sid_dim, dis_isid);
// 2) get the neighbour element K2 to K1 across S
size_type dis_ie2 = S.master (0);
size_type first_dis_ie = _gs.ownership_by_dimension[map_dim].first_index();
size_type dis_ie1 = first_dis_ie + ie1;
if (dis_ie2 != dis_ie1) return dis_ie2;
return S.master(1);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M) \
template void geo_base_rep<T,M>::init_neighbour() const; \
template geo_base_rep<T,M>::size_type \
geo_base_rep<T,M>::neighbour (size_type, size_type) const;
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
} // namespace
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
rheolef::geo_base_rep< T, distributed >::size_type
base::size_type size_type
Definition: geo.h:533
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::geo_element::master
size_type master(bool i) const
Definition: geo_element.h:165
rheolef::index_set::inplace_intersection
void inplace_intersection(const index_set &b)
Definition: index_set_body.icc:126
rheolef::geo_abstract_base_rep::size_type
geo_element_hack::size_type size_type
Definition: geo.h:260
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::geo_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::index_set
Definition: index_set_header.icc:58
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::geo_element::get_side_orientation
orientation_type get_side_orientation(const geo_element &S) const
Definition: geo_element.cc:240
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
mkgeo_ball.variant
variant
Definition: mkgeo_ball.sh:149
vertex
const point vertex[n_vertex]
Definition: edge.icc:67
rheolef::geo_element::n_node
size_type n_node() const
Definition: geo_element.h:170
rheolef::geo_base_rep::init_neighbour
void init_neighbour() const
Definition: geo_neighbour.cc:83
rheolef::geo_iterator
geo iterator
Definition: geo.h:193
rheolef::geo_element::orientation_type
geo_element_indirect::orientation_type orientation_type
Definition: geo_element.h:134
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::geo_element::set_master
void set_master(bool i, size_type dis_ie) const
Definition: geo_element.h:174
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::space_numbering::dis_inod
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
Definition: space_numbering.cc:177
Float
see the Float page for the full documentation
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::distributor::is_owned
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
Definition: distributor.h:213
rheolef::distributor::first_index
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition: distributor.h:151
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
rheolef::geo_element::face
size_type face(size_type i) const
Definition: geo_element.h:210
rheolef::distributed
distributed
Definition: asr.cc:228
rheolef::geo_element::edge
size_type edge(size_type i) const
Definition: geo_element.h:209
rheolef::geo_base_rep::neighbour
size_type neighbour(size_type ie, size_type loc_isid) const
Definition: geo_neighbour.cc:234