geo_subdivide.cc
Go to the documentation of this file.
//
// build a new mesh by subdividing each mesh edge in "nsub" sub-edges
//
// TODO:
// * add qTPH elts
// * add domains
// * distributed: initialize _inod2ios_dis_inod all all geo_rep<dis> data
#include "rheolef/geo.h"
#include "rheolef/space.h"
namespace rheolef {
// ----------------------------------------------------------------------------
// subdivide one element
// ----------------------------------------------------------------------------
template <class M>
static
void
build_edge (
const std::vector<size_t>& new_dis_inod,
size_t nsub,
hack_array<geo_element_hack,M>& ge,
size_t& ie)
{
typedef typename hack_array<geo_element_hack,M>::size_type size_type;
typedef point_basic<size_type> ilat;
for (size_type i = 0; i < nsub; i++) {
size_type loc_inod0 = reference_element_e::ilat2loc_inod (nsub, ilat(i));
size_type loc_inod1 = reference_element_e::ilat2loc_inod (nsub, ilat(i+1));
ge[ie][0] = new_dis_inod[loc_inod0];
ge[ie][1] = new_dis_inod[loc_inod1];
++ie;
}
}
template <class M>
static
void
build_triangle (
const std::vector<size_t>& new_dis_inod,
size_t nsub,
std::array<hack_array<geo_element_hack,M>,reference_element::max_variant>&
ge,
std::array<size_t, reference_element::max_variant>&
count_by_variant)
{
constexpr size_t e = reference_element::e;
constexpr size_t t = reference_element::t;
size_t& ie = count_by_variant[t];
size_t& e_ie = count_by_variant[e];
typedef typename hack_array<geo_element_hack,M>::size_type size_type;
typedef point_basic<size_type> ilat;
for (size_type i = 0; i < nsub; i++) {
for (size_type j = 0; i+j < nsub; j++) {
// add two newly created triangles
size_type loc_inod00 = reference_element_t::ilat2loc_inod (nsub, ilat(i, j));
size_type loc_inod10 = reference_element_t::ilat2loc_inod (nsub, ilat(i+1, j));
size_type loc_inod01 = reference_element_t::ilat2loc_inod (nsub, ilat(i, j+1));
ge[t][ie][0] = new_dis_inod[loc_inod00];
ge[t][ie][1] = new_dis_inod[loc_inod10];
ge[t][ie][2] = new_dis_inod[loc_inod01];
++ie;
if (i+j+1 >= nsub) continue;
size_type loc_inod11 = reference_element_t::ilat2loc_inod (nsub, ilat(i+1, j+1));
ge[t][ie][0] = new_dis_inod[loc_inod10];
ge[t][ie][1] = new_dis_inod[loc_inod11];
ge[t][ie][2] = new_dis_inod[loc_inod01];
++ie;
// add also newly createe internal edges
ge[e][e_ie][0] = new_dis_inod[loc_inod10];
ge[e][e_ie][1] = new_dis_inod[loc_inod11];
++e_ie;
ge[e][e_ie][0] = new_dis_inod[loc_inod11];
ge[e][e_ie][1] = new_dis_inod[loc_inod01];
++e_ie;
ge[e][e_ie][0] = new_dis_inod[loc_inod01];
ge[e][e_ie][1] = new_dis_inod[loc_inod10];
++e_ie;
}
}
}
template <class M>
static
void
build_element (
const reference_element& hat_K,
const std::vector<size_t>& new_dis_inod,
size_t nsub,
std::array<hack_array<geo_element_hack,M>,reference_element::max_variant>&
ge,
std::array<size_t, reference_element::max_variant>&
count_by_variant)
{
switch (hat_K.variant()) {
case reference_element::e:
build_edge (new_dis_inod, nsub, ge[reference_element::e], count_by_variant[reference_element::e]);
break;
case reference_element::t:
build_triangle (new_dis_inod, nsub, ge, count_by_variant);
break;
default:
error_macro ("unexpected element variant");
}
}
// ----------------------------------------------------------------------------
// mesh subdivide
// ----------------------------------------------------------------------------
template <class T, class M>
void
geo_build_by_subdividing (
geo_rep<T,M>& new_omega,
const geo_basic<T,M>& old_omega,
typename geo_rep<T,M>::size_type nsub)
{
typedef typename geo_rep<T,M>::size_type size_type;
// ------------------------------------------------------------------------
// 1) store header infos in geo
// ------------------------------------------------------------------------
new_omega._version = 4;
new_omega._have_connectivity = true;
new_omega._name = old_omega.name() + "-P" + itos(nsub);
new_omega._dimension = old_omega.dimension();
new_omega._gs._map_dimension = old_omega.map_dimension();
new_omega._sys_coord = old_omega.coordinate_system();
new_omega._serial_number = 0;
new_omega._piola_basis = basis_basic<T>("P1");
// P1: subdivided mesh has order=1 for simplicity,
// even when original old_omega is curved, i.e. Pl, l>=1
// ------------------------------------------------------------------------
// 2) count & build nodes
// ------------------------------------------------------------------------
space_basic<T,M> old_Xh_sub (old_omega, "P"+itos(nsub));
new_omega._node = old_Xh_sub.get_xdofs();
new_omega._gs.ownership_by_variant[0]
= new_omega._gs.ownership_by_dimension[0]
= new_omega._gs.node_ownership
= new_omega._node.ownership();
// ------------------------------------------------------------------------
// 3) count elements
// ------------------------------------------------------------------------
size_type map_d = old_omega.map_dimension();
communicator comm = old_omega.comm();
std::array<size_type, reference_element::max_variant> size_by_variant;
std::fill (size_by_variant.begin(), size_by_variant.end(), 0);
for (size_type d = 1; d <= map_d; ++d) {
for (size_type variant = reference_element::first_variant_by_dimension(d);
variant < reference_element:: last_variant_by_dimension(d); variant++) {
size_type loc_nge = pow(nsub,d);
size_type old_nge = old_omega.sizes().ownership_by_variant [variant].size();
size_by_variant[variant] = loc_nge*old_nge;
if (variant == reference_element::t) {
// count also newly locally created internal edges
size_type loc_nedg = 3*nsub*(nsub-1)/2;
size_by_variant[reference_element::e] += loc_nedg*old_nge;
}
// TODO: add internal also for q, T, P, H
}
}
for (size_type d = 1; d <= map_d; ++d) {
size_type ne = 0, dis_ne = 0;
for (size_type variant = reference_element::first_variant_by_dimension(d);
variant < reference_element:: last_variant_by_dimension(d); variant++) {
geo_element::parameter_type param (variant, 1);
new_omega._gs.ownership_by_variant [variant].resize (distributor::decide, comm, size_by_variant[variant]);
new_omega._geo_element [variant].resize (new_omega._gs.ownership_by_variant [variant], param);
ne += size_by_variant[variant];
}
new_omega._gs.ownership_by_dimension [d].resize (distributor::decide, comm, ne);
}
// ------------------------------------------------------------------------
// 4) compute elements
// ------------------------------------------------------------------------
std::array<size_type, reference_element::max_variant> count_by_variant;
std::fill (count_by_variant.begin(), count_by_variant.end(), 0);
std::vector<size_type> dis_inod;
for (size_type d = 1; d <= map_d; ++d) {
for (size_type ie = 0, ne = old_omega.size(d); ie < ne; ++ie) {
const geo_element& old_K = old_omega.get_geo_element (d, ie);
old_Xh_sub.dis_idof (old_K, dis_inod);
size_type variant = old_K.variant();
build_element (old_K, dis_inod, nsub, new_omega._geo_element, count_by_variant);
}
}
// ------------------------------------------------------------------------
// 5) create vertex-element (0d elements)
// ------------------------------------------------------------------------
geo_element::parameter_type param (reference_element::p, 1);
new_omega._geo_element [reference_element::p].resize (new_omega._gs.ownership_by_dimension[0], param);
size_type first_iv = new_omega._gs.ownership_by_dimension[0].first_index();
{
size_type iv = 0;
for (auto iter = new_omega.begin(0), last = new_omega.end(0); iter != last; iter++, iv++) {
geo_element& P = *iter;
P[0] = first_iv + iv;
P.set_dis_ie (first_iv + iv); // TODO: P[0] & dis_ie redundant for `p'
P.set_ios_dis_ie (first_iv + iv);
}
}
// ------------------------------------------------------------------------
// 5b) set faces & edge in new_omega
// ------------------------------------------------------------------------
if (new_omega._gs._map_dimension > 0) {
for (size_type side_dim = new_omega._gs._map_dimension - 1; side_dim >= 1; side_dim--) {
#ifdef TO_CLEAN
size_type nsid = new_omega._gs.ownership_by_dimension[side_dim].dis_size();
new_omega._gs.ownership_by_dimension [side_dim] = distributor (nsid, new_omega.comm(), nsid);
#endif // TO_CLEAN
size_type isid = 0;
for (typename geo_rep<T,M>::iterator iter = new_omega.begin(side_dim), last = new_omega.end(side_dim); iter != last; iter++, isid++) {
geo_element& S = *iter;
S.set_ios_dis_ie (isid);
S.set_dis_ie (isid);
}
}
}
// ------------------------------------------------------------------------
// 6) K.set_dis_ie
// ------------------------------------------------------------------------
for (size_type d = 0; d <= map_d; ++d) {
size_type first_dis_ie = new_omega._gs.ownership_by_dimension[d].first_index();
for (size_type ie = 0, ne = new_omega.size(d); ie < ne; ++ie) {
geo_element& K = new_omega.get_geo_element(d,ie);
K.set_dis_ie (first_dis_ie + ie);
K.set_ios_dis_ie (first_dis_ie + ie);
}
}
#ifdef TODO
// ------------------------------------------------------------------------
// 8) get domain, until end-of-file
// ------------------------------------------------------------------------
// TODO: build domin in new_omega from those in old_omega
vector<index_set> ball [4];
domain_indirect_basic<sequential> dom;
while (dom.get (ips, *this, ball)) {
base::_domains.push_back (dom);
}
#endif // TODO
// ------------------------------------------------------------------------
// 9) set indexes on faces and edges of elements, for P2 approx
// ------------------------------------------------------------------------
new_omega.set_element_side_index (1);
new_omega.set_element_side_index (2);
// ------------------------------------------------------------------------
// 10) bounding box: _xmin, _xmax
// ------------------------------------------------------------------------
new_omega.compute_bbox();
}
#define _RHEOLEF_geo_build_by_subdividing(M) \
template <class T> \
void \
geo_rep<T,M>::build_by_subdividing (const geo_basic<T,M>& omega, size_type nsub) \
{ \
geo_build_by_subdividing (*this, omega, nsub); \
}
_RHEOLEF_geo_build_by_subdividing(sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_geo_build_by_subdividing(distributed)
#endif // _RHEOLEF_HAVE_MPI
#undef _RHEOLEF_geo_build_by_subdividing
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciate(T,M) \
template void geo_rep<Float,M>::build_by_subdividing ( \
const geo_basic<Float,M>& omega, size_t nsub);
_RHEOLEF_instanciate(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciate(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
#undef _RHEOLEF_instanciate
} // namespace rheolef
size_type
field::size_type size_type
Definition: branch.cc:425
Float
see the Float page for the full documentation
rheolef::basis_basic
Definition: basis.h:532
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_basic< sequential >
Definition: domain_indirect.h:342
rheolef::domain_indirect_basic< sequential >::get
idiststream & get(idiststream &ips, const geo_rep< T, sequential > &omega, std::vector< index_set > *ball)
Definition: domain_indirect.h:535
rheolef::geo_basic< T, M >
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_element::set_ios_dis_ie
void set_ios_dis_ie(size_type ios_dis_ie)
Definition: geo_element.h:173
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::geo_element::set_dis_ie
void set_dis_ie(size_type dis_ie)
Definition: geo_element.h:172
rheolef::geo_rep
sequential mesh representation
Definition: geo.h:778
rheolef::hack_array::size_type
std::vector< T, A >::size_type size_type
Definition: hack_array.h:349
rheolef::reference_element_e::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:350
rheolef::reference_element_t::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:413
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element::max_variant
static const variant_type max_variant
Definition: reference_element.h:82
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::space_basic< T, M >
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
mkgeo_ball.variant
int variant
Definition: mkgeo_ball.sh:149
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
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
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::_RHEOLEF_instanciate
_RHEOLEF_instanciate(Float, sequential) _RHEOLEF_instanciate(Float
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
rheolef::_RHEOLEF_geo_build_by_subdividing
_RHEOLEF_geo_build_by_subdividing(sequential) _RHEOLEF_geo_build_by_subdividing(distributed) _RHEOLEF_instanciate(Float
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
rheolef::geo_build_by_subdividing
void geo_build_by_subdividing(geo_rep< T, M > &new_omega, const geo_basic< T, M > &old_omega, typename geo_rep< T, M >::size_type k)
Definition: geo_subdivide.cc:129
rheolef::distributed
Definition: distributed.h:88
rheolef::geo_element::parameter_type
Definition: geo_element.h:136
rheolef::sequential
Definition: distributed.h:28