AFEPack
MPI_Migration.h
浏览该文件的文档。
00001 
00011 #ifndef __MPI_Migration_h__
00012 #define __MPI_Migration_h__
00013 
00014 #include <mpi.h>
00015 #include <AFEPack/Miscellaneous.h>
00016 #include <AFEPack/Migration.details.h>
00017 #include <AFEPack/HGeometry.h>
00018 #include <AFEPack/mpi/MPI_HGeometry.h>
00019 #include <AFEPack/mpi/MPI_SyncProp.h>
00020 #include <AFEPack/FEMSpace.h>
00021 
00022 namespace Migration {
00023   namespace details {
00025 
00028     template <class GEO> void
00029       clear_geo_data_buffer(GEO& geo) {
00030       geo.buffer.clear();
00031       for (u_int i = 0;i < geo.n_boundary;++ i) {
00032         clear_geo_data_buffer(*geo.boundary[i]);
00033       }
00034       if (geo.isRefined()) {
00035         for (u_int i = 0;i < geo.n_child;++ i) {
00036           clear_geo_data_buffer<GEO>(*geo.child[i]);
00037         }
00038       }
00039     }
00040     template <int DOW> void
00041       clear_geo_data_buffer(HGeometry<0,DOW>& geo) {
00042       geo.buffer.clear();
00043     }
00044     template <int DOW> void
00045       clear_geo_data_buffer(HGeometry<1,DOW>& geo) {
00046       geo.buffer.clear();
00047       for (u_int i = 0;i < geo.n_vertex;++ i) {
00048         clear_geo_data_buffer(*geo.vertex[i]);
00049       }
00050       if (geo.isRefined()) {
00051         for (u_int i = 0;i < geo.n_child;++ i) {
00052           clear_geo_data_buffer(*geo.child[i]);
00053         }
00054       }
00055     }
00057 
00059 
00062     struct buffer_trait {
00063       buffer_t * p_buf;
00064 
00065       friend AFEPack::ostream<>&
00066       operator<<(AFEPack::ostream<>& os,
00067                  const buffer_trait& bt) {
00068         const buffer_t& buf = *(bt.p_buf);
00069         u_int n_data = buf.size();
00070         os << n_data;
00071         buffer_t::const_iterator
00072           the_data = buf.begin(),
00073           end_data = buf.end();
00074         for (;the_data != end_data;++ the_data) {
00075           os << the_data->first << the_data->second;
00076         }
00077         return os;
00078       }
00079       friend AFEPack::istream<>&
00080       operator>>(AFEPack::istream<>& is,
00081                  buffer_trait& bt) {
00082         buffer_t& buf = *(bt.p_buf);
00083         data_buffer_t spool;
00084         u_int n_data;
00085         is >> n_data;
00086         for (u_int i = 0;i < n_data;++ i) {
00087           data_id_t id;
00088           is >> id;
00089           if (buf.find(id) == buf.end()) {
00090             is >> buf[id]; 
00091           } else { 
00092             is >> spool;
00093           }
00094         }
00095         return is;
00096       }
00097     };
00098 
00099     template <int D, class HTREE> void
00100       sync_data_buffer_prepare(const HTREE& htree,
00101                                const property_id_t<buffer_trait>& pid) {
00102       enum {dow = HTREE::dow};
00103       typedef HGeometry<D,dow> GEO;
00104       const MPI::Shared_ptr_list<GEO>& so = htree.template get_shared_list<D>();
00105       typename MPI::Shared_ptr_list<GEO>::const_iterator
00106         the_entry = so.begin(), end_entry = so.end();
00107       for (;the_entry != end_entry;++ the_entry) {
00108         GEO * p_geo = the_entry->local_pointer();
00109         if (p_geo->get_property(pid) == NULL) {
00110           buffer_trait * p_bt = p_geo->new_property(pid);
00111           p_bt->p_buf = &(p_geo->buffer);
00112         }
00113       }
00114     }
00115 
00116     template <class HTREE> void
00117       sync_data_buffer_prepare(const HTREE& htree,
00118                                const property_id_t<buffer_trait>& pid,
00119                                int D) {
00120       switch (D) {
00121       case 0: sync_data_buffer_prepare<0,HTREE>(htree, pid); break;
00122       case 1: sync_data_buffer_prepare<1,HTREE>(htree, pid); break;
00123       case 2: sync_data_buffer_prepare<2,HTREE>(htree, pid); break;
00124       case 3: sync_data_buffer_prepare<3,HTREE>(htree, pid); break;
00125       }
00126     }
00127 
00128     template <class HTREE> void 
00129       sync_data_buffer(const HTREE& htree) {
00130       enum {dim = HTREE::dim};
00131       property_id_t<buffer_trait> pid;
00132       new_property_id(pid);
00133 
00134       MPI::PropSyncer<HTREE,buffer_trait> syncer(htree, pid);
00135       for (int i = 0;i <= dim;++ i) {
00136         sync_data_buffer_prepare(htree, pid, i);
00137         syncer.sync(i);
00138       }
00139       free_property_id(pid);
00140     }
00141 
00143 
00144 
00146 
00152     template <class HGEO, class T>
00153       void export_property(const HGEO& geo,
00154                            const data_id_t& data_id,
00155                            const property_id_t<T>& pid,
00156                            const property_id_t<>& flag) {
00157       if (geo.has_property(flag)) return; 
00158       if (! geo.has_property(pid)) return; 
00159 
00160       geo.new_property(flag); 
00161 
00162       T * p_prp = geo.get_property(pid);
00163       AFEPack::ostream<> os;
00164       os.set_buffer(get_buffer(geo, data_id, true));
00165       os << *p_prp;
00166 
00167       for (u_int i = 0;i < geo.n_vertex;++ i) {
00168         export_property(*geo.vertex[i], data_id, pid, flag);
00169       }
00170 
00171       for (u_int i = 0;i < geo.n_boundary;++ i) {
00172         export_property(*geo.boundary[i], data_id, pid, flag);
00173       }
00174 
00175       if (geo.isRefined()) {
00176         for (u_int i = 0;i < geo.n_child;++ i) {
00177           export_property(*geo.child[i], data_id, pid, flag);
00178         }
00179       }
00180     }
00181 
00182     template <class HGEO, class T>
00183       void import_property(const HGEO& geo,
00184                            const data_id_t& data_id,
00185                            const property_id_t<T>& pid,
00186                            const property_id_t<>& flag) {
00187       if (geo.has_property(flag)) return; 
00188 
00189       BinaryBuffer<>& buf = get_buffer(geo, data_id, false);
00190       if (&buf == NULL) return; 
00191 
00192       geo.new_property(flag); 
00193 
00194       T * p_prp = geo.get_property(pid);
00195       if (p_prp == NULL) {
00196         p_prp = geo.new_property(pid);
00197       }
00198 
00199       AFEPack::istream<> is;
00200       is.set_buffer(buf);
00201       is >> *p_prp;
00202 
00203       for (u_int i = 0;i < geo.n_vertex;++ i) {
00204         import_property(*geo.vertex[i], data_id, pid, flag);
00205       }
00206 
00207       for (u_int i = 0;i < geo.n_boundary;++ i) {
00208         import_property(*geo.boundary[i], data_id, pid, flag);
00209       }
00210 
00211       if (geo.isRefined()) {
00212         for (u_int i = 0;i < geo.n_child;++ i) {
00213           import_property(*geo.child[i], data_id, pid, flag);
00214         }
00215       }
00216     }
00217 
00218     template <class HTREE, class T>
00219       void export_property(const HTREE& tree,
00220                            const data_id_t& data_id,
00221                            const property_id_t<T>& pid) {
00222       property_id_t<> flag;
00223       new_property_id(flag);
00224 
00225       enum {dim = HTREE::dim, dow = HTREE::dow};
00226       typename HTREE::ConstRootIterator 
00227         the_ele = tree.beginRootElement(),
00228         end_ele = tree.endRootElement();
00229       for (;the_ele != end_ele;++ the_ele) {
00230         const HGeometry<dim,dow> * p_geo = &(*the_ele);
00231         do {
00232           if (p_geo->parent != NULL) {
00233             p_geo = p_geo->parent;
00234           } else break;
00235         } while (true);
00236         export_property(*p_geo, data_id, pid, flag);
00237       }
00238     }
00239 
00240     template <class HTREE, class T>
00241       void import_property(const HTREE& tree,
00242                            const data_id_t& data_id,
00243                            const property_id_t<T>& pid) {
00244       property_id_t<> flag;
00245       new_property_id(flag);
00246 
00247       enum {dim = HTREE::dim, dow = HTREE::dow};
00248       typename HTREE::ConstRootIterator 
00249         the_ele = tree.beginRootElement(),
00250         end_ele = tree.endRootElement();
00251       for (;the_ele != end_ele;++ the_ele) {
00252         const HGeometry<dim,dow> * p_geo = &(*the_ele);
00253         do {
00254           if (p_geo->parent != NULL) {
00255             p_geo = p_geo->parent;
00256           } else break;
00257         } while (true);
00258         import_property(*p_geo, data_id, pid, flag);
00259       }
00260     }
00261 
00263   }
00264 
00270   template <class HTREE> void
00271     clear_data_buffer(HTREE& tree) {
00272     enum {dim = HTREE::dim, dow = HTREE::dow};
00273     typename HTREE::RootIterator 
00274       the_ele = tree.beginRootElement(),
00275       end_ele = tree.endRootElement();
00276     for (;the_ele != end_ele;++ the_ele) {
00277       HGeometry<dim,dow> * p_geo = &(*the_ele);
00278       do {
00279         if (p_geo->parent != NULL) {
00280           p_geo = p_geo->parent;
00281         } else break;
00282       } while (true);
00283       details::clear_geo_data_buffer(*p_geo);
00284     }
00285   }
00286 
00287   void initialize(MPI_Comm comm); 
00288   data_id_t register_data_name(const data_name_t& dn); 
00289   void load_config(const std::string& filename); 
00290   void save_config(const std::string& filename); 
00291 
00292   void ensured_open_fstream(const std::string& filename,
00293                             std::ifstream& is);
00294   void ensured_open_filtering_stream(const std::string& filename,
00295                                      filtering_istream& is);
00296 
00304   template <class FUNC>
00305     void export_fe_func(const FUNC& fun,
00306                         const data_id_t& data_id) {
00307     typedef FUNC fe_func_t;
00308     typedef typename fe_func_t::fe_space_t fe_space_t;
00309     typedef typename fe_space_t::dof_info_t dof_info_t;
00310     typedef RegularMesh<fe_space_t::dim,fe_space_t::dow> mesh_t;
00311     const fe_space_t& sp = fun.femSpace();
00312     const mesh_t& mesh = dynamic_cast<const mesh_t&>(sp.mesh());
00313     AFEPack::ostream<> os;
00314     u_int n_dof = sp.n_dof();
00315     for (u_int i = 0;i < n_dof;++ i) {
00316       const dof_info_t& di = sp.dofInfo(i);
00317       const DOFIndex& dof_idx = sp.dofIndex(i);
00318       get_export_stream(mesh,
00319                         data_id,
00320                         dof_idx.dimension,
00321                         dof_idx.geometry_index,
00322                         os);
00323       os << di.identity
00324          << di.interp_point 
00325          << fun(i);
00326     }
00327   }
00337   template <class FUNC>
00338     void export_fe_func(const std::vector<FUNC>& fun,
00339                         const data_id_t& data_id) {
00340     typedef FUNC fe_func_t;
00341     typedef typename fe_func_t::fe_space_t fe_space_t;
00342     typedef typename fe_space_t::dof_info_t dof_info_t;
00343     typedef RegularMesh<fe_space_t::dim,fe_space_t::dow> mesh_t;
00344     const fe_space_t& sp = fun[0].femSpace();
00345     const mesh_t& mesh = dynamic_cast<const mesh_t&>(sp.mesh());
00346     AFEPack::ostream<> os;
00347     u_int n_dof = sp.n_dof();
00348     for (u_int i = 0;i < n_dof;++ i) {
00349       const dof_info_t& di = sp.dofInfo(i);
00350       const DOFIndex& dof_idx = sp.dofIndex(i);
00351       get_export_stream(mesh,
00352                         data_id,
00353                         dof_idx.dimension,
00354                         dof_idx.geometry_index,
00355                         os);
00356       os << di.identity
00357          << di.interp_point;
00358       for(u_int j=0; j<fun.size(); j++)
00359         os<< fun[j](i);
00360     }
00361   }
00362 
00370   template <class FUNC>
00371     void import_fe_func(FUNC& fun,
00372                         const data_id_t& data_id) {
00373     typedef FUNC fe_func_t;
00374     typedef typename fe_func_t::fe_space_t fe_space_t;
00375     typedef typename fe_space_t::dof_info_t dof_info_t;
00376     typedef RegularMesh<fe_space_t::dim,fe_space_t::dow> mesh_t;
00377     const fe_space_t& sp = fun.femSpace();
00378     const mesh_t& mesh = dynamic_cast<const mesh_t&>(sp.mesh());
00379     AFEPack::istream<> is;
00380     u_int n_dof = sp.n_dof();
00381     for (u_int i = 0;i < n_dof;++ i) {
00382       const dof_info_t& di = sp.dofInfo(i);
00383       const DOFIndex& dof_idx = sp.dofIndex(i);
00384       double D = di.interp_point.length();
00385       get_import_stream(mesh,
00386                         data_id,
00387                         dof_idx.dimension,
00388                         dof_idx.geometry_index,
00389                         is);
00390       do {
00391         dof_info_t dof_info;
00392         is >> dof_info.identity
00393            >> dof_info.interp_point 
00394            >> fun(i);
00395         if (!(dof_info.identity == di.identity)) continue;
00396         double D1 = dof_info.interp_point.length() + D;
00397         if (distance(dof_info.interp_point, 
00398                      di.interp_point) <= 1.0e-08*D1) break;
00399       } while (1);
00400     }
00401   }
00402 
00412   template <class FUNC>
00413     void import_fe_func(std::vector<FUNC> & fun,
00414                         const data_id_t& data_id) {
00415     typedef FUNC fe_func_t;
00416     typedef typename fe_func_t::fe_space_t fe_space_t;
00417     typedef typename fe_space_t::dof_info_t dof_info_t;
00418     typedef RegularMesh<fe_space_t::dim,fe_space_t::dow> mesh_t;
00419     const fe_space_t& sp = fun[0].femSpace();
00420     const mesh_t& mesh = dynamic_cast<const mesh_t&>(sp.mesh());
00421     AFEPack::istream<> is;
00422     u_int n_dof = sp.n_dof();
00423     for (u_int i = 0;i < n_dof;++ i) {
00424       const dof_info_t& di = sp.dofInfo(i);
00425       const DOFIndex& dof_idx = sp.dofIndex(i);
00426       double D = di.interp_point.length();
00427       get_import_stream(mesh,
00428                         data_id,
00429                         dof_idx.dimension,
00430                         dof_idx.geometry_index,
00431                         is);
00432       do {
00433         dof_info_t dof_info;
00434         is >> dof_info.identity
00435            >> dof_info.interp_point;
00436         for(u_int j=0; j<fun.size(); j++)
00437           is >> fun[j](i);
00438         if (!(dof_info.identity == di.identity)) continue;
00439         double D1 = dof_info.interp_point.length() + D;
00440         if (distance(dof_info.interp_point, 
00441                      di.interp_point) <= 1.0e-08*D1) break;
00442       } while (1);
00443     }
00444   }
00445 
00446   template <class HTREE, class T>
00447     void export_property(const HTREE& htree,
00448                          const data_id_t& data_id,
00449                          const property_id_t<T>& pid) {
00450     details::export_property(htree, data_id, pid);
00451   }
00452 
00453   template <class HTREE, class T>
00454     void import_property(const HTREE& htree,
00455                          const data_id_t& data_id,
00456                          const property_id_t<T>& pid) {
00457     details::import_property(htree, data_id, pid);
00458   }
00459 
00460   template <class HTREE>
00461     void sync_data_buffer(const HTREE& htree) {
00462     details::sync_data_buffer(htree);
00463   }
00464 
00465 }
00466 
00467 #endif // __MPI_Migration_h__
00468