MeshKit
1.0
|
00001 #include "meshkit/ParallelMesher.hpp" 00002 #include "meshkit/ModelEnt.hpp" 00003 #include "meshkit/ParExchangeMesh.hpp" 00004 00005 #include "MKDefines.h" 00006 #include "MBTagConventions.hpp" 00007 #include "MBEntityType.h" 00008 #include "moab/Range.hpp" 00009 #include "moab/Types.hpp" 00010 00011 #include "CGMParallelComm.hpp" 00012 #include "RefEntity.hpp" 00013 #include "RefVolume.hpp" 00014 #include "RefFace.hpp" 00015 #include "RefEdge.hpp" 00016 #include "RefVertex.hpp" 00017 #include "CADefines.hpp" 00018 #include "GeometryQueryTool.hpp" 00019 00020 const bool debug_parallelmesher = false; 00021 00022 namespace MeshKit 00023 { 00024 // static registration of this mesh scheme 00025 moab::EntityType ParallelMesher_tps[] = {moab::MBVERTEX, moab::MBTRI, moab::MBTET, moab::MBMAXTYPE}; 00026 const moab::EntityType* ParallelMesher::output_types() 00027 { return ParallelMesher_tps; } 00028 00029 ParallelMesher::ParallelMesher(MKCore *mkcore, const MEntVector &me_vec) 00030 : MeshScheme(mkcore, me_vec) 00031 { 00032 // get information related to MOAB parallel communication 00033 m_mpcomm = moab::ParallelComm::get_pcomm(mk_core()->moab_instance(), 0); 00034 if (NULL == m_mpcomm) m_mpcomm = new moab::ParallelComm(mk_core()->moab_instance(), MPI_COMM_WORLD); 00035 m_rank = m_mpcomm->proc_config().proc_rank(); 00036 00037 // create tag 00038 iMesh::Error err = mk_core()->imesh_instance()->createTag("PARALLEL_UNIQUE_ID", 1, iBase_INTEGER, m_mPuniqueIDTag); 00039 IBERRCHK(err, "Trouble create a parallel unique id tag handle."); 00040 00041 m_sEntity.resize(10); 00042 00043 if (debug_parallelmesher) m_mpcomm->set_debug_verbosity(5); 00044 } 00045 00046 ParallelMesher::~ParallelMesher() 00047 { 00048 delete m_mpcomm; 00049 } 00050 00051 MeshOp *ParallelMesher::get_mesher(PARALLEL_OP_TYPE type) 00052 { 00053 MeshOpProxy* proxy = NULL; 00054 std::vector<MeshOpProxy *> proxies; 00055 if (type == MESH_VERTEX) { 00056 mk_core()->meshop_by_mesh_type(moab::MBVERTEX, proxies); 00057 proxy = *proxies.begin(); } 00058 else if (type == MESH_EDGE) { 00059 mk_core()->meshop_by_mesh_type(moab::MBEDGE, proxies); 00060 proxy = *proxies.begin(); 00061 } 00062 else if (type == MESH_INTER_SURF || type == MESH_NINTER_SURF) { 00063 //mk_core()->meshop_by_mesh_type(moab::MBTRI, proxies); 00064 proxy = MKCore::meshop_proxy("CAMALTriAdvance"); 00065 } 00066 else if (type == MESH_VOLUME) { 00067 //mk_core()->meshop_by_mesh_type(moab::MBTET, proxies); 00068 proxy = MKCore::meshop_proxy("CAMALTetMesher"); 00069 } 00070 else if (type == EXCHANGE_VERTEX || type == EXCHANGE_EDGE) { 00071 proxy = MKCore::meshop_proxy("ParExchangeMesh"); 00072 } 00073 else if (type == POST_RECV) { 00074 proxy = MKCore::meshop_proxy("ParPostRecv"); 00075 } 00076 else if (type == SEND_POST_SURF_MESH) { 00077 proxy = MKCore::meshop_proxy("ParSendPostSurfMesh"); 00078 } 00079 else if (type == RECV_SURF_MESH) { 00080 proxy = MKCore::meshop_proxy("ParRecvSurfMesh"); 00081 } 00082 00083 if (proxy == NULL) throw Error(MK_FAILURE, "Couldn't find a MeshOp capable of producing the given mesh type."); 00084 00085 return mk_core()->construct_meshop(proxy); 00086 } 00087 00088 void ParallelMesher::setup_this() 00089 { 00090 // make partitioned and send and receive model entity sets 00091 for (MEntSelection::iterator mit = mentSelection.begin(); 00092 mit != mentSelection.end(); mit++) { 00093 ModelEnt *me_vol = (*mit).first; 00094 if (me_vol->dimension() != 3) throw Error(MK_BAD_INPUT, "Parallel mesher assigned to an entity with dimension != 3."); 00095 00096 RefEntity* vol = reinterpret_cast<RefEntity*> (me_vol->geom_handle()); 00097 TDParallel *td_par_vol = (TDParallel *) vol->get_TD(&TDParallel::is_parallel); 00098 if (td_par_vol == NULL) ECERRCHK(MK_FAILURE, "Volume should have partitioned information."); 00099 00100 unsigned int charge_proc = td_par_vol->get_charge_proc(); 00101 MEntVector children; 00102 00103 // partititoned vols 00104 if (m_rank == charge_proc) { 00105 m_sEntity[MESH_VOLUME].insert(me_vol); // volume 00106 if (debug_parallelmesher) print_geom_info(me_vol, 3, true); 00107 00108 // get non-interface surfaces 00109 children.clear(); 00110 me_vol->get_adjacencies(2, children); 00111 for (MEntVector::iterator cit = children.begin(); cit != children.end(); cit++) { 00112 RefEntity* child = reinterpret_cast< RefEntity* > ((*cit)->geom_handle()); 00113 TDParallel *td_par_child = (TDParallel *) child->get_TD(&TDParallel::is_parallel); 00114 if (td_par_child == NULL) { 00115 m_sEntity[MESH_NINTER_SURF].insert(*cit); 00116 } 00117 } 00118 } 00119 00120 // get child interface entities 00121 for (int i = 2; i > -1; i--) { 00122 children.clear(); 00123 me_vol->get_adjacencies(i, children); 00124 00125 for (MEntVector::iterator cit = children.begin(); cit != children.end(); cit++) { 00126 RefEntity* child = reinterpret_cast< RefEntity* > ((*cit)->geom_handle()); 00127 TDParallel *td_par_child = (TDParallel *) child->get_TD(&TDParallel::is_parallel); 00128 if (td_par_child != NULL) { 00129 check_partition(td_par_child, *cit, i); 00130 if (debug_parallelmesher) print_geom_info(*cit, i, m_rank == td_par_child->get_charge_proc()); 00131 } 00132 } 00133 } 00134 } 00135 00136 add_parallel_mesh_op(MESH_VERTEX); 00137 add_parallel_mesh_op(EXCHANGE_VERTEX); 00138 add_parallel_mesh_op(MESH_EDGE); 00139 add_parallel_mesh_op(EXCHANGE_EDGE); 00140 add_parallel_mesh_op(MESH_INTER_SURF); 00141 add_parallel_mesh_op(SEND_POST_SURF_MESH); 00142 add_parallel_mesh_op(MESH_NINTER_SURF); 00143 add_parallel_mesh_op(RECV_SURF_MESH); 00144 add_parallel_mesh_op(MESH_VOLUME); 00145 00146 if (debug_parallelmesher) { 00147 for (PARALLEL_OP_TYPE type = MESH_VERTEX; type <= MESH_VOLUME;) { 00148 std::cout << "# of parallel mesh type " << type << "=" 00149 << m_sEntity[type].size() << std::endl; 00150 type = static_cast<PARALLEL_OP_TYPE> (type + 1); 00151 } 00152 } 00153 } 00154 00155 void ParallelMesher::check_partition(TDParallel* td_par, ModelEnt* me, int dim) 00156 { 00157 bool b_partitioned = false; 00158 unsigned int charge_p = td_par->get_charge_proc(); 00159 if (m_rank == charge_p) { // charge processor 00160 if (dim == 2) { 00161 m_sEntity[MESH_INTER_SURF].insert(me); 00162 m_sEntity[SEND_POST_SURF_MESH].insert(me); 00163 m_sEntity[RECV_SURF_MESH].insert(me); 00164 } 00165 else if (dim == 1) { 00166 m_sEntity[MESH_EDGE].insert(me); 00167 m_sEntity[EXCHANGE_EDGE].insert(me); 00168 } 00169 else if (dim == 0) { 00170 m_sEntity[MESH_VERTEX].insert(me); 00171 m_sEntity[EXCHANGE_VERTEX].insert(me); 00172 } 00173 b_partitioned = true; 00174 } 00175 else { 00176 DLIList<int>* shared_procs = td_par->get_shared_proc_list(); 00177 int n_shared = shared_procs->size(); 00178 for (int i = 0; i < n_shared; i++) { 00179 if (m_rank == (unsigned int)shared_procs->get_and_step()) { // shared processor 00180 if (dim == 2) { 00181 m_sEntity[SEND_POST_SURF_MESH].insert(me); 00182 m_sEntity[RECV_SURF_MESH].insert(me); 00183 } 00184 else if (dim == 1) { 00185 m_sEntity[EXCHANGE_EDGE].insert(me); 00186 } 00187 else if (dim == 0) { 00188 m_sEntity[EXCHANGE_VERTEX].insert(me); 00189 } 00190 b_partitioned = true; 00191 break; 00192 } 00193 } 00194 } 00195 00196 // set geometry unique id to corresponding surface entityset 00197 if (b_partitioned) { 00198 unsigned int unique_id = td_par->get_unique_id(); 00199 iBase_EntitySetHandle entityset = reinterpret_cast<iBase_EntitySetHandle> ((me)->mesh_handle()); 00200 iMesh::Error err = mk_core()->imesh_instance()->setEntSetIntData(entityset, m_mPuniqueIDTag, unique_id); 00201 IBERRCHK(err, "Couldn't set geometry unique id to corresponding mesh entityset."); 00202 } 00203 } 00204 00205 void ParallelMesher::print_geom_info(ModelEnt* me, const int dim, 00206 const bool local) 00207 { 00208 RefEntity* ent = reinterpret_cast<RefEntity*> (me->geom_handle()); 00209 CubitVector edge_geom = ent->center_point(); 00210 moab::EntityHandle entityset = me->mesh_handle(); 00211 TDParallel *td_par = (TDParallel *) ent->get_TD(&TDParallel::is_parallel); 00212 std::cout << "Partitioned_dim=" << dim; 00213 if (local) std::cout << ",local"; 00214 else std::cout << ",remote"; 00215 std::cout << ",geom_center=" << edge_geom.x() 00216 << " "<< edge_geom.y() << " " << edge_geom.z() << ", meshset=" 00217 << entityset << ",uid=" << td_par->get_unique_id() 00218 << ",id=" << ent->id() << std::endl; 00219 } 00220 00221 void ParallelMesher::add_parallel_mesh_op(PARALLEL_OP_TYPE type, bool after) 00222 { 00223 bool inserted = false; 00224 MeshOp *mesher = NULL; 00225 std::vector<MeshOp*> meshops; 00226 MEntSet::iterator iter = m_sEntity[type].begin(); 00227 MEntSet::iterator end_iter = m_sEntity[type].end(); 00228 00229 for (; iter != end_iter; iter++) { 00230 bool found = false; 00231 if (!mesher) mesher = get_mesher(type); // get a mesher if we haven't already 00232 meshops.clear(); 00233 (*iter)->get_meshops(meshops); 00234 int n_meshops = meshops.size(); 00235 00236 if (n_meshops > 0) { 00237 for (int i = 0; i < n_meshops; i++) { 00238 if (meshops[i] == mesher) { 00239 found = true; 00240 break; 00241 } 00242 } 00243 } 00244 00245 if (!found) { // if no specified meshop 00246 // add this entity to it, and if first, make sure it's added upstream 00247 mesher->add_modelent(*iter); 00248 if (!inserted) { 00249 if (after) mk_core()->insert_node(mesher, mk_core()->leaf_node(), this); 00250 else mk_core()->insert_node(mesher, this); 00251 inserted = true; 00252 } 00253 } 00254 } 00255 } 00256 00257 void ParallelMesher::execute_this() 00258 { 00259 } 00260 00261 void ParallelMesher::print_mesh() 00262 { 00263 moab::ErrorCode rval; 00264 int tmp_procs[MAX_SHARING_PROCS]; 00265 moab::EntityHandle tmp_hs[MAX_SHARING_PROCS]; 00266 unsigned char pstat; 00267 int num_ps; 00268 moab::Range entities; 00269 00270 for (moab::EntityType type = MBVERTEX; type != MBMAXTYPE; type++) { 00271 entities.clear(); 00272 rval = mk_core()->moab_instance()->get_entities_by_type(0, type, entities); 00273 MBERRCHK(rval, mk_core()->moab_instance()); 00274 00275 for (moab::Range::iterator rit = entities.begin(); rit != entities.end(); rit++) { 00276 rval = m_mpcomm->get_sharing_data(*rit, tmp_procs, tmp_hs, pstat, num_ps); 00277 MBERRCHK(rval, mk_core()->moab_instance()); 00278 00279 if (type != MBENTITYSET) { 00280 std::cout << "ParallelMesher::entity=" << *rit << ", type=" << type; 00281 if (type == MBVERTEX) { 00282 double coord[3]; 00283 rval = mk_core()->moab_instance()->get_coords(&(*rit), 1, coord); 00284 MBERRCHK(rval, mk_core()->moab_instance()); 00285 std::cout << ", coords=" << coord[0] << " " << coord[1] << " " << coord[2]; 00286 } 00287 else { 00288 std::vector<moab::EntityHandle> connect; 00289 rval = mk_core()->moab_instance()->get_connectivity(&(*rit), 1, connect); 00290 MBERRCHK(rval, mk_core()->moab_instance()); 00291 int n_conn = connect.size(); 00292 std::cout << ", connect="; 00293 for (int j = 0; j < n_conn; j++) { 00294 std::cout << connect[j] << " "; 00295 } 00296 } 00297 00298 std::cout << ", shared_info="; 00299 for (int ii = 0; ii < num_ps; ii++) { 00300 std::cout << tmp_procs[ii] << ":" << tmp_hs[ii] << ", "; 00301 } 00302 std::cout << std::endl; 00303 } 00304 } 00305 } 00306 } 00307 } // namespace MeshKit