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