AFEPack
MPI_DOF.templates.h
浏览该文件的文档。
00001 
00011 #define TEMPLATE template <class FOREST, class FESPACE>
00012 #define THIS GlobalIndex<FOREST,FESPACE>
00013 
00014 TEMPLATE int
00015 THIS::dof_primary_rank(u_int i) const {
00016   typedef RegularMesh<fe_space_t::dim,fe_space_t::dow> mesh_t;
00017   const mesh_t& mesh = dynamic_cast<const mesh_t&>(_fe_sp->mesh());
00018 
00019   const DOFIndex& dof_idx = _fe_sp->dofIndex(i);
00020   const int& dim = dof_idx.dimension;
00021   const int& gid = dof_idx.geometry_index;
00022   const HBuffer * p_buf = mesh.h_geometry(dim, gid);
00023 
00024   int result;
00025 #define MPI_DOF_PRIMARY_GEOMETRY(D)                                     \
00026   case D:                                                               \
00027     result = _forest->template primary_rank<D>                          \
00028     (*dynamic_cast<const HGeometry<D,fe_space_t::dow> *>(p_buf));       \
00029     break;
00030   switch (dim) {
00031     MPI_DOF_PRIMARY_GEOMETRY(0)
00032     MPI_DOF_PRIMARY_GEOMETRY(1)
00033     MPI_DOF_PRIMARY_GEOMETRY(2)
00034     MPI_DOF_PRIMARY_GEOMETRY(3)
00035   }
00036 #undef MPI_DOF_PRIMARY_GEOMETRY
00037   return result;
00038 }
00039 
00040 TEMPLATE bool
00041 THIS::_is_dof_on_primary_geometry(u_int i) const {
00042   typedef RegularMesh<fe_space_t::dim,fe_space_t::dow> mesh_t;
00043   const mesh_t& mesh = dynamic_cast<const mesh_t&>(_fe_sp->mesh());
00044 
00045   const DOFIndex& dof_idx = _fe_sp->dofIndex(i);
00046   const int& dim = dof_idx.dimension;
00047   const int& gid = dof_idx.geometry_index;
00048   const HBuffer * p_buf = mesh.h_geometry(dim, gid);
00049 
00050   bool result;
00051 #define MPI_DOF_IS_ON_PRIMARY_GEOMETRY(D)                               \
00052   case D:                                                               \
00053     result = _forest->template is_primary_geometry<D>                   \
00054     (*dynamic_cast<const HGeometry<D,fe_space_t::dow> *>(p_buf));       \
00055     break;
00056   switch (dim) {
00057     MPI_DOF_IS_ON_PRIMARY_GEOMETRY(0)
00058     MPI_DOF_IS_ON_PRIMARY_GEOMETRY(1)
00059     MPI_DOF_IS_ON_PRIMARY_GEOMETRY(2)
00060     MPI_DOF_IS_ON_PRIMARY_GEOMETRY(3)
00061   }
00062 #undef MPI_DOF_IS_ON_PRIMARY_GEOMETRY
00063   return result;
00064 }
00065 
00066 TEMPLATE bool
00067 THIS::is_dof_on_primary_geometry(u_int i) const {
00068   return _idopg[i];
00069 }
00070 
00072 TEMPLATE void
00073 THIS::sync_idx() {
00074   u_int n_dof = _fe_sp->n_dof();
00075   typedef typename fe_space_t::value_t value_t;
00076   typedef typename fe_space_t::dof_info_t dof_info_t;
00077   typedef RegularMesh<fe_space_t::dim,fe_space_t::dow> mesh_t;
00078   const mesh_t& mesh = dynamic_cast<const mesh_t&>(_fe_sp->mesh());
00079 
00081   new_property_id(_pid_global_idx);
00082   BinaryBuffer<> * p_buf;
00083   for (u_int i = 0;i < n_dof;++ i) {
00084     if ((*this)(i) == UNUSED_IDX) continue;
00085 
00086     const dof_info_t& di = _fe_sp->dofInfo(i);
00087     const DOFIndex& dof_idx = _fe_sp->dofIndex(i);
00088     const int& dim = dof_idx.dimension;
00089     const int& gid = dof_idx.geometry_index;
00090     if (_forest->is_geometry_shared(mesh, dim, gid)) {
00091       p_buf = mesh.get_property(dim, gid, _pid_global_idx);
00092       if (p_buf == NULL) {
00093         p_buf = mesh.new_property(dim, gid, _pid_global_idx);
00094       }
00095       AFEPack::ostream<> os(*p_buf);
00096       os << di.interp_point << di.identity << (*this)(i);
00097     }
00098   }
00099 
00101   MPI_Comm comm = _forest->communicator();
00102 #define MPI_DOF_SYNC_DATA_DIM(D)                                        \
00103   if (forest_t::dim >= D) {                                             \
00104     sync_data(comm, _forest->template get_shared_list<D>(), *this,      \
00105               &this_t::template pack_global_idx<D>,                     \
00106               &this_t::template unpack_global_idx<D>);                  \
00107   }
00108 
00109   MPI_DOF_SYNC_DATA_DIM(0);
00110   MPI_DOF_SYNC_DATA_DIM(1);
00111   MPI_DOF_SYNC_DATA_DIM(2);
00112   MPI_DOF_SYNC_DATA_DIM(3);
00113 
00114 #undef MPI_DOF_SYNC_DATA_DIM
00115 
00117   const typename forest_t::matcher_t& matcher = _forest->matcher();
00118 
00119   std::vector<bool> flag(n_dof, false);
00120   typename fe_space_t::ConstElementIterator
00121     the_ele = _fe_sp->beginElement(),
00122     end_ele = _fe_sp->endElement();
00123   for (;the_ele != end_ele;++ the_ele) {
00124     const std::vector<int>& ele_dof = the_ele->dof();
00125     u_int n_ele_dof = ele_dof.size();
00126     for (u_int j = 0;j < n_ele_dof;++ j) {
00127       const int& i = ele_dof[j];
00128       if (flag[i] || (*this)(i) != UNUSED_IDX) continue;
00129 
00130       const dof_info_t& di = _fe_sp->dofInfo(i);
00131       const DOFIndex& dof_idx = _fe_sp->dofIndex(i);
00132       const int& dim = dof_idx.dimension;
00133       const int& gid = dof_idx.geometry_index;
00134       if (_forest->is_geometry_shared(mesh, dim, gid)) {
00135         flag[i] = true;
00136 
00137         const GeometryBM& geo = the_ele->geometry();
00138         double h = distance(mesh.point(mesh.geometry(0, geo.vertex(0)).vertex(0)),
00139                             mesh.point(mesh.geometry(0, geo.vertex(1)).vertex(0)));
00140 
00141         p_buf = mesh.get_property(dim, gid, _pid_global_idx);
00142         assert (p_buf != NULL);
00143         AFEPack::istream<> is(*p_buf);
00144         do {
00145           dof_info_t dof_info;
00146           is >> dof_info.interp_point >> dof_info.identity >> (*this)(i);
00147           if (!(dof_info.identity == di.identity)) continue;
00148           if (matcher.value(dof_info.interp_point, di.interp_point, 1.0e-04*h) >= 0) break;
00149         } while (1);
00150       }
00151     }
00152   }
00153 
00154   free_property_id(_pid_global_idx);
00155 }
00156 
00157 TEMPLATE void
00158 THIS::build() {
00159   u_int n_dof = _fe_sp->n_dof();
00160   typedef typename fe_space_t::value_t value_t;
00161   typedef typename fe_space_t::dof_info_t dof_info_t;
00162   typedef RegularMesh<fe_space_t::dim,fe_space_t::dow> mesh_t;
00163   const mesh_t& mesh = dynamic_cast<const mesh_t&>(_fe_sp->mesh());
00164 
00166   base_t::resize(n_dof, UNUSED_IDX); 
00167   _idopg.resize(n_dof, false);
00168   int idx = 0;
00169   for (u_int i = 0;i < n_dof;++ i) {
00170     if (_is_dof_on_primary_geometry(i)) {
00171       idx += 1;
00172       (*this)(i) = 0;
00173       _idopg[i] = true;
00174     }
00175   }
00176   _n_primary_dof = idx;
00177 
00179   int global_idx;
00180   MPI_Comm comm = _forest->communicator();
00181   u_int n_rank = _forest->n_rank();
00182   MPI_Scan(&idx, &global_idx, 1, MPI_INT, MPI_SUM, comm);
00183   _n_global_dof = global_idx; 
00184   MPI_Bcast(&_n_global_dof, 1, MPI_UNSIGNED, n_rank - 1, comm);
00185   if (_forest->rank() == 0) {
00186     std::cerr << "Total #DOF: " << _n_global_dof << std::endl;
00187   }
00188 
00189   _first_idx = global_idx - idx; 
00190 
00191   std::vector<int>::iterator
00192     the_idx = base_t::begin(),
00193     end_idx = base_t::end();
00194   for (idx = _first_idx;the_idx != end_idx;++ the_idx) {
00195     if (*the_idx == UNUSED_IDX) continue;
00196     *the_idx = idx ++;
00197   }
00198 
00199   sync_idx(); 
00200 }
00201 
00202 TEMPLATE void
00203 THIS::build_primary_index(int * idx_ptr) const {
00204   std::vector<int>::const_iterator
00205     the_idx = base_t::begin(),
00206     end_idx = base_t::end();
00207   for (u_int i = 0;the_idx != end_idx;++ the_idx) {
00208     if (is_dof_on_primary_geometry(i ++)) {
00209       *idx_ptr = *the_idx;
00210       ++ idx_ptr;
00211     }
00212   }
00213 }
00214 
00215 TEMPLATE 
00216 template <class MAP> void
00217 THIS::build_epetra_map(MAP& map) const {
00218   int * idx_ptr = map.MyGlobalElements();
00219   std::vector<int>::const_iterator
00220     the_idx = base_t::begin(),
00221     end_idx = base_t::end();
00222   for (u_int i = 0;the_idx != end_idx;++ the_idx) {
00223     if (is_dof_on_primary_geometry(i ++)) {
00224       *idx_ptr = *the_idx;
00225       ++ idx_ptr;
00226     }
00227   }
00228 }
00229 
00230 TEMPLATE 
00231 template <class MAP, class INVMAP> void
00232 THIS::build_epetra_map(MAP& map, INVMAP& inv_map) const {
00233   int * idx_ptr = map.MyGlobalElements();
00234   std::vector<int>::const_iterator
00235     the_idx = base_t::begin(),
00236     end_idx = base_t::end();
00237   for (u_int i = 0, j = 0;the_idx != end_idx;++ the_idx) {
00238     if (is_dof_on_primary_geometry(i)) {
00239       *idx_ptr = *the_idx;
00240       inv_map[j ++] = i ++;
00241       ++ idx_ptr;
00242     }
00243   }
00244 }
00245 
00246 TEMPLATE
00247 template <class LC, class GC> void
00248   THIS::translate(const LC& lc, GC& gc) const {
00249   u_int n = lc.size();
00250   gc.resize(n);
00251   for (u_int i = 0;i < n;++ i) {
00252     gc[i] = (*this)(lc[i]);
00253   }
00254 }
00255 
00256 TEMPLATE
00257 template <int GDIM> void
00258 THIS::pack_global_idx(HGeometry<GDIM,fe_space_t::dow> * geo,
00259                       int remote_rank,
00260                       AFEPack::ostream<>& os) {
00261   BinaryBuffer<> * p_buf = geo->get_property(_pid_global_idx);
00262   u_int n = 0;
00263   if (p_buf != NULL) {
00264     n = p_buf->size();
00265     os << n;
00266     os.encode_binary(&(*p_buf)[0], n);
00267   } else {
00268     os << n;
00269   }
00270 }
00271 
00272 TEMPLATE
00273 template <int GDIM> void
00274 THIS::unpack_global_idx(HGeometry<GDIM,fe_space_t::dow> * geo,
00275                         int remote_rank,
00276                         AFEPack::istream<>& is) {
00277   u_int n = 0;
00278   is >> n;
00279   if (n > 0) {
00280     BinaryBuffer<> * p_buf;;
00281     assert ((p_buf  = geo->get_property(_pid_global_idx)) == NULL);
00282     p_buf = geo->new_property(_pid_global_idx);
00283     p_buf->resize(n);
00284     is.decode_binary(&(*p_buf)[0], n);
00285   }
00286 }
00287 
00288 TEMPLATE
00289 template <class INVEC, class OUTVEC>
00290   void THIS::scatter_hypre_vector(INVEC& iv,
00291                                   OUTVEC& ov) const {
00303 
00304   const int& n_rank = _forest->n_rank();
00305   int n_ldof = n_local_dof();
00306   std::vector<std::list<std::pair<int,int> > > rdof(n_rank);
00307   for (int i = 0;i < n_ldof;++ i) {
00308     if (is_dof_on_primary_geometry(i)) continue;
00309     int remote_rank = dof_primary_rank(i);
00310     rdof[remote_rank].push_back(std::pair<int,int>(i, (*this)(i)));
00311   }
00312 
00314   std::vector<int> target(n_rank);
00315   std::vector<BinaryBuffer<> > obuf(n_rank), ibuf(n_rank);
00316   for (int i = 0;i < n_rank;++ i) {
00317     target[i] = i;
00318     AFEPack::ostream<> os(obuf[i]);
00319 
00320     u_int n_item = rdof[i].size();
00321     os << n_item;
00322     typename std::list<std::pair<int,int> >::iterator
00323       the_dof = rdof[i].begin(),
00324       end_dof = rdof[i].end();
00325     for (;the_dof != end_dof;++ the_dof) {
00326       os << the_dof->first << the_dof->second;
00327     }
00328   }
00329 
00331   MPI_Barrier(_forest->communicator());
00332   sendrecv_data(_forest->communicator(), n_rank, obuf.begin(),
00333                 ibuf.begin(), target.begin());
00334 
00336   obuf.clear();
00337   obuf.resize(n_rank);
00338   for (int i = 0;i < n_rank;++ i) {
00339     AFEPack::istream<> is(ibuf[i]);
00340     AFEPack::ostream<> os(obuf[i]);
00341 
00342     u_int n_item;
00343     is >> n_item;
00344     os << n_item;
00345     if (n_item > 0) {
00346       std::vector<int> ldof(n_item);
00347       std::vector<int> gdof(n_item);
00348       std::vector<double> value(n_item);
00349       for (u_int j = 0;j < n_item;++ j) {
00350         is >> ldof[j] >> gdof[j];
00351       }
00352       HYPRE_IJVectorGetValues(iv, n_item, &gdof[0], &value[0]);
00353       for (u_int j = 0;j < n_item;++ j) {
00354         os << ldof[j] << value[j];
00355       }
00356     }
00357   }
00358 
00360   ibuf.clear();
00361   ibuf.resize(n_rank);
00362   MPI_Barrier(_forest->communicator());
00363   sendrecv_data(_forest->communicator(), n_rank, obuf.begin(),
00364                 ibuf.begin(), target.begin());
00365 
00367   for (int i = 0;i < n_rank;++ i) {
00368     AFEPack::istream<> is(ibuf[i]);
00369 
00370     u_int n_item;
00371     is >> n_item;
00372 
00373     int ldof;
00374     double value;
00375 
00376     for (u_int j = 0;j < n_item;++ j) {
00377       is >> ldof >> value;
00378       ov(ldof) = value;
00379     }
00380   }
00381 
00386   for (u_int i = 0, j = 0;i < n_ldof;++ i) {
00387     if (is_dof_on_primary_geometry(i)) {
00388       HYPRE_IJVectorGetValues(iv, 1, &(*this)(i), &ov(i));
00389     }
00390   }
00391 }
00392 
00393 #undef THIS
00394 #undef TEMPLATE 
00395