Loading...
Searching...
No Matches
space_constitution_assembly.cc
Go to the documentation of this file.
//
// helper for dof numbering with space product and component access
// continuation of space_constitution.cc
//
// contents:
// 8. dos_idof, for assembly
//
#include "rheolef/geo_domain.h"
#include "rheolef/space.h"
#include "rheolef/space_numbering.h"
#include "rheolef/space_mult.h"
#include "rheolef/piola_util.h"
#include "rheolef/basis_get.h"
#include <sstream>
namespace rheolef {
// ============================================================================
// 8. dis_idof, for assembly
// ============================================================================
template <class T, class M>
typename space_constitution_rep<T,M>::size_type
space_constitution_rep<T,M>::loc_ndof (const reference_element& hat_K) const
{
// lazy evaluated on the fly:
if ( _loc_ndof [hat_K.variant()] != std::numeric_limits<size_type>::max()) {
return _loc_ndof [hat_K.variant()];
}
if (! is_hierarchical()) {
size_type loc_comp_ndof = get_basis().ndof (hat_K);
size_type n_comp = (!_is_hier ? 1 : size());
_loc_ndof [hat_K.variant()] = n_comp*loc_comp_ndof;
return _loc_ndof [hat_K.variant()];
}
// mixed multi-component:
size_type loc_ndof = 0;
for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
const space_constitution<T,M>& constit = *iter;
loc_ndof += constit.loc_ndof(hat_K); // recursive call
}
_loc_ndof [hat_K.variant()] = loc_ndof;
return loc_ndof;
}
// DG:
// integrate ("internal_sides", jump(u)*average(v));
// => in that case, have double dofs on internal sides
template <class T, class M>
static
bool is_dg_not_broken (
const space_constitution_rep<T,M>& constit,
const geo_basic<T,M>& omega_K,
const geo_element& K)
{
bool is_dg = (!constit.get_basis().is_continuous()
&& omega_K.map_dimension() + 1 == constit.get_geo().map_dimension()
&& omega_K.map_dimension() == K.dimension());
// K is a side on omega_K and geo_space is_broken: two domains have the "is_broken" flag: "sides" and "internal_sides"
bool is_broken = constit.get_geo().is_broken(); // e.g. HDG lagrange multipliers
bool res = (is_dg && !is_broken);
return res;
}
// HDG:
// space Mh(omega["sides"],"P1d"); trial lambda(Mh);
// integrate (omega, on_local_sides(lambda));
// => in that case, have to loop on sides
template <class T, class M>
static
bool is_broken (
const space_constitution_rep<T,M>& constit,
const geo_basic<T,M>& omega_K)
{
// only two domains have the "is_broken" flag: "sides" and "internal_sides"
bool res = constit.get_geo().is_broken() && (omega_K.map_dimension() == constit.get_geo().map_dimension() + 1);
return res;
}
// HDG-POST:
// space Mhs(omega,"P1d[sides]"); trial lambda_s(Mhs);
// integrate (omega, on_local_sides(lambda_s));
// => in that case, have to loop on sides with bi-valued lambda_s
template <class T, class M>
static
bool is_dg_restricted_to_interface (
const space_constitution_rep<T,M>& constit,
const geo_basic<T,M>& omega_K,
const geo_element& K)
{
return constit.get_basis().option().is_trace_n()
&& constit.get_geo().get_background_geo().map_dimension() == omega_K.map_dimension()
&& omega_K.map_dimension() == K.dimension() + 1;
}
// assembly loop on geo_elements K on a domain omega_K
// give its corresponding element K1 in space.geo
template <class T, class M>
const geo_element&
assembly2space_geo_element (
const geo_basic<T,M>& space_geo,
const geo_basic<T,M>& omega_K,
const geo_element& K_in)
{
bool is_broken = space_geo.is_broken() && (omega_K.map_dimension() == space_geo.map_dimension() + 1);
bool use_dom2bgd = (omega_K.variant() == geo_abstract_base_rep<T>::geo_domain &&
space_geo.variant() != geo_abstract_base_rep<T>::geo_domain);
bool use_bgd2dom = (omega_K.variant() != geo_abstract_base_rep<T>::geo_domain &&
space_geo.variant() == geo_abstract_base_rep<T>::geo_domain);
if (!is_broken && use_bgd2dom) { // get from background mesh:
return space_geo.bgd2dom_geo_element (K_in);
} else if (!is_broken && use_dom2bgd) { // get on background mesh:
return omega_K.dom2bgd_geo_element (K_in);
} else { // usual case: use K directly
return K_in;
}
}
// assembly_loc_ndof
// assembly_dis_idof : take into account DG on internal sides with doubled loc_ndof
template <class T, class M>
typename space_constitution_rep<T,M>::size_type
space_constitution_rep<T,M>::assembly_loc_ndof (
const geo_basic<T,M>& omega_K,
const geo_element& K_in) const
{
if (! is_hierarchical()) {
trace_macro("loc_ndof: omega_K="<<omega_K.name()<<", K_in="<<K_in.name()<<K_in.dis_ie()<<", space="<<name()<<"...");
const geo_element& K = assembly2space_geo_element (get_geo(), omega_K, K_in);
size_type loc_ndof = 0;
if (is_dg_not_broken (*this, omega_K, K_in)) {
trace_macro("loc_ndof: is_dg_not_broken...");
// dg: check for internal sides with doubled dofs
size_type map_d = K.dimension();
size_type dis_ie0 = K.master(0),
dis_ie1 = K.master(1);
check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
"unexpected isolated side K_in="<<K_in.name()<<K_in.dis_ie() << " in omega_K=" << omega_K.name()
<< " associated to K="<<K.name()<<K.dis_ie() << " in space_geo=" << get_geo().name());
if (dis_ie1 == std::numeric_limits<size_type>::max()) {
// boundary side K
const geo_element& L = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
loc_ndof = this->loc_ndof(L);
} else {
// internal side K
const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
loc_ndof = this->loc_ndof(L0) + this->loc_ndof(L1);
}
} else if (is_dg_restricted_to_interface (*this, omega_K, K_in)) {
trace_macro("loc_ndof: is_dg_restricted_to_interface...");
// HDG-POST multiplier: check for boundary sides dofs
// example: K=e, b="P1d[sides]", omega_K="square" & get_geo="square"
size_type map_d = K.dimension();
size_type dis_ie0 = K.master(0),
dis_ie1 = K.master(1);
check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
"unexpected isolated mesh side K.dis_ie="<<K.dis_ie());
if (dis_ie1 == std::numeric_limits<size_type>::max()) {
// boundary side K
const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
side_information_type sid0;
L0.get_side_informations (K, sid0);
loc_ndof = get_basis().local_ndof_on_side (L0, sid0);
} else {
// internal side K
const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
side_information_type sid0, sid1;
L0.get_side_informations (K, sid0);
L1.get_side_informations (K, sid1);
loc_ndof = get_basis().local_ndof_on_side (L0, sid0)
+ get_basis().local_ndof_on_side (L1, sid1);
// TODO: how to manage e.g. lambda_hs["interface"] = 1; it involves bi-valued side-based space (Pkd(S))^2
fatal_macro ("HDG-POST multipliplier: cannot be extracted on an internal interface (bi-valued)");
}
} else if (is_broken (*this, omega_K)) {
trace_macro("loc_ndof: is_broken...");
// HDG multiplier: check for boundary sides dofs
// example: omega_K="square" & get_geo="square[sides]"
size_type space_map_d = get_geo().map_dimension();
check_macro (space_map_d+1 == K.dimension(), "internal error: unexpected element dimension");
reference_element hat_K = K;
for (size_type isid = 0, nsid = hat_K.n_side(); isid < nsid; ++isid) {
reference_element hat_S = hat_K.side(isid);
loc_ndof += get_basis().ndof(hat_S);
}
} else {
trace_macro("loc_ndof: usual...");
// usual case:
loc_ndof = this->loc_ndof (K);
}
trace_macro("loc_ndof: omega_K="<<omega_K.name()<<", K_in="<<K_in.name()<<K_in.dis_ie()<<", space="<<name()<<": result="<<loc_ndof<<" done");
return loc_ndof;
}
// hybrid multi-component:
trace_macro("loc_ndof_hier: omega_K="<<omega_K.name()<<", K_in="<<K_in.name()<<K_in.dis_ie()<<", space="<<name()<<"...");
check_macro (_is_hier, "unexpected hierarchical space");
size_type loc_ndof = 0;
for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
const space_constitution<T,M>& constit = *iter;
loc_ndof += constit.assembly_loc_ndof(omega_K, K_in); // recursive call
}
trace_macro("loc_ndof_hier: omega_K="<<omega_K.name()<<", K_in="<<K_in.name()<<K_in.dis_ie()<<", space="<<name()<<": result="<<loc_ndof<<" done");
return loc_ndof;
}
template <class T, class M>
void
space_constitution_rep<T,M>::_assembly_dis_idof_recursive (
const geo_basic<T,M>& omega_K,
const geo_element& K_in,
typename std::vector<geo_element::size_type>::iterator& dis_idof_t,
const distributor& hier_ownership,
const std::vector<distributor>& start_by_flattened_component,
size_type& i_flat_comp) const
{
size_type space_map_d = get_geo().map_dimension();
size_type map_d = K_in.dimension();
if (! is_hierarchical()) {
const geo_element& K = assembly2space_geo_element (get_geo(), omega_K, K_in);
size_type loc_flat_comp_ndof = 0;
if (is_dg_not_broken (*this, omega_K, K_in)) {
// DG case:
size_type dis_ie0 = K.master(0),
dis_ie1 = K.master(1);
check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
"unexpected isolated mesh side K.dis_ie="<<K.dis_ie());
if (dis_ie1 == std::numeric_limits<size_type>::max()) {
// boundary side K with one adjacent element
const geo_element& L = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
space_numbering::dis_idof (get_basis(), get_geo().sizes(), L, dis_idof_t);
loc_flat_comp_ndof = get_basis().ndof (L);
} else {
// internal side K with two adjacent elements
const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
space_numbering::dis_idof (get_basis(), get_geo().sizes(), L0, dis_idof_t);
size_type loc_flat_comp_ndof0 = get_basis().ndof (L0);
space_numbering::dis_idof (get_basis(), get_geo().sizes(), L1, dis_idof_t + loc_flat_comp_ndof0);
loc_flat_comp_ndof = loc_flat_comp_ndof0 + get_basis().ndof (L1);
}
} else if (is_dg_restricted_to_interface (*this, omega_K, K_in)) {
// HDG-POST multiplier: check for boundary sides dofs
// example: K=e, b="P1d[sides]", omega_K="square" & get_geo="square"
size_type dis_ie0 = K.master(0),
dis_ie1 = K.master(1);
check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
"unexpected isolated mesh side K.dis_ie="<<K.dis_ie());
if (dis_ie1 == std::numeric_limits<size_type>::max()) {
// boundary side K with one adjacent element
const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
side_information_type sid0;
L0.get_side_informations (K, sid0);
std::vector<size_type> dis_idof_L0;
space_numbering::dis_idof (get_basis(), get_geo().sizes(), L0, dis_idof_L0);
Eigen::Matrix<size_type,Eigen::Dynamic,1> loc_idof_L0;
get_basis().local_idof_on_side (L0, sid0, loc_idof_L0);
for (size_type loc_sid_idof = 0, loc_sid_ndof = loc_idof_L0.size(); loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
size_type loc_idof = loc_idof_L0 [loc_sid_idof];
dis_idof_t [loc_sid_idof] = dis_idof_L0 [loc_idof];
}
loc_flat_comp_ndof = loc_idof_L0.size();
} else {
// internal side K with two adjacent elements
const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
side_information_type sid0, sid1;
L0.get_side_informations (K, sid0);
L1.get_side_informations (K, sid1);
std::vector<size_type> dis_idof_L0, dis_idof_L1;
space_numbering::dis_idof (get_basis(), get_geo().sizes(), L0, dis_idof_L0);
space_numbering::dis_idof (get_basis(), get_geo().sizes(), L1, dis_idof_L1);
Eigen::Matrix<size_type,Eigen::Dynamic,1> loc_idof_L0, loc_idof_L1;
get_basis().local_idof_on_side (L0, sid0, loc_idof_L0);
get_basis().local_idof_on_side (L1, sid1, loc_idof_L1);
check_macro (loc_idof_L0.size() == loc_idof_L1.size(), "unexpected bi-valued side: unmatch sizes");
for (size_type loc_sid_idof = 0, loc_sid_ndof = loc_idof_L0.size(); loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
size_type loc_idof = loc_idof_L0 [loc_sid_idof];
dis_idof_t [loc_sid_idof] = dis_idof_L0 [loc_idof];
dis_idof_t [loc_sid_ndof+loc_sid_idof] = dis_idof_L1 [loc_idof];
}
loc_flat_comp_ndof = loc_idof_L0.size() + loc_idof_L1.size();
fatal_macro ("HDG-POST multipliplier: cannot be extracted on an internal interface (bi-valued)");
}
} else if (is_broken (*this, omega_K)) {
// HDG multiplier: check for boundary sides dofs
// example: omega_K="square" & get_geo="square[sides]"
check_macro (space_map_d+1 == K.dimension(), "internal error: unexpected element dimension");
for (size_type isid = 0, nsid = K_in.n_subgeo(space_map_d); isid < nsid; ++isid) {
size_type dis_isid_in = (K_in.dimension() == 1) ? K_in[isid] : (K_in.dimension() == 2) ? K_in.edge(isid) : K_in.face(isid);
const geo_element& S_in = omega_K.dis_get_geo_element (space_map_d, dis_isid_in);
bool have_bgd_dom = (get_geo().variant() == geo_abstract_base_rep<T>::geo_domain && get_geo().name() != omega_K.name());
const geo_element* S_ptr = (!have_bgd_dom) ? &S_in : &(get_geo().bgd2dom_geo_element(S_in));
const geo_element& S = *S_ptr;
space_numbering::dis_idof (get_basis(), get_geo().sizes(), S, dis_idof_t + loc_flat_comp_ndof);
loc_flat_comp_ndof += get_basis().ndof (S);
}
} else {
// non-broken cases: dg or usual
// usual: non-DG case:
space_numbering::dis_idof (get_basis(), get_geo().sizes(), K, dis_idof_t);
loc_flat_comp_ndof = get_basis().ndof (K);
}
// shift dis_idof depending on i_flat_comp
for (size_type loc_flat_comp_idof = 0; loc_flat_comp_idof < loc_flat_comp_ndof; ++loc_flat_comp_idof) {
size_type dis_flat_comp_idof = dis_idof_t [loc_flat_comp_idof];
distributor comp_ownership = ownership();
size_type iproc = comp_ownership.find_owner (dis_flat_comp_idof);
size_type first_dis_flat_comp_idof = comp_ownership.first_index(iproc);
size_type first_dis_idof = hier_ownership.first_index(iproc);
check_macro (dis_flat_comp_idof >= first_dis_flat_comp_idof, "unexpected dis_comp_idof");
size_type flat_comp_idof = dis_flat_comp_idof - first_dis_flat_comp_idof;
size_type start_flat_comp_idof = start_by_flattened_component [i_flat_comp].size(iproc);
size_type idof = start_flat_comp_idof + flat_comp_idof;
size_type dis_idof = first_dis_idof + idof;
dis_idof_t [loc_flat_comp_idof] = dis_idof;
}
dis_idof_t += loc_flat_comp_ndof;
i_flat_comp++;
return;
}
// hierarchical case
for (size_type i_comp = 0, n_comp = size(); i_comp < n_comp; i_comp++) {
const space_constitution<T,M>& comp_constit = operator[] (i_comp);
comp_constit.data()._assembly_dis_idof_recursive (omega_K, K_in, dis_idof_t, hier_ownership, start_by_flattened_component, i_flat_comp); // recursive call
}
}
template <class T, class M>
void
space_constitution_rep<T,M>::assembly_dis_idof (
const geo_basic<T,M>& omega_K,
const geo_element& K_in,
std::vector<geo_element::size_type>& dis_idof_t) const
{
size_type i_flat_comp = 0;
size_type loc_ndof = assembly_loc_ndof (omega_K, K_in);
dis_idof_t.resize (loc_ndof);
typename std::vector<geo_element::size_type>::iterator first_dis_idof_t = dis_idof_t.begin();
_assembly_dis_idof_recursive (omega_K, K_in, first_dis_idof_t, ownership(), _start_by_flattened_component, i_flat_comp);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class space_constitution_rep<Float,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class space_constitution_rep<Float,distributed>;
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::distributor::find_owner
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
Definition distributor.cc:106
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:158
rheolef::geo_abstract_base_rep
abstract base interface class
Definition geo.h:248
rheolef::geo_abstract_base_rep::name
virtual std::string name() const =0
rheolef::geo_basic
generic mesh with rerefence counting
Definition geo.h:1089
rheolef::geo_element
see the geo_element page for the full documentation
Definition geo_element.h:102
rheolef::geo_element::edge
size_type edge(size_type i) const
Definition geo_element.h:209
rheolef::geo_element::face
size_type face(size_type i) const
Definition geo_element.h:210
rheolef::geo_element::dimension
size_type dimension() const
Definition geo_element.h:167
rheolef::geo_element::master
size_type master(bool i) const
Definition geo_element.h:165
rheolef::geo_element::get_side_informations
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
Definition geo_element.cc:185
rheolef::geo_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition geo_element.h:212
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition geo_element.h:163
rheolef::geo_element::name
char name() const
Definition geo_element.h:169
rheolef::reference_element
see the reference_element page for the full documentation
Definition reference_element.h:66
rheolef::reference_element::side
reference_element side(size_type loc_isid) const
Definition reference_element.h:119
rheolef::reference_element::variant
variant_type variant() const
Definition reference_element.h:99
rheolef::reference_element::n_side
size_type n_side() const
Definition reference_element.h:104
rheolef::smart_pointer_base::data
const T & data() const
Definition smart_pointer.h:266
rheolef::space_constitution_rep
Definition space_constitution.h:208
rheolef::space_constitution_rep::get_basis
const basis_basic< T > & get_basis() const
Definition space_constitution.cc:247
rheolef::space_constitution_rep::assembly_dis_idof
void assembly_dis_idof(const geo_basic< T, M > &dom, const geo_element &bgd_K, std::vector< geo_element::size_type > &dis_idof) const
Definition space_constitution_assembly.cc:344
rheolef::space_constitution_rep::_assembly_dis_idof_recursive
void _assembly_dis_idof_recursive(const geo_basic< T, M > &dom, const geo_element &bgd_K, typename std::vector< geo_element::size_type >::iterator &dis_idof_t, const distributor &hier_ownership, const std::vector< distributor > &start_by_flattened_component, size_type &i_flat_comp) const
Definition space_constitution_assembly.cc:223
rheolef::space_constitution_rep::loc_ndof
size_type loc_ndof(const reference_element &hat_K) const
Definition space_constitution_assembly.cc:43
rheolef::space_constitution_rep::const_iterator
hierarchy_type::const_iterator const_iterator
Definition space_constitution.h:216
rheolef::space_constitution_rep::size_type
hierarchy_type::size_type size_type
Definition space_constitution.h:214
rheolef::space_constitution_rep::get_geo
const geo_basic< T, M > & get_geo() const
Definition space_constitution.cc:209
rheolef::space_constitution_rep::assembly_loc_ndof
size_type assembly_loc_ndof(const geo_basic< T, M > &dom, const geo_element &bgd_K) const
Definition space_constitution_assembly.cc:137
rheolef::space_constitution
Definition space_constitution.h:405
rheolef::space_constitution::loc_ndof
size_type loc_ndof(const reference_element &hat_K) const
Definition space_constitution.h:443
rheolef::space_constitution::assembly_loc_ndof
size_type assembly_loc_ndof(const geo_basic< T, M > &dom, const geo_element &bgd_K) const
Definition space_constitution.h:444
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
fatal_macro
#define fatal_macro(message)
Definition dis_macros.h:33
get_geo
void get_geo(istream &in, my_geo &omega)
Definition field2gmsh_pos.cc:126
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::space_numbering::dis_idof
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
Definition space_numbering.cc:187
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::assembly2space_geo_element
const geo_element & assembly2space_geo_element(const geo_basic< T, M > &space_geo, const geo_basic< T, M > &omega_K, const geo_element &K_in)
Definition space_constitution_assembly.cc:115
rheolef::side_information_type
Definition reference_element_face_transformation.h:37