MPI_HGeometry.templates.h

浏览该文件的文档。

#define TEMPLATE template <int DIM, int DOW, class MATCHER>
#define THIS HGeometryForest<DIM,DOW,MATCHER>

TEMPLATE void 
THIS::readMesh(const std::string& filename) {

  const int& rank = this->_rank;
  const int& n_rank = this->_n_rank;

  char filename_buf[256];
  if (n_rank > 1) {
    sprintf(filename_buf, "%s%d.mesh", filename.c_str(), rank);
  } else {
    sprintf(filename_buf, "%s.mesh", filename.c_str());
  }
  std::cerr << "Reading in mesh data file " 
            << filename_buf 
            << " as geometry tree root ..." 
            << std::endl;
  std::ifstream is(filename_buf);

  u_int n_point;
  is >> n_point;
  std::cerr << "\t# points: " << n_point << std::endl;
  std::vector<Point<this_t::dow> > point(n_point);
  for (u_int i = 0;i < n_point;i ++) is >> point[i];
  is >> n_point;
  std::vector<HGeometry<0,this_t::dow> *> geo_0d(n_point, (HGeometry<0,this_t::dow> *)NULL);
  for (u_int i = 0;i < n_point;i ++) {
    u_int j, k;
    is >> j; geo_0d[j] = new HGeometry<0,this_t::dow>();
    is >> k >> k; *dynamic_cast<Point<this_t::dow> *>(geo_0d[j]) = point[k];
    is >> k >> k >> geo_0d[j]->bmark;
  }
  point.clear();
        
#define GDIM 1
  u_int n_geo_1d;
  std::vector<HGeometry<GDIM,this_t::dow>*> geo_1d;
  if (DIM >= 1) {
    is >> n_geo_1d;
    std::cerr << "\t# 1D-geometry: " << n_geo_1d << std::endl;
    geo_1d.resize(n_geo_1d, NULL);
    for (u_int i = 0;i < n_geo_1d;i ++) {
      u_int j, k, l;
      is >> j >> k; geo_1d[j] = new HGeometry<GDIM,this_t::dow>();
      for (k = 0;k < GDIM + 1;k ++) {
        is >> l; geo_1d[j]->vertex[k] = geo_0d[l];
      }
      is >> k;
      for (k = 0;k < GDIM + 1;k ++) {
        is >> l;
      }
      is >> geo_1d[j]->bmark;
    }
  }
#undef GDIM

#define GDIM 2
  u_int n_geo_2d;
  std::vector<HGeometry<GDIM,this_t::dow>*> geo_2d;
  if (DIM >= 2) {
    is >> n_geo_2d;
    std::cerr << "\t# 2D-geometry: " << n_geo_2d << std::endl;
    geo_2d.resize(n_geo_2d, NULL);
    for (u_int i = 0;i < n_geo_2d;i ++) {
      u_int j, k, l;
      is >> j >> k; geo_2d[j] = new HGeometry<GDIM,this_t::dow>();
      for (k = 0;k < GDIM + 1;k ++) {
        is >> l; geo_2d[j]->vertex[k] = geo_0d[l];
      }
      is >> k;
      for (k = 0;k < GDIM + 1;k ++) {
        is >> l; geo_2d[j]->boundary[k] = geo_1d[l];
      }
      is >> geo_2d[j]->bmark;
    }
  }
#undef GDIM

#define GDIM 3
  u_int n_geo_3d;
  std::vector<HGeometry<GDIM,this_t::dow>*> geo_3d;
  if (DIM >= 3) { 
    is >> n_geo_3d;
    std::cerr << "\t# 3D-geometry: " << n_geo_3d << std::endl;
    geo_3d.resize(n_geo_3d, NULL);
    for (u_int i = 0;i < n_geo_3d;i ++) {
      u_int j, k, l;
      is >> j >> k; geo_3d[j] = new HGeometry<GDIM,this_t::dow>();
      for (k = 0;k < GDIM + 1;k ++) {
        is >> l; geo_3d[j]->vertex[k] = geo_0d[l];
      }
      is >> k;
      for (k = 0;k < GDIM + 1;k ++) {
        is >> l; geo_3d[j]->boundary[k] = geo_2d[l];
      }
      is >> geo_3d[j]->bmark;
    }
#undef GDIM
  }
  is.close();

  if (DIM == 1) {
    for (u_int i = 0;i < n_geo_1d;++ i) {
      this->rootElement().push_back((HGeometry<DIM,this_t::dow> *)geo_1d[i]);
    }
  } else if (DIM == 2) {
    for (u_int i = 0;i < n_geo_2d;++ i) {
      this->rootElement().push_back((HGeometry<DIM,this_t::dow> *)geo_2d[i]);
    }
  } else if (DIM == 3) {
    for (u_int i = 0;i < n_geo_3d;++ i) {
      this->rootElement().push_back((HGeometry<DIM,this_t::dow> *)geo_3d[i]);
    }
  }

  if (n_rank == 1) return; 

  sprintf(filename_buf, "%s%d.share", filename.c_str(), rank);
  is.open(filename_buf); 
  std::cerr << "Reading in shared data file " 
            << filename_buf
            << " as shared information ..." 
            << std::endl;

  std::vector<bool> proc_flag(n_rank, false);
  std::vector<BinaryBuffer<> > proc_buf_in(n_rank), proc_buf_out(n_rank);
  std::vector<AFEPack::ostream<> > proc_os(n_rank);
  for (u_int i = 0;i < n_rank;++ i) {
    proc_os[i].set_buffer(proc_buf_out[i]);
  }

  int n_item;
  std::vector<std::pair<std::vector<std::pair<int,int> >, int> > 
    shared_geo_0d, shared_geo_1d, shared_geo_2d, shared_geo_3d;        
  if (DIM >= 0) {
    std::vector<int> n_rank_shared_0d(n_rank, 0);
    int n_shared_geo_0d;
    is >> n_shared_geo_0d >> n_shared_geo_0d;
    std::cerr << "\t# 0d-shared geometry: " << n_shared_geo_0d << std::endl;
    shared_geo_0d.resize(n_shared_geo_0d);
    for (int i = 0;i < n_shared_geo_0d;++ i) {
      is >> n_item; 
      shared_geo_0d[i].first.resize(n_item);
      for (int j = 0;j < n_item;++ j) {
        is >> shared_geo_0d[i].first[j].first   
           >> shared_geo_0d[i].first[j].second; 
        n_rank_shared_0d[shared_geo_0d[i].first[j].first] += 1;
      }
      is >> shared_geo_0d[i].second; 
    }

    for (int i = 0;i < n_rank;++ i) {
      proc_os[i] << n_rank_shared_0d[i];
      if (n_rank_shared_0d[i] > 0) {
        proc_flag[i] = true;
      }
    }
    for (int i = 0;i < n_shared_geo_0d;++ i) {
      n_item = shared_geo_0d[i].first.size();
      int local_idx;
      for (int j = 0;j < n_item;++ j) {
        int rnk = shared_geo_0d[i].first[j].first;
        int idx = shared_geo_0d[i].first[j].second;
        if (rnk == rank) {
          local_idx = idx;
          break;
        }
      }
      for (int j = 0;j < n_item;++ j) {
        int rnk = shared_geo_0d[i].first[j].first;
        int idx = shared_geo_0d[i].first[j].second;
        if (rnk != rank) {
          proc_os[rnk] << idx << geo_0d[local_idx];
        }
      }
    }
  }

  if (DIM >= 1) {
    std::vector<int> n_rank_shared_1d(n_rank, 0);
    int n_shared_geo_1d;
    is >> n_shared_geo_1d >> n_shared_geo_1d;
    std::cerr << "\t# 1d-shared geometry: " << n_shared_geo_1d << std::endl;
    shared_geo_1d.resize(n_shared_geo_1d);        
    for (int i = 0;i < n_shared_geo_1d;++ i) {
      is >> n_item; 
      shared_geo_1d[i].first.resize(n_item);
      for (int j = 0;j < n_item;++ j) {
        is >> shared_geo_1d[i].first[j].first   
           >> shared_geo_1d[i].first[j].second; 
        n_rank_shared_1d[shared_geo_1d[i].first[j].first] += 1;
      }
      is >> shared_geo_1d[i].second; 
    }

    for (int i = 0;i < n_rank;++ i) {
      proc_os[i] << n_rank_shared_1d[i];
      if (n_rank_shared_1d[i] > 0) {
        proc_flag[i] = true;
      }
    }
    for (int i = 0;i < n_shared_geo_1d;++ i) {
      n_item = shared_geo_1d[i].first.size();
      int local_idx;
      for (int j = 0;j < n_item;++ j) {
        int rnk = shared_geo_1d[i].first[j].first;
        int idx = shared_geo_1d[i].first[j].second;
        if (rnk == rank) {
          local_idx = idx;
          break;
        }
      }
      for (int j = 0;j < n_item;++ j) {
        int rnk = shared_geo_1d[i].first[j].first;
        int idx = shared_geo_1d[i].first[j].second;
        if (rnk != rank) {
          proc_os[rnk] << idx << geo_1d[local_idx];
        }
      }
    }
  }

  if (DIM >= 2) {
    std::vector<int> n_rank_shared_2d(n_rank, 0);
    int n_shared_geo_2d;
    is >> n_shared_geo_2d >> n_shared_geo_2d;
    std::cerr << "\t# 2d-shared geometry: " << n_shared_geo_2d << std::endl;
    shared_geo_2d.resize(n_shared_geo_2d);        
    for (int i = 0;i < n_shared_geo_2d;++ i) {
      is >> n_item; 
      shared_geo_2d[i].first.resize(n_item);
      for (int j = 0;j < n_item;++ j) {
        is >> shared_geo_2d[i].first[j].first   
           >> shared_geo_2d[i].first[j].second; 
        n_rank_shared_2d[shared_geo_2d[i].first[j].first] += 1;
      }
      is >> shared_geo_2d[i].second; 
    }

    for (int i = 0;i < n_rank;++ i) {
      proc_os[i] << n_rank_shared_2d[i];
      if (n_rank_shared_2d[i] > 0) {
        proc_flag[i] = true;
      }
    }
    for (int i = 0;i < n_shared_geo_2d;++ i) {
      n_item = shared_geo_2d[i].first.size();
      int local_idx;
      for (int j = 0;j < n_item;++ j) {
        int rnk = shared_geo_2d[i].first[j].first;
        int idx = shared_geo_2d[i].first[j].second;
        if (rnk == rank) {
          local_idx = idx;
          break;
        }
      }
      for (int j = 0;j < n_item;++ j) {
        int rnk = shared_geo_2d[i].first[j].first;
        int idx = shared_geo_2d[i].first[j].second;
        if (rnk != rank) {
          proc_os[rnk] << idx << geo_2d[local_idx];
        }
      }
    }
  }

  if (DIM >= 3) {
    std::vector<int> n_rank_shared_3d(n_rank, 0);
    int n_shared_geo_3d;
    is >> n_shared_geo_3d >> n_shared_geo_3d;
    std::cerr << "\t# 3d-shared geometry: " << n_shared_geo_3d << std::endl;
    shared_geo_3d.resize(n_shared_geo_3d);        
    for (int i = 0;i < n_shared_geo_3d;++ i) {
      is >> n_item; 
      shared_geo_3d[i].first.resize(n_item);
      for (int j = 0;j < n_item;++ j) {
        is >> shared_geo_3d[i].first[j].first   
           >> shared_geo_3d[i].first[j].second; 
        n_rank_shared_3d[shared_geo_3d[i].first[j].first] += 1;
      }
      is >> shared_geo_3d[i].second; 
    }

    for (int i = 0;i < n_rank;++ i) {
      proc_os[i] << n_rank_shared_3d[i];
      if (n_rank_shared_3d[i] > 0) {
        proc_flag[i] = true;
      }
    }
    for (int i = 0;i < n_shared_geo_3d;++ i) {
      n_item = shared_geo_3d[i].first.size();
      int local_idx;
      for (int j = 0;j < n_item;++ j) {
        int rnk = shared_geo_3d[i].first[j].first;
        int idx = shared_geo_3d[i].first[j].second;
        if (rnk == rank) {
          local_idx = idx;
          break;
        }
      }
      for (int j = 0;j < n_item;++ j) {
        int rnk = shared_geo_3d[i].first[j].first;
        int idx = shared_geo_3d[i].first[j].second;
        if (rnk != rank) {
          proc_os[rnk] << idx << geo_3d[local_idx];
        }
      }
    }
  }

  std::cerr << "\ttransfering data ..." << std::endl;
  int tag = get_random_tag(_comm), n_req = 0;
  proc_flag[rank] = false; 
  MPI_Request request[2*n_rank];
  MPI_Status status[2*n_rank];
  for (int i = 0;i < n_rank;++ i) {
    if (proc_flag[i] == false) continue;
    MPI_Isend(proc_buf_out[i].start_address(), proc_buf_out[i].size(), MPI_CHAR,
              i, tag, _comm, &request[n_req ++]);
  }
  for (int i = 0;i < n_rank;++ i) {
    MPI_Status mpi_status;
    if (proc_flag[i] == false) continue;
    MPI_Probe(i, tag, _comm, &mpi_status);
    int msg_size;
    MPI_Get_count(&mpi_status, MPI_CHAR, &msg_size);
    proc_buf_in[i].resize(msg_size);
    MPI_Irecv(proc_buf_in[i].start_address(), msg_size, MPI_CHAR,
              i, tag, _comm, &request[n_req ++]);
  }
  MPI_Waitall(n_req, request, status);

  std::cerr << "\tdecoding data received ..." << std::endl;
  for (int i = 0;i < n_rank;++ i) { 
    if (proc_flag[i] == false) continue;
    AFEPack::istream<> proc_is(proc_buf_in[i]);
    int n_shared_item, local_idx;

    if (DIM >= 0) { 
#define GDIM 0
#define GEO HGeometry<GDIM,this_t::dow>
#define OBJ Shared_object<GEO>
      proc_is >> n_shared_item;
      for (int j = 0;j < n_shared_item;++ j) {
        GEO * remote_obj, * local_obj;
        proc_is >> local_idx >> remote_obj;
        local_obj = geo_0d[local_idx];
        OBJ * p_info = get_shared_info(*local_obj);
        if (p_info == NULL) p_info = new_shared_info(*local_obj);
        p_info->add_clone(i, remote_obj);
      }
#undef OBJ
#undef GEO
#undef GDIM
    }

    if (DIM >= 1) {
#define GDIM 1
#define GEO HGeometry<GDIM,this_t::dow>
#define OBJ Shared_object<GEO>
      proc_is >> n_shared_item;
      for (int j = 0;j < n_shared_item;++ j) {
        GEO * remote_obj, * local_obj;
        proc_is >> local_idx >> remote_obj;
        local_obj = geo_1d[local_idx];
        OBJ * p_info = get_shared_info(*local_obj);
        if (p_info == NULL) p_info = new_shared_info(*local_obj);
        p_info->add_clone(i, remote_obj);
      }
#undef OBJ
#undef GEO
#undef GDIM
    }

    if (DIM >= 2) {
#define GDIM 2
#define GEO HGeometry<GDIM,this_t::dow>
#define OBJ Shared_object<GEO>
      proc_is >> n_shared_item;
      for (int j = 0;j < n_shared_item;++ j) {
        GEO * remote_obj, * local_obj;
        proc_is >> local_idx >> remote_obj;
        local_obj = geo_2d[local_idx];
        OBJ * p_info = get_shared_info(*local_obj);
        if (p_info == NULL) p_info = new_shared_info(*local_obj);
        p_info->add_clone(i, remote_obj);
      }
#undef OBJ
#undef GEO
#undef GDIM
    }

    if (DIM >= 3) {
#define GDIM 3
#define GEO HGeometry<GDIM,this_t::dow>
#define OBJ Shared_object<GEO>
      proc_is >> n_shared_item;
      for (int j = 0;j < n_shared_item;++ j) {
        GEO * remote_obj, * local_obj;
        proc_is >> local_idx >> remote_obj;
        local_obj = geo_3d[local_idx];
        OBJ * p_info = get_shared_info(*local_obj);
        if (p_info == NULL) p_info = new_shared_info(*local_obj);
        p_info->add_clone(i, remote_obj);
      }
#undef OBJ
#undef GEO
#undef GDIM
    }
  }
}

TEMPLATE
void THIS::renumerateRootElement(void (*f)(double, double, double,
                                           double&,double&,double&))
{
  std::cerr << "Renumerating element of the mesh using Hibert space curve filling ..." 
            << std::flush;
  int n_ele = this->n_rootElement();
  std::vector<HGeometry<DIM,this_t::dow>*> tmp_ele(n_ele);
  std::vector<double> x(n_ele, 0.0), y(n_ele, 0.0), z(n_ele, 0.0);
  typename std::list<HGeometry<DIM,this_t::dow>*>::iterator
    the_ele = this->rootElement().begin();
  for (int i = 0;i < n_ele;++ i, ++ the_ele) {
    tmp_ele[i] = *the_ele;
    HGeometry<DIM,this_t::dow>& ele = **the_ele;
    const int& n_vtx = ele.n_vertex;
    for (int j = 0;j < n_vtx;++ j) {
      Point<this_t::dow>& pnt = *ele.vertex[j];
      x[i] += pnt[0]; y[i] += pnt[1];
      if (this_t::dow == 3) z[i] += pnt[2];
    }
    x[i] /= n_vtx; y[i] /= n_vtx;
    if (this_t::dow == 3) z[i] /= n_vtx;
  }
  std::vector<int> new_index(n_ele);
  hsfc_renumerate(n_ele, &x[0], &y[0], &z[0], &new_index[0], f);

  the_ele = this->rootElement().begin();
  for (int i = 0;i < n_ele;++ i, ++ the_ele) {
    *the_ele = tmp_ele[new_index[i]];
  }
  std::cerr << " OK!" << std::endl;
}


TEMPLATE
void THIS::eraseRootElement(u_int level) {
  for (u_int i = 0;i < level;++ i) {
    this->eraseRootElementOneLevel();
  }
}

TEMPLATE
void THIS::eraseRootElementOneLevel() {
  std::list<HGeometry<DIM,this_t::dow>*> new_root;
  typename std::list<HGeometry<DIM,this_t::dow>*>::iterator
    the_ele = this->rootElement().begin(),
    end_ele = this->rootElement().end();
  for (;the_ele != end_ele;++ the_ele) {
    for (u_int i = 0;i < HGeometry<DIM,this_t::dow>::n_child;++ i) {
      HGeometry<DIM,this_t::dow> * p_geo = (*the_ele)->child[i];
      if (p_geo == NULL) {
        std::cerr << "Macro element is not refined. The root elements are not erased!" 
                  << std::endl;
        abort();
      }
      new_root.push_back(p_geo);
      this->nullParent(*p_geo);
    }
  }
  new_root.swap(this->rootElement());

  property_id_t<bool> pid;
  new_property_id(pid);
  end_ele = new_root.end();
  the_ele = new_root.begin();
  for (;the_ele != end_ele;++ the_ele) {
    this->tryDeleteGeometry(*the_ele, pid);
  }
  free_property_id(pid);
  the_ele = new_root.begin();
  for (;the_ele != end_ele;++ the_ele) {
    this->deleteGeometry(*the_ele);
  }
}

#undef THIS
#undef TEMPLATE


#define TEMPLATE template <class FOREST>
#define THIS HGeometryMatcher<FOREST>

TEMPLATE
bool THIS::match_geometry() {
  typedef THIS this_t;

  std::cerr << "Rank " << _forest->rank()
            << ": Matching h-geometry among partitions ..." << std::endl;
  bool is_operated = false;
  do {
    _is_operated = false;

    if (this_t::dim >= 3)
      sync_data(_forest->_comm, _forest->_shared_list_3d, *this, 
                &this_t::template pack_match_geometry<3>,
                &this_t::template unpack_match_geometry<3>,
                &this_t::template is_pack_match_geometry<3>);
    if (this_t::dim >= 2)
      sync_data(_forest->_comm, _forest->_shared_list_2d, *this,
                &this_t::template pack_match_geometry<2>,
                &this_t::template unpack_match_geometry<2>,
                &this_t::template is_pack_match_geometry<2>);
    if (this_t::dim >= 1)
      sync_data(_forest->_comm, _forest->_shared_list_1d, *this,
                &this_t::template pack_match_geometry<1>,
                &this_t::template unpack_match_geometry<1>,
                &this_t::template is_pack_match_geometry<1>);

    int src = _is_operated?1:0, dst = 0; 
    MPI_Allreduce(&src, &dst, 1, MPI_INT, MPI_SUM, _forest->_comm);
    if (_is_operated = (dst > 0)) {
      is_operated = true;
    }
  } while (_is_operated);
  return is_operated;
}

TEMPLATE
template <int GDIM> bool 
THIS::is_pack_match_geometry(HGeometry<GDIM,dow> * geo) const {
  u_int n_unsent_vtx = 0;
  u_int n_unsent_bnd = 0;
  if (GDIM == 1) { 

    u_int n_vertex = geo->n_vertex;
    for (u_int i = 0;i < n_vertex;++ i) {
      if (! is_shared_info_sent(*geo->vertex[i])) {
        n_unsent_vtx += 1;
      }
    }
  } else { 

    u_int n_boundary = geo->n_boundary;
    for (u_int i = 0;i < n_boundary;++ i) {
      HGeometry<GDIM-1,dow> * bnd = geo->boundary[i];
      if (bnd->parent == NULL && 
          ! is_shared_info_sent(*bnd)) {
        n_unsent_bnd += 1;
      }
    }
  }

  u_int n_child = 0;
  if (geo->isRefined() && 
      (! _forest->is_shared_info_sent(*geo->child[0]))) {
    n_child = geo->n_child;
  }
  return ((n_unsent_vtx > 0) ||
          (n_unsent_bnd > 0) ||
          (n_child > 0));
}

TEMPLATE
template <int GDIM> void 
THIS::pack_match_geometry(HGeometry<GDIM,dow> * geo,
                          int remote_rank,
                          AFEPack::ostream<>& os) {
  Shared_object<HGeometry<GDIM,dow> > * so = get_shared_info(*geo);
  Point<dow> grp; 

  bool is_last_entry = (so->rbegin()->first == remote_rank);

  if (GDIM == 1) { 

    u_int n_vertex = geo->n_vertex;
    u_int n_unsent_vtx = 0;
    for (u_int i = 0;i < n_vertex;++ i) {
      if (! is_shared_info_sent(*geo->vertex[i])) {
        n_unsent_vtx += 1;
      }
    }
    if (n_unsent_vtx > 0) { 
      this->_is_operated = true;
    }

    os << n_unsent_vtx; 
    for (u_int i = 0;i < n_vertex;++ i) {
      HGeometry<0,dow> * vtx = geo->vertex[i];
      if (! is_shared_info_sent(*vtx)) {
        os << vtx; 
        for (u_int n = 0;n < this_t::dow;++ n) {
          os << (*vtx)[n]; 
        }
        if (is_last_entry) { 
          set_shared_info_sent(*vtx); 
        }
      }
    }
  } else { 

    u_int n_boundary = geo->n_boundary;
    u_int n_unsent_bnd = 0;
    for (u_int i = 0;i < n_boundary;++ i) {
      HGeometry<GDIM-1,dow> * bnd = geo->boundary[i];
      if (bnd->parent == NULL && 
          ! is_shared_info_sent(*bnd)) {
        n_unsent_bnd += 1;
      }
    }

    if (n_unsent_bnd > 0) { 
      this->_is_operated = true;
    }

    os << n_unsent_bnd; 
    for (u_int i = 0;i < n_boundary;++ i) {
      HGeometry<GDIM-1,dow> * bnd = geo->boundary[i];
      if (bnd->parent == NULL && 
          ! is_shared_info_sent(*bnd)) {
        geometry_reference_point(*bnd, grp);
        os << bnd << grp;
        if (is_last_entry) { 
          set_shared_info_sent(*bnd); 
        }
      }
    }
  }

  if (geo->isRefined() && 
      (! _forest->is_shared_info_sent(*geo->child[0]))) {
    u_int n_child = geo->n_child;
    os << n_child; 
    for (u_int i = 0;i < n_child;++ i) {
      HGeometry<GDIM,dow> * chd = geo->child[i]; 
      Shared_object<HGeometry<GDIM,dow> > * p_info = get_shared_info(*chd);
      if (p_info == NULL) p_info = new_shared_info(*chd); 
      geometry_reference_point(*chd, grp); 
      os << chd << grp; 
      if (is_last_entry) { 
        set_shared_info_sent(*chd); 
      } 
    }

    this->_is_operated = true; 
  } else { 
    u_int n_child = 0;
    os << n_child;
  }
}

TEMPLATE
template <int GDIM> void 
THIS::unpack_match_geometry(HGeometry<GDIM,dow> * geo,
                            int remote_rank,
                            AFEPack::istream<>& is) {
  double h = distance(*geo->vertex[0], *geo->vertex[1]);

  if (GDIM == 1) { 

    u_int n_unsent_vtx = 0;
    is >> n_unsent_vtx;
    std::vector<HGeometry<0,dow>*> remote_vtx(n_unsent_vtx);
    std::vector<Point<dow> > pnt(n_unsent_vtx);
    for (u_int i = 0;i < n_unsent_vtx;++ i) {
      is >> remote_vtx[i] >> pnt[i];
    }

    u_int n_vertex = geo->n_vertex;
    for (u_int i = 0;i < n_unsent_vtx;++ i) {
      for (u_int j = 0;j < n_vertex;++ j) {
        HGeometry<0,dow> * vtx = geo->vertex[j];
        int type = matcher().value(pnt[i], *vtx, 1.0e-4*h);
        if (type >= 0) {
          Shared_object<HGeometry<0,dow> > * 
            p_info = get_shared_info(*vtx);
          if (p_info == NULL) p_info = new_shared_info(*vtx);
          p_info->add_clone(remote_rank, type, remote_vtx[i]);
          break;
        }
      }
    }
  } else { 

    u_int n_unsent_bnd = 0;
    is >> n_unsent_bnd;
    std::vector<HGeometry<GDIM-1,dow>*> remote_bnd(n_unsent_bnd);
    std::vector<Point<dow> > remote_grp(n_unsent_bnd);
    for (u_int i = 0;i < n_unsent_bnd;++ i) {
      is >> remote_bnd[i] >> remote_grp[i];
    }

    Point<dow> grp;
    u_int n_boundary = geo->n_boundary;
    for (u_int i = 0;i < n_unsent_bnd;++ i) {
      for (u_int j = 0;j < n_boundary;++ j) {
        HGeometry<GDIM-1,dow> * bnd = geo->boundary[j];
        if (bnd->parent != NULL) continue;
        geometry_reference_point(*bnd, grp);
        int type = matcher().value(remote_grp[i], grp, 1.0e-4*h);
        if (type >= 0) {
          Shared_object<HGeometry<GDIM-1,dow> > *
            p_info = get_shared_info(*bnd);
          if (p_info == NULL) p_info = new_shared_info(*bnd);
          p_info->add_clone(remote_rank, type, remote_bnd[i]);
          break;
        }
      }
    }
  }
      
  u_int n_child;
  is >> n_child;
  if (n_child != 0) { 
    assert (n_child == geo->n_child); 
        
    std::vector<HGeometry<GDIM,dow>*> remote_chd(n_child);
    std::vector<Point<dow> > remote_grp(n_child);
    for (u_int i = 0;i < n_child;++ i) {
      is >> remote_chd[i] >> remote_grp[i];
    }

    if (_forest->is_dummy(*geo)) return;

    if (! geo->isRefined()) geo->refine();

    Point<dow> grp;
    u_int n_matched_child = 0;
    for (u_int i = 0;i < n_child;++ i) {
      HGeometry<GDIM,dow> * chd = geo->child[i];
      geometry_reference_point(*chd, grp);
      for (u_int j = 0;j < n_child;++ j) {
        if (remote_chd[j] == NULL) continue; 

        int type = matcher().value(grp, remote_grp[j], 1.0e-4*h);
        if (type >= 0) {
          Shared_object<HGeometry<GDIM,dow> > * p_info = get_shared_info(*chd); 
          if (p_info == NULL) p_info = new_shared_info(*chd);
          p_info->add_clone(remote_rank, type, remote_chd[j]);
          remote_chd[j] = NULL;
          n_matched_child += 1;
          break;
        }
      }
    }
    assert (n_matched_child == n_child);
  }
}

TEMPLATE
bool THIS::sync_is_used() {
  typedef THIS this_t;
  std::cerr << "Rank " << _forest->rank()
            << ": Sync is_used flag ..." << std::endl;
  _is_operated = false;
  if (this_t::dim >= 1)
    sync_data(_forest->_comm, _forest->_shared_list_1d, *this, 
              &this_t::template pack_sync_is_used<1>,
              &this_t::template unpack_sync_is_used<1>);
  if (this_t::dim >= 2)
    sync_data(_forest->_comm, _forest->_shared_list_2d, *this, 
              &this_t::template pack_sync_is_used<2>,
              &this_t::template unpack_sync_is_used<2>);
  if (this_t::dim >= 3)
    sync_data(_forest->_comm, _forest->_shared_list_3d, *this, 
              &this_t::template pack_sync_is_used<3>,
              &this_t::template unpack_sync_is_used<3>);

  int src = _is_operated?1:0, dst = 0; 
  MPI_Allreduce(&src, &dst, 1, MPI_INT, MPI_SUM, _forest->_comm);
  return (dst > 0);
}

TEMPLATE
template <int GDIM> void 
THIS::pack_sync_is_used(HGeometry<GDIM,dow> * geo,
                        int remote_rank,
                        AFEPack::ostream<>& os) {
  bool is_used = _tools.isGeometryUsed(*geo);
  os << is_used;
}

TEMPLATE
template <int GDIM> void 
THIS::unpack_sync_is_used(HGeometry<GDIM,dow> * geo,
                          int remote_rank,
                          AFEPack::istream<>& is) {
  bool is_used;
  is >> is_used;

  if (is_used) {
    if (! _tools.isGeometryUsed(*geo)) {
      _tools.setGeometryUsed(*geo);
      _is_operated = true;
    }
  }
}

#undef THIS


#define THIS BirdView<FOREST>

TEMPLATE void 
THIS::semiregularize() {
  forest_t& forest = this->getForest();

  if (! forest.lock()) { // 对几何遗传树进行加锁
    std::cerr << "The hierarchy geometry tree is locked, aborting ...";
    abort();
  }

  if (forest.rank() == 0)
    std::cerr << "Semiregularizing the mesh ...  " << std::flush;
  const char * timer = "-/|\\";
  int n_element_refined = 0;
  this->prepareSemiregularize();

  int round = 0;
  bool is_operated;
  HGeometryMatcher<forest_t> matcher(forest);
  do {
    bool flag;
    do {
      flag = false;
      this->semiregularizeHelper(flag, n_element_refined);
    } while (flag);

    if (forest.rank() == 0)
      std::cerr << "\b" << timer[round] << std::flush;
    round = (round + 1)%4;

    is_operated = false;
    if (1) { // forest.n_rank() > 1) {
      is_operated |= matcher.match_geometry();
      is_operated |= matcher.sync_is_used();
    }
  } while (is_operated);

  if (forest.rank() == 0) 
    std::cerr << "\bOK!\n"
              << "\t" << n_element_refined 
              << " elements refined in semiregularization."
              << std::endl;
}

TEMPLATE
void THIS::eraseRootElement(u_int level) {
  for (u_int i = 0;i < level;++ i) {
    this->eraseRootElementOneLevel();
  }
}

TEMPLATE
void THIS::eraseRootElementOneLevel() {
  typedef typename base_t::container_t container_t;
  typedef typename base_t::element_t element_t;
  container_t new_root;
  typename container_t::iterator
    the_ele = this->rootElement().begin(),
    end_ele = this->rootElement().end();
  for (;the_ele != end_ele;++ the_ele) {
    for (u_int i = 0;i < element_t::n_child;++ i) {
      element_t * p_geo = (*the_ele)->child[i];
      if (p_geo == NULL) {
        std::cerr << "Macro element is not refined. The root elements CAN'T be erased!" 
                  << std::endl;
        abort();
      }
      new_root.push_back(p_geo);
      p_geo->parent = NULL;
    }
  }
  new_root.swap(this->rootElement());

  the_ele = new_root.begin();
  end_ele = new_root.end();
  for (;the_ele != end_ele;++ the_ele) {
    delete *the_ele;
  }
}

#undef THIS
#undef TEMPLATE