AFEPack
|
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