MPI_HGeometry.h

浏览该文件的文档。

#ifndef __MPI_HGeometry_h__
#define __MPI_HGeometry_h__

#include <set>
#include <AFEPack/HGeometry.h>
#include "MPI.h"

namespace MPI {

  template <int DIM, int DOW, class MATCHER> class HGeometryForest;
  template <class FOREST> class BirdView;
  template <class FOREST> class BirdViewSet;
  template <class FOREST> class HGeometryMatcher;
  template <class FOREST> class HLoadBalance;

  template <int DOW>
    struct PointDistance {
      int value(const Point<DOW>& p0,
                const Point<DOW>& p1,
                double tol) const {
        return (distance(p0, p1) < tol)?0:-1;
      }
    };

  template <int DIM, int DOW=DIM, class MATCHER=PointDistance<DOW> >
    class HGeometryForest : public HGeometryTree<DIM,DOW> {
  private:
  typedef HGeometryForest<DIM,DOW,MATCHER> this_t;
  public:
  typedef HGeometryTree<DIM,DOW> tree_t;
  typedef MATCHER matcher_t;

  template <int D>
  const Shared_ptr_list<HGeometry<D,DOW> >&
  get_shared_list() const {
    const void * p_list;
    switch (D) {
    case 0: p_list = &_shared_list_0d; break;
    case 1: p_list = &_shared_list_1d; break;
    case 2: p_list = &_shared_list_2d; break;
    case 3: p_list = &_shared_list_3d; break;
    }
    return *((const Shared_ptr_list<HGeometry<D,DOW> > *)p_list);
  }

  template <int GDIM>
  bool is_on_primary_rank(const HGeometry<GDIM,DOW>& geo) const {
    typedef HGeometry<GDIM,DOW> geo_t;
    typedef Shared_object<geo_t> obj_t;
    obj_t * p_obj = this->get_shared_info(geo);
    if (p_obj == NULL) {
      return true;
    } else {
      return p_obj->is_on_primary_rank(_rank);
    }
  }

  template <int GDIM>
  int primary_rank(const HGeometry<GDIM,DOW>& geo) const {
    typedef HGeometry<GDIM,DOW> geo_t;
    typedef Shared_object<geo_t> obj_t;
    obj_t * p_obj = this->get_shared_info(geo);
    if (p_obj == NULL) {
      return _rank;
    } else {
      return p_obj->primary_rank(_rank);
    }
  }

  template <int GDIM>
  bool is_primary_geometry(const HGeometry<GDIM,DOW>& geo) const {
    typedef HGeometry<GDIM,DOW> geo_t;
    typedef Shared_object<geo_t> obj_t;
    obj_t * p_obj = this->get_shared_info(geo);
    if (p_obj == NULL) {
      return true;
    } else {
      return p_obj->is_primary_object(_rank);
    }
  }

  bool is_geometry_shared(const RegularMesh<DIM,DOW>& mesh,
                          int dim, int idx) const {
    switch(dim) {
    case 0: return (get_shared_info(*mesh.template h_geometry<0>(idx)) != NULL);
    case 1: return (get_shared_info(*mesh.template h_geometry<1>(idx)) != NULL);
    case 2: return (get_shared_info(*mesh.template h_geometry<2>(idx)) != NULL);
    case 3: return (get_shared_info(*mesh.template h_geometry<3>(idx)) != NULL);
    }
  }

  const matcher_t& matcher() const { return _matcher; }
  matcher_t& matcher() { return _matcher; }

  private:
  MPI_Comm _comm;
  int _rank;
  int _n_rank;
  matcher_t _matcher;

  property_id_t<> _pid_shared_info_sent; 
  template <class GEO> bool
  is_shared_info_sent(GEO& geo) const {
    return (geo.get_property(_pid_shared_info_sent) != NULL);
  }
  template <class GEO> void
  set_shared_info_sent(GEO& geo) const {
    geo.new_property(_pid_shared_info_sent);
  }
  template <class GEO> void
  clear_shared_info_sent(GEO& geo) const {
    geo.free_property(_pid_shared_info_sent);
  }

  property_id_t<> _pid_dummy;
  template <class GEO> bool
  is_dummy(GEO& geo) const {
    return (geo.get_property(_pid_dummy) != NULL);
  }
  template <class GEO> void
  dummy(GEO& geo) const {
    geo.new_property(_pid_dummy);
  }
  template <class GEO> void
  undummy(GEO& geo) const {
    geo.free_property(_pid_dummy);
  }

#define GDIM 0
#define GEO HGeometry<GDIM,DOW>
#define OBJ Shared_object<GEO>
  private:
  property_id_t<OBJ> _pid_so_0d;
  Shared_ptr_list<GEO> _shared_list_0d;

  public:
  OBJ * get_shared_info(const GEO& geo) const {
    return geo.get_property(_pid_so_0d);
  }
  OBJ * new_shared_info(GEO& geo) {
    OBJ * p_obj = geo.new_property(_pid_so_0d);
    p_obj->local_pointer() = &geo;
    _shared_list_0d.push_back(p_obj);
    return p_obj;
  }
  void erase_shared_info(GEO& geo) {
    OBJ * p_obj = geo.get_property(_pid_so_0d);
    if (p_obj == NULL) return;
    _shared_list_0d.erase(std::find(_shared_list_0d.begin_ptr(), 
                                    _shared_list_0d.end_ptr(), 
                                    p_obj));
  }
#undef OBJ
#undef GEO
#undef GDIM

#define GDIM 1
#define GEO HGeometry<GDIM,DOW>
#define OBJ Shared_object<GEO>
  private:
  property_id_t<OBJ> _pid_so_1d;
  Shared_ptr_list<GEO> _shared_list_1d;

  public:
  OBJ * get_shared_info(const GEO& geo) const {
    return geo.get_property(_pid_so_1d);
  }
  OBJ * new_shared_info(GEO& geo) {
    OBJ * p_obj = geo.new_property(_pid_so_1d);
    p_obj->local_pointer() = &geo;
    _shared_list_1d.push_back(p_obj);
    return p_obj;
  }
  void erase_shared_info(GEO& geo) {
    OBJ * p_obj = geo.get_property(_pid_so_1d);
    if (p_obj == NULL) return;
    _shared_list_1d.erase(std::find(_shared_list_1d.begin_ptr(), 
                                    _shared_list_1d.end_ptr(), 
                                    p_obj));
  }
#undef OBJ
#undef GEO
#undef GDIM

#define GDIM 2
#define GEO HGeometry<GDIM,DOW>
#define OBJ Shared_object<GEO>
  private:
  property_id_t<OBJ> _pid_so_2d;
  Shared_ptr_list<GEO> _shared_list_2d;

  public:
  OBJ * get_shared_info(const GEO& geo) const {
    return geo.get_property(_pid_so_2d);
  }
  OBJ * new_shared_info(GEO& geo) {
    OBJ * p_obj = geo.new_property(_pid_so_2d);
    p_obj->local_pointer() = &geo;
    _shared_list_2d.push_back(p_obj);
    return p_obj;
  }
  void erase_shared_info(GEO& geo) {
    OBJ * p_obj = geo.get_property(_pid_so_2d);
    if (p_obj == NULL) return;
    _shared_list_2d.erase(std::find(_shared_list_2d.begin_ptr(), 
                                    _shared_list_2d.end_ptr(), 
                                    p_obj));
  }
#undef OBJ
#undef GEO
#undef GDIM

#define GDIM 3
#define GEO HGeometry<GDIM,DOW>
#define OBJ Shared_object<GEO>
  private:
  property_id_t<OBJ> _pid_so_3d;
  Shared_ptr_list<GEO> _shared_list_3d;

  public:
  OBJ * get_shared_info(const GEO& geo) const {
    return geo.get_property(_pid_so_3d);
  }
  OBJ * new_shared_info(GEO& geo) {
    OBJ * p_obj = geo.new_property(_pid_so_3d);
    p_obj->local_pointer() = &geo;
    _shared_list_3d.push_back(p_obj);
    return p_obj;
  }
  void erase_shared_info(GEO& geo) {
    OBJ * p_obj = geo.get_property(_pid_so_3d);
    if (p_obj == NULL) return;
    _shared_list_3d.erase(std::find(_shared_list_3d.begin_ptr(), 
                                    _shared_list_3d.end_ptr(), 
                                    p_obj));
  }
#undef OBJ
#undef GEO
#undef GDIM

  bool lock() { return HGeometryTree<DIM,DOW>::lock(); }
  void unlock() { HGeometryTree<DIM,DOW>::unlock(); }

  public:
  HGeometryForest() { new_property(); }
  HGeometryForest(MPI_Comm comm) {
    set_communicator(comm);
    new_property();
  }
  virtual ~HGeometryForest() { this->clear(); }

  void clear() {
    _shared_list_0d.clear();
    _shared_list_1d.clear();
    _shared_list_2d.clear();
    _shared_list_3d.clear();

    free_property();
    HGeometryTree<DIM,DOW>::clear();
    new_property();
  }

  void new_property() {
    new_property_id(_pid_shared_info_sent);
    new_property_id(_pid_dummy);
    if (DIM >= 0) new_property_id(_pid_so_0d);
    if (DIM >= 1) new_property_id(_pid_so_1d);
    if (DIM >= 2) new_property_id(_pid_so_2d);
    if (DIM >= 3) new_property_id(_pid_so_3d);
  }

  void free_property() {
    free_property_id(_pid_shared_info_sent);
    free_property_id(_pid_dummy);
    if (DIM >= 0) free_property_id(_pid_so_0d);
    if (DIM >= 1) free_property_id(_pid_so_1d);
    if (DIM >= 2) free_property_id(_pid_so_2d);
    if (DIM >= 3) free_property_id(_pid_so_3d);
  }

  void set_communicator(MPI_Comm comm) {
    _comm = comm;
    MPI_Comm_rank(_comm, &_rank);
    MPI_Comm_size(_comm, &_n_rank);
  }
  MPI_Comm communicator() const { return _comm; }
  int rank() const { return _rank; }
  int n_rank() const { return _n_rank; }
  void readMesh(const std::string&);

  void eraseRootElement(u_int level = 1);
  void renumerateRootElement(void (*f)(double, double, double,
                                       double&,double&,double&) = NULL);

  private:
  void eraseRootElementOneLevel();
  template <class GEO>
  void nullParent(GEO& geo) {
    geo.parent = NULL;
    for (u_int i = 0;i < GEO::n_boundary;++ i) {
      this->nullParent(*geo.boundary[i]);
    }
  }
  template <class GEO>
  void tryDeleteGeometry(GEO * p_geo, 
                         const property_id_t<bool>& pid) {
    for (u_int i = 0;i < GEO::n_boundary;++ i) {
      if (p_geo->boundary[i]->get_property(pid) != NULL) {
        p_geo->boundary[i] = NULL;
      } else {
        this->tryDeleteGeometry(p_geo->boundary[i], pid);
      }
    }
    bool * p_prp = p_geo->get_property(pid);
    assert (p_prp == NULL);
    p_geo->new_property(pid);
    this->erase_shared_info(*p_geo); 
  }
  template <class GEO>
  void deleteGeometry(GEO * p_geo) {
    for (u_int i = 0;i < GEO::n_boundary;++ i) {
      if (p_geo->boundary[i] != NULL) {
        this->deleteGeometry(p_geo->boundary[i]);
      }
    }
    delete p_geo;
  }

  friend class BirdView<this_t>;
  friend class HGeometryMatcher<this_t>;
  friend class HLoadBalance<this_t>;

  public:
  typedef BirdView<this_t> birdview_t;
  typedef BirdViewSet<this_t> birdview_set_t;
  typedef HLoadBalance<this_t> load_balancer_t;

  }; // HGeometryForest

  template <class FOREST>
    class BirdView : public IrregularMesh<FOREST::dim,FOREST::dow> {
  public:
    enum {dim = FOREST::dim, dow = FOREST::dow};
    typedef IrregularMesh<dim, dow> base_t;
    typedef FOREST forest_t;
  public:
  BirdView() : base_t() {}
    explicit BirdView(forest_t& forest) : base_t(forest) {}
    virtual ~BirdView() {}

    forest_t& getForest() const {
      return dynamic_cast<forest_t&>(this->geometryTree());
    }
    virtual void semiregularize();
    void eraseRootElement(u_int level = 1);
  private:
    void eraseRootElementOneLevel();
  };

  template <class FOREST>
    class BirdViewSet : public std::set<BirdView<FOREST>*> {
  private:
    typedef BirdView<FOREST> birdview_t;
    typedef std::set<birdview_t*> base_t;
  public:
    typedef _Deref_iterator<typename base_t::iterator, birdview_t> iterator;
    typedef _Deref_iterator<typename base_t::const_iterator, const birdview_t> const_iterator;
  public:
    void add(birdview_t& bv) { base_t::insert(&bv); }
    void erase(birdview_t& bv) { base_t::erase(&bv); }

    iterator begin() { return base_t::begin(); }
    iterator end() { return base_t::end(); }

    const_iterator begin() const { return base_t::begin(); }
    const_iterator end() const { return base_t::end(); }

    typename base_t::iterator begin_ptr() { return base_t::begin(); }
    typename base_t::iterator end_ptr() { return base_t::end(); }

    typename base_t::const_iterator begin_ptr() const { return base_t::begin(); }
    typename base_t::const_iterator end_ptr() const { return base_t::end(); }
  };

  template <class FOREST>
    BirdViewSet<FOREST> 
    make_set(BirdView<FOREST>& mesh0, 
             BirdView<FOREST>& mesh1) {
    BirdViewSet<FOREST> bvs;
    bvs.add(mesh0);
    bvs.add(mesh1);
    return bvs;
  }
  template <class FOREST>
    BirdViewSet<FOREST> 
    make_set(BirdView<FOREST>& mesh0, 
             BirdView<FOREST>& mesh1,
             BirdView<FOREST>& mesh2) {
    BirdViewSet<FOREST> bvs;
    bvs.add(mesh0);
    bvs.add(mesh1);
    bvs.add(mesh2);
    return bvs;
  }
  template <class FOREST>
    BirdViewSet<FOREST> 
    make_set(BirdView<FOREST>& mesh0, 
             BirdView<FOREST>& mesh1,
             BirdView<FOREST>& mesh2,
             BirdView<FOREST>& mesh3) {
    BirdViewSet<FOREST> bvs;
    bvs.add(mesh0);
    bvs.add(mesh1);
    bvs.add(mesh2);
    bvs.add(mesh3);
    return bvs;
  }

  template <class FOREST>
    class HGeometryMatcher {
  public:
    enum {dim = FOREST::dim, dow = FOREST::dow};
    typedef FOREST forest_t;
  private:
    typedef HGeometryMatcher<forest_t> this_t;
    typedef typename forest_t::matcher_t matcher_t;

    HTools _tools;
    forest_t * _forest;
    bool _is_operated;

    const matcher_t& matcher() const { return _forest->matcher(); }

    template <class GEO> bool 
      is_shared_info_sent(GEO& geo) const {
      return _forest->is_shared_info_sent(geo);
    }
    template <class GEO> void 
      set_shared_info_sent(GEO& geo) const {
      _forest->set_shared_info_sent(geo);
    }
    template <class GEO> Shared_object<GEO> * 
      new_shared_info(GEO& geo) const {
      return _forest->new_shared_info(geo);
    }
    template <class GEO> Shared_object<GEO> *
      get_shared_info(GEO& geo) const {
      return _forest->get_shared_info(geo);
    }

  public:
  HGeometryMatcher(forest_t& forest) : _forest(&forest) {}

    bool match_geometry();
    template <int GDIM> bool 
      is_pack_match_geometry(HGeometry<GDIM,dow> * geo) const;
    template <int GDIM> void 
      pack_match_geometry(HGeometry<GDIM,dow> * geo,
                          int remote_rank,
                          AFEPack::ostream<>& os);
    template <int GDIM> void 
      unpack_match_geometry(HGeometry<GDIM,dow> * geo,
                            int remote_rank,
                            AFEPack::istream<>& is);

  public:
    void set_matcher(const matcher_t& _matcher) { matcher = _matcher; }

    bool sync_is_used();
    template <int GDIM> void 
      pack_sync_is_used(HGeometry<GDIM,dow> * geo,
                        int remote_rank,
                        AFEPack::ostream<>& os);
    template <int GDIM> void 
      unpack_sync_is_used(HGeometry<GDIM,dow> * geo,
                          int remote_rank,
                          AFEPack::istream<>& is);

  private:
    template <class GEO>
      void geometry_reference_point(GEO& geo, 
                                    Point<dow>& pnt) const {
      int n_vtx = geo.n_vertex;
      for (int n = 0;n < dow;++ n) {
        pnt[n] = 0.0;
        for (int i = 0;i < n_vtx;++ i) {
          pnt[n] += (*geo.vertex[i])[n];
        }
        pnt[n] /= n_vtx;
      }
    }

  };


#include "MPI_HGeometry.templates.h"

} // namespace MPI

#endif // __MPI_HGeometry_h__