AFEPack
MPI_PeriodHandler.h
浏览该文件的文档。
00001 
00057 #ifndef __MPI_PeriodHandler_h__
00058 #define __MPI_PeriodHandler_h__
00059 
00060 namespace MPI {
00061   namespace Periodic {
00062 
00067     template <int DOW>
00068       struct PointDistance {
00069         int bmark;
00070         std::vector<double> period;
00074         int value(const Point<DOW>& p0,
00075                   const Point<DOW>& p1,
00076                   double tol) const {
00077           int type = 0;
00078           for (int i = 0;i < DOW && (type >= 0);++ i) {
00079             double a = fabs(p0[i] - p1[i]);
00080             if (a > tol) {
00081               a = fabs(a - period[i]);
00082               type = (a < tol)?1:-1;
00083             }
00084           }
00085           return type;
00086         }
00100          void readConfigFile(const std::string& file) {
00101           filtering_istream is;
00102           OpenFilteredStream(file, is);
00103           is >> bmark;
00104           period.resize(DOW);
00105           for (int i = 0;i < DOW;++ i) {
00106             is >> period[i];
00107           }
00108         }
00109       };
00110 
00111     namespace details {
00112 
00116       template <int DIM, int DOW>
00117         void collectGeometry(HGeometry<DIM,DOW>& geo,
00118                              int bnd_mark,
00119                              property_id_t<bool>& pid,
00120                              std::vector<std::vector<void *> >& bnd_geo) {
00121         if (geo.get_property(pid) == NULL) {
00122           if (geo.bmark == bnd_mark) {
00123             geo.new_property(pid);
00124             bnd_geo[DIM].push_back((void *)(&geo));
00125             if (geo.isRefined()) { 
00126               for (u_int i = 0;i < geo.n_child;++ i) {
00127                 collectGeometry(*(geo.child[i]), bnd_mark, pid, bnd_geo);
00128               }
00129             }
00130           }
00131           if (DIM > 1) { 
00132             for (u_int i = 0;i < geo.n_boundary;++ i) {
00133               collectGeometry(*(geo.boundary[i]), bnd_mark, pid, bnd_geo);
00134             }
00135           } else if (DIM == 1) { 
00136             for (u_int i = 0;i < geo.n_vertex;++ i) {
00137               collectGeometry(*(geo.vertex[i]), bnd_mark, pid, bnd_geo);
00138             }
00139           }
00140         }
00141       }
00142 
00146       template <class GEO>
00147         void geometry_reference_point(const GEO& geo, 
00148                                       Point<GEO::dow>& pnt) {
00149         int n_vtx = geo.n_vertex;
00150         for (int n = 0;n < GEO::dow;++ n) {
00151           pnt[n] = 0.0;
00152           for (int i = 0;i < n_vtx;++ i) {
00153             pnt[n] += (*geo.vertex[i])[n];
00154           }
00155           pnt[n] /= n_vtx;
00156         }
00157       }
00158       template <int DOW>
00159         void geometry_reference_point(const HGeometry<0,DOW>& geo, 
00160                                       Point<DOW>& pnt) {
00161         pnt = geo;
00162       }
00163 
00164       struct matchTolerence {
00165         double operator()() const {
00166           return 1.0e-08; 
00167         }
00168       };
00169 
00184       template <int DIM, class FOREST>
00185         void matchRankGeometryBuffer(FOREST& forest,
00186                                      int local_rank,
00187                                      std::vector<void *>& local_geo,
00188                                      std::vector<Point<FOREST::dow> >& local_grp,
00189                                      int remote_rank,
00190                                      std::vector<void *>& remote_geo,
00191                                      std::vector<Point<FOREST::dow> >& remote_grp) {
00192         typedef HGeometry<DIM,FOREST::dow> geometry_t;
00193         matchTolerence h;
00194 
00195         u_int n_remote_geo = remote_geo.size();
00196         u_int n_local_geo = local_geo.size();
00197         for (u_int i = 0;i < n_local_geo;++ i) {
00198           geometry_t * p_geo_i = (geometry_t *)(local_geo[i]);
00199           Point<FOREST::dow>& pi = local_grp[i];
00200           for (u_int j = 0;j < n_remote_geo;++ j) {
00201             geometry_t * p_geo_j = (geometry_t *)(remote_geo[j]);
00202             Point<FOREST::dow>& pj = remote_grp[j];
00203 
00204             int type = forest.matcher().value(pi, pj, h());
00205             if (type > 0) { 
00206               Shared_object<geometry_t> * 
00207                 p_info = forest.get_shared_info(*p_geo_i);
00208               if (p_info == NULL) p_info = forest.new_shared_info(*p_geo_i);
00209               p_info->add_clone(remote_rank, type, p_geo_j);
00210             }
00215           }
00216         }
00217       }
00218 
00223       template <int DIM, class FOREST>
00224         void matchRankGeometry(FOREST& forest,
00225                                int shift,
00226                                std::vector<void *>& geo,
00227                                std::vector<Point<FOREST::dow> >& grp,
00228                                BinaryBuffer<>& buf, 
00229                                int tag) {
00230         if (shift == 0) {
00231           matchRankGeometryBuffer<DIM,FOREST>(forest, shift, geo, grp, forest.rank(), geo, grp);
00232         } else {
00233           BinaryBuffer<> ibuf;
00234           int rank = forest.rank(), n_rank = forest.n_rank();
00235           int src = rank - shift;
00236           int dst = rank + shift;
00237           if (src < 0) src += n_rank;
00238           if (dst >= n_rank) dst -= n_rank;
00239           
00240           MPI_Comm comm = forest.communicator();
00241           int send_size = buf.size(), recv_size;
00242           MPI_Status status;
00243           MPI_Sendrecv(&send_size, 1, MPI_INT, dst, tag, 
00244                        &recv_size, 1, MPI_INT, src, tag, 
00245                        comm, &status);
00246 
00247           ibuf.resize(recv_size);
00248           MPI_Sendrecv(buf.start_address(), send_size, MPI_CHAR, dst, tag,
00249                        ibuf.start_address(), recv_size, MPI_CHAR, src, tag, 
00250                        comm, &status);
00251 
00252           AFEPack::istream<BinaryBuffer<> > is(ibuf);
00253           u_int n_remote_geo;
00254           is >> n_remote_geo;
00255           if (n_remote_geo > 0) {
00256             std::vector<void *> remote_geo(n_remote_geo);
00257             std::vector<Point<FOREST::dow> > remote_grp(n_remote_geo);
00258             for (u_int i = 0;i < n_remote_geo;++ i) {
00259               is >> remote_geo[i] >> remote_grp[i];
00260             }
00261             matchRankGeometryBuffer<DIM,FOREST>(forest, rank, geo, grp, 
00262                                                 src, remote_geo, remote_grp);
00263           }
00264         }
00265       }
00266 
00270       template <int DIM, class FOREST>
00271         void matchGeometry(FOREST& forest,
00272                            std::vector<void *>& geo) {
00273         typedef HGeometry<DIM,FOREST::dow> geometry_t;
00274 
00278         BinaryBuffer<> buf;
00279         AFEPack::ostream<BinaryBuffer<> > os(buf);
00280         u_int n = geo.size();
00281         std::vector<Point<FOREST::dow> > grp(n);
00282         os << n;
00283         for (u_int i = 0;i < n;++ i) {
00284           geometry_t * p_geo = (geometry_t *)(geo[i]);
00285           geometry_reference_point(*p_geo, grp[i]);
00286           os << geo[i] << grp[i];
00287         }
00288 
00289         int tag = 9996;
00290 
00291         int n_rank = forest.n_rank();
00292         for (int i = 0; i < n_rank;++ i) {
00293           matchRankGeometry<DIM,FOREST>(forest, i, geo, grp, buf, tag);
00294         }
00295       }
00296 
00297     }
00298 
00302     template <class FOREST>
00303       void readConfigFile(FOREST& forest,
00304                           const std::string& file) {
00305       forest.matcher().readConfigFile(file);
00306     }
00307 
00329     template <class FOREST>
00330       void periodizeForest(FOREST& forest) {
00331       int bnd_mark = forest.matcher().bmark;
00332 
00334       std::vector<std::vector<void *> > bnd_geo(FOREST::dim + 1);
00335       property_id_t<bool> pid;
00336       new_property_id(pid);
00337       typename FOREST::RootIterator 
00338         the_geo = forest.beginRootElement(),
00339         end_geo = forest.endRootElement();
00340       for (;the_geo != end_geo;++ the_geo) {
00341         details::collectGeometry(*the_geo, bnd_mark,
00342                                  pid, bnd_geo);
00343       }
00344       free_property_id(pid);
00345 
00347       for (int i = 0;i <= FOREST::dim;++ i) {
00348         switch (i) {
00349         case 0: details::matchGeometry<0,FOREST>(forest, bnd_geo[i]); break;
00350         case 1: details::matchGeometry<1,FOREST>(forest, bnd_geo[i]); break;
00351         case 2: details::matchGeometry<2,FOREST>(forest, bnd_geo[i]); break;
00352         case 3: details::matchGeometry<3,FOREST>(forest, bnd_geo[i]); break;
00353         }
00354       }
00355     }
00356 
00361     template <class FOREST>
00362       void periodizeForest(FOREST& forest,
00363                            const std::string& file) {
00365       readConfigFile(forest, file);
00366       periodizeForest(forest);
00367     }
00368 
00369   }
00370 }
00371 
00372 #endif // __MPI_PeriodHandler_h__
00373