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" + std::to_string(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"+std::to_string(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