|
MeshKit
1.0
|
|
MeshKit
1.0
|
00001 #include "meshkit/CoreGen.hpp" 00002 #include "moab/verdict/VerdictWrapper.hpp" 00003 #include "moab/NestedRefine.hpp" 00004 00005 #define ERROR(a) {if (iBase_SUCCESS != err) std::cerr << a << std::endl;} 00006 #define ERRORR(a,b) {if (iBase_SUCCESS != err) {std::cerr << a << std::endl; return b;}} 00007 namespace MeshKit 00008 { 00009 // static registration of this mesh scheme 00010 moab::EntityType CoreGen_tps[] = { moab::MBVERTEX, 00011 moab::MBEDGE, 00012 moab::MBTRI, 00013 moab::MBHEX, 00014 moab::MBMAXTYPE}; 00015 const moab::EntityType* CoreGen::output_types() 00016 { return CoreGen_tps; } 00017 00018 CoreGen::CoreGen( MKCore *mk, const MEntVector &me_vec) 00019 : MeshScheme( mk, me_vec), 00020 igeom(mk->igeom_instance()), imesh(mk->imesh_instance()), 00021 mb (mk->moab_instance()) 00022 { 00023 err = 0; 00024 run_flag = 1; 00025 UNITCELL_DUCT = 0; 00026 ASSY_TYPES = 1; 00027 pack_type = 1; 00028 symm = 1; 00029 z_height = 1; 00030 z_divisions = 2; 00031 set_DIM = 3; // default is 3D 00032 PII = acos(-1.0); 00033 comment = "!"; 00034 MAXCHARS = 2000; 00035 compute_meshtogeom = false; 00036 extrude_flag = false; 00037 mem_tflag = false; 00038 global_ids = true; 00039 merge_tol = 1.0e-4; 00040 do_merge = 1; 00041 update_sets = 0; 00042 merge_tag = NULL; 00043 nst_flag = false; 00044 nsb_flag = false; 00045 nss_flag = false; 00046 nssall_flag = false; 00047 have_hex27 = false; 00048 umr_flag = false; 00049 nst_Id = 9997; 00050 nsb_Id = 9998; 00051 prob_type = "mesh"; 00052 savefiles = "one"; 00053 num_nsside = 0; 00054 linenumber = 0; 00055 info = "off"; 00056 minfo = "off"; 00057 nringsx = 0; 00058 nringsy = 0; 00059 nDegree = 0; 00060 00061 // initialize more memory/time related variables 00062 ctload = 0, ctcopymove = 0, ctmerge = 0, ctextrude = 0, ctns = 0, ctgid = 0, ctsave = 0; 00063 tload = 0, tcopymove = 0, tmerge = 0, textrude = 0, tns = 0, tgid = 0, tsave = 0; 00064 ld_t = 0, ld_tload = 0, ld_tcopymove = 0, ld_tsave = 0, ld_tgid = 0, ld_tmerge = 0, ld_tns = 0; 00065 mem1 = 0, mem2 = 0, mem3 = 0, mem4 = 0, mem5 = 0, mem6 = 0, mem7 = 0; 00066 00067 00068 } 00069 00070 CoreGen::~CoreGen() 00071 {} 00072 00073 bool CoreGen::add_modelent(ModelEnt *model_ent) 00074 { 00075 return MeshOp::add_modelent(model_ent); 00076 } 00077 00078 void CoreGen::setup_this() 00079 { 00080 if(rank == 0) 00081 logfile << "Setting-up in CoreGen meshop.." << std::endl; 00082 if(run_flag != 0){ 00083 double ctload = 0; 00084 clock_t tload = 0; 00085 CClock ld_time; 00086 int ld_t = 0, ld_tload = 0; 00087 unsigned long long mem1 = 0; 00088 if (prob_type == "mesh") { 00089 if (procs == 1) { 00090 err = load_meshes(); 00091 if(err !=0) {logfile << "load meshes failed!" << std::endl; exit(1);} 00092 } 00093 else { 00094 #ifdef USE_MPI 00095 if(procs < (int) files.size()){ 00096 err = load_meshes_parallel(rank, procs); 00097 if(err !=0) {logfile << "failed to load meshes in parallel!" << std::endl; exit(1);} 00098 } 00099 else{ 00100 // if there are more procs than files distribute the copy/move work on each proc 00101 err = distribute_mesh(rank, procs); 00102 if(err !=0) {logfile << "distribute meshes failed!" << std::endl; exit(1);} 00103 00104 MPI_Barrier(MPI_COMM_WORLD); 00105 00106 err = load_meshes_more_procs(rank, procs); 00107 if(err !=0) {logfile << "load m meshes failed!" << std::endl; exit(1);} 00108 } 00109 //Get a pcomm object 00110 pc = new moab::ParallelComm(mk_core()->moab_instance(), MPI_COMM_WORLD, &err); 00111 #endif 00112 } 00113 00114 mk_core()->moab_instance()->estimated_memory_use(0, 0, 0, &mem1); 00115 ld_tload = ld_time.DiffTime(); 00116 tload = clock(); 00117 ctload = (double) (tload - sTime)/(60*CLOCKS_PER_SEC); 00118 00119 if (mem_tflag == true && rank == 0) { 00120 logfile << "\n" << " Clock time taken to load mesh files = " << ld_tload 00121 << " seconds" << std::endl; 00122 logfile << " CPU time = " << ctload << " mins" << std::endl; 00123 logfile << " Memory used: " << mem1/1e6 << " Mb\n For rank 0\n" << std::endl; 00124 } 00125 } 00126 00127 /*********************************************/ 00128 // load geometry files 00129 /*********************************************/ 00130 else if (prob_type == "geometry" && procs == 1) { 00131 err = load_geometries(); 00132 if(err !=0) {logfile << "load geometry failed!" << std::endl; exit(1);} 00133 } 00134 else if(prob_type == "geometry" && procs > 1){ 00135 logfile << " Parallel mode not supported for problem-type: Geometry " << std::endl; 00136 exit(1); 00137 } 00138 00139 /*********************************************/ 00140 // copy move 00141 /*********************************************/ 00142 CClock ld_cm; 00143 err = copymove(rank, procs); 00144 if (prob_type == "mesh"){ 00145 mk_core()->moab_instance()->estimated_memory_use(0, 0, 0, &mem2); 00146 ld_tcopymove = ld_cm.DiffTime(); 00147 tcopymove = clock(); 00148 ctcopymove = (double) (tcopymove - tload)/(60*CLOCKS_PER_SEC); 00149 00150 if (mem_tflag == true && (strcmp(prob_type.c_str(), "mesh") == 0) && rank == 0) { 00151 logfile << "\n" << " Clock time taken to copy/move mesh files = " << ld_tcopymove 00152 << " seconds" << std::endl; 00153 logfile << " CPU time = " << ctcopymove << " mins" << std::endl; 00154 logfile << " Memory used: " << mem2/1e6 << " Mb\n For rank 0\n" << std::endl; 00155 } 00156 } 00157 00158 for (unsigned int i = 0; i < assys.size(); i++) { 00159 if(prob_type =="mesh") 00160 cm[i]->setup_called(true); 00161 if(prob_type =="geometry") 00162 cg[i]->setup_called(true); 00163 } 00164 00165 if (prob_type == "mesh") { 00166 /*********************************************/ 00167 // merge 00168 /*********************************************/ 00169 CClock ld_mm; 00170 if (procs == 1){ 00171 // merge mesh now 00172 //std::vector<iBase_EntityHandle> ents; 00173 moab::Range ents; 00174 //int dim = imesh->getGeometricDimension(); 00175 mb->get_entities_by_dimension(0, set_DIM, ents); 00176 logfile << " Merging nodes.."<< std::endl; 00177 moab::MergeMesh mm(mb); 00178 moab::ErrorCode err = mm.merge_entities(ents, merge_tol, true); 00179 if (moab::MB_SUCCESS != err) { 00180 std::cerr << "Error in MergeMesh during merging entities" << std::endl; 00181 exit(2); 00182 } 00183 } 00184 else if(procs > 1){ 00185 if (rank == 0) { 00186 logfile << "Merging nodes in parallel. " << std::endl; 00187 } 00188 00189 #ifdef USE_MPI 00190 //Call the resolve parallel function 00191 moab::ParallelMergeMesh pm(pc, merge_tol); 00192 err = pm.merge(); 00193 if (err != moab::MB_SUCCESS) { 00194 std::cerr << "Merge Failed" << std::endl; 00195 //MPI_Abort(MPI_COMM_WORLD, 1); 00196 } 00197 #endif 00198 } 00199 mk_core()->moab_instance()->estimated_memory_use(0, 0, 0, &mem3); 00200 ld_tmerge = ld_mm.DiffTime(); 00201 tmerge = clock(); 00202 ctmerge = (double) (tmerge - tcopymove)/(60*CLOCKS_PER_SEC); 00203 00204 if (mem_tflag == true && rank == 0 ) { 00205 logfile << "\n" << " Clock time taken to merge nodes = " << ld_tmerge 00206 << " seconds" << std::endl; 00207 logfile << " CPU time = " << ctmerge << " mins" << std::endl; 00208 logfile << " Memory used: " << mem3/1e6 << " Mb\n For rank 0\n" << std::endl; 00209 } 00210 #ifdef USE_MPI 00211 MPI_Barrier(MPI_COMM_WORLD); 00212 #endif 00213 /*********************************************/ 00214 // extrude 00215 /*********************************************/ 00216 if(procs == 1){ 00217 if (extrude_flag == true) { 00218 CClock ld_em; 00219 err = extrude(); 00220 if(err !=0) {logfile << "extrusion failed!" << std::endl; exit(1);} 00221 00222 mk_core()->moab_instance()->estimated_memory_use(0, 0, 0, &mem4); 00223 ld_t = ld_em.DiffTime(); 00224 textrude = clock(); 00225 ctextrude = (double) (textrude - tmerge)/(60*CLOCKS_PER_SEC); 00226 00227 if (mem_tflag == true && rank == 0) { 00228 logfile << "\n" << " Clock time taken to extrude = " << ld_t 00229 << " seconds" << std::endl; 00230 logfile << " CPU time = " << ctextrude << " mins" << std::endl; 00231 logfile << " Memory used: " << mem4/1e6 << " Mb\n For rank 0\n" 00232 << std::endl; 00233 } 00234 } 00235 } 00236 if(extrude_flag == true) 00237 em->setup_called(true); 00238 #ifdef USE_MPI 00239 MPI_Barrier(MPI_COMM_WORLD); 00240 #endif 00241 } 00242 } 00243 } 00244 00245 00246 void CoreGen::execute_this() 00247 { 00248 if(rank == 0) 00249 logfile << "Executing in CoreGen meshop.." << std::endl; 00250 00251 if(run_flag != 0 && prob_type != "geometry"){ 00252 for (unsigned int i = 0; i < assys.size(); i++) { 00253 cm[i]->execute_called(true); 00254 } 00255 if(extrude_flag == true) 00256 em->execute_called(true); 00257 /*********************************************/ 00258 // assign gids 00259 /*********************************************/ 00260 CClock ld_gid; 00261 mk_core()->moab_instance()->estimated_memory_use(0, 0, 0, &mem5); 00262 ld_tgid = ld_gid.DiffTime(); 00263 tgid = clock(); 00264 ctgid = (double) (tgid-tmerge)/(60*CLOCKS_PER_SEC); 00265 00266 if (mem_tflag == true && rank == 0) { 00267 logfile << "\n" << " Clock time taken to assign gids = " << ld_tgid 00268 << " seconds" << std::endl; 00269 logfile << " CPU time = " << ctgid << " mins" << std::endl; 00270 logfile << " Memory used: " << mem5/1e6 << " Mb\n For rank 0\n" << std::endl; 00271 } 00272 /*********************************************/ 00273 // create neumann sets on the core model 00274 /*********************************************/ 00275 if((nss_flag == true || nsb_flag == true 00276 || nst_flag == true) && procs == 1){ 00277 CClock ld_ns; 00278 err = create_neumannset(); 00279 if(err !=0) {logfile << "create NeumannSet failed!" << std::endl; exit(1);} 00280 00281 mk_core()->moab_instance()->estimated_memory_use(0, 0, 0, &mem6); 00282 ld_tns = ld_ns.DiffTime(); 00283 tns = clock(); 00284 ctns = (double) (tns-tgid)/(60*CLOCKS_PER_SEC); 00285 if (mem_tflag == true && rank == 0) { 00286 logfile << "\n" << " Clock time taken to create neumann sets = " << ld_tns 00287 << " seconds" << std::endl; 00288 logfile << " CPU time = " << ctns << " mins" << std::endl; 00289 logfile << " Memory used: " << mem6/1e6 << " Mb\n For rank 0\n" << std::endl; 00290 } 00291 } 00292 if(umr_flag == true){ 00293 #ifdef USE_MPI 00294 err = refine_coremodel(); 00295 #endif 00296 } 00297 } 00298 if (prob_type == "mesh") { 00299 /*********************************************/ 00300 // save 00301 /*********************************************/ 00302 CClock ld_sv; 00303 if (procs == 1) { 00304 err = save_mesh(); 00305 if(err !=0) {logfile << "save mesh failed!" << std::endl; exit(1);} 00306 } else { 00307 if(savefiles != "one" && (savefiles == "multiple" || savefiles == "both")){ 00308 err = save_mesh(rank); // uncomment to save the meshes with each proc 00309 if(err !=0) {logfile << "save mesh failed!" << std::endl; exit(1);} 00310 } 00311 if(savefiles != "multiple"){ 00312 #ifdef USE_MPI 00313 double write_time = MPI_Wtime(); 00314 err = save_mesh_parallel(rank, procs); 00315 if(err !=0) {logfile << "save parallel mesh failed!" << std::endl; exit(1);} 00316 write_time = MPI_Wtime() - write_time; 00317 if (rank == 0) 00318 logfile << "Parallel write time = " << write_time/60.0 << " mins" << std::endl; 00319 #endif 00320 } 00321 } 00322 00323 mk_core()->moab_instance()->estimated_memory_use(0, 0, 0, &mem7); 00324 ld_tsave = ld_sv.DiffTime(); 00325 tsave = clock(); 00326 ctsave = (double) (tsave - tgid)/(60*CLOCKS_PER_SEC); 00327 00328 if (mem_tflag == true && rank == 0 ) { 00329 logfile << "\n" << " Clock time taken to save = " << ld_tsave << " seconds" 00330 << std::endl; 00331 logfile << " CPU time = " << ctsave << " mins" << std::endl; 00332 logfile << " Memory used: " << mem7/1e6 << " Mb\n For rank 0\n" << std::endl; 00333 } 00334 } 00335 /*********************************************/ 00336 // geometry operations 00337 /*********************************************/ 00338 else if (prob_type == "geometry") { 00339 err = save_geometry(); 00340 if(err !=0) {logfile << "save geometry failed!" << std::endl; exit(1);} 00341 } 00342 00343 /*********************************************/ 00344 // print memory and timing if using mpi 00345 /*********************************************/ 00346 mem1/=1e6; 00347 mem2/=1e6; 00348 mem3/=1e6; 00349 mem5/=1e6; 00350 mem7/=1e6; 00351 00352 #ifdef USE_MPI 00353 unsigned long max_mem7 = 1.0; 00354 MPI_Reduce( &mem7, &max_mem7, 1, MPI_UNSIGNED_LONG, MPI_MAX, 0, MPI_COMM_WORLD); 00355 #endif 00356 00357 #ifdef USE_MPI 00358 if (mem_tflag == true) { 00359 00360 unsigned long max_mem1 = 1.0, max_mem2 = 1.0, max_mem3 = 1.0, max_mem5 = 1.0; 00361 00362 MPI_Reduce( &mem1, &max_mem1, 1, MPI_UNSIGNED_LONG, MPI_MAX, 0, MPI_COMM_WORLD); 00363 MPI_Reduce( &mem2, &max_mem2, 1, MPI_UNSIGNED_LONG, MPI_MAX, 0, MPI_COMM_WORLD); 00364 MPI_Reduce( &mem3, &max_mem3, 1, MPI_UNSIGNED_LONG, MPI_MAX, 0, MPI_COMM_WORLD); 00365 MPI_Reduce( &mem5, &max_mem5, 1, MPI_UNSIGNED_LONG, MPI_MAX, 0, MPI_COMM_WORLD); 00366 00367 double max_ctload = -1.0, max_ctcopymove = -1.0, max_ctgid = -1.0, max_ctsave = -1.0, max_ctmerge = -1.0; 00368 MPI_Reduce( &ctload, &max_ctload, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); 00369 MPI_Reduce( &ctcopymove, &max_ctcopymove, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); 00370 MPI_Reduce( &ctmerge, &max_ctmerge, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); 00371 MPI_Reduce( &ctgid, &max_ctgid, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); 00372 MPI_Reduce( &ctsave, &max_ctsave, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); 00373 00374 int max_tload = -1.0, max_tcopymove = -1.0, max_tgid = -1.0, max_tsave = -1.0, max_tmerge = -1.0; 00375 MPI_Reduce( &ld_tload, &max_tload, 1, MPI_INT, MPI_MAX, 0, MPI_COMM_WORLD); 00376 MPI_Reduce( &ld_tcopymove, &max_tcopymove, 1, MPI_INT, MPI_MAX, 0, MPI_COMM_WORLD); 00377 MPI_Reduce( &ld_tmerge, &max_tmerge, 1, MPI_INT, MPI_MAX, 0, MPI_COMM_WORLD); 00378 MPI_Reduce( &ld_tgid, &max_tgid, 1, MPI_INT, MPI_MAX, 0, MPI_COMM_WORLD); 00379 MPI_Reduce( &ld_tsave, &max_tsave, 1, MPI_INT, MPI_MAX, 0, MPI_COMM_WORLD); 00380 00381 if(rank == 0 && procs > 1){ 00382 logfile << "\nMAXIMUM TIME TAKEN OVER ALL PROCS\nCLOCK TIME:-"; 00383 logfile << "\n**r = " << rank<< " Time taken to load mesh files = " << max_tload 00384 << " secs" << std::endl; 00385 logfile << "***r = : " << rank<< " Memory used: " << max_mem1 << " Mb" << std::endl; 00386 00387 // copymove 00388 logfile << "\n**r = " << rank<< " Time taken to copy/move mesh files = " << max_tcopymove 00389 << " secs" << std::endl; 00390 logfile << "***r = " << rank<< " Memory used: " << max_mem2 << " Mb" << std::endl; 00391 00392 // merge 00393 logfile << "\n**r = " << rank<< " Time taken to merge nodes = " << max_tmerge 00394 << " secs" << std::endl; 00395 logfile << "***r = " << rank<< " Memory used: " << max_mem3 << " kb" << std::endl; 00396 00397 // assign gid 00398 logfile << "\n**r = " << rank<< " Time taken to assign gids = " << max_tgid 00399 << " secs" << std::endl; 00400 logfile << "*** r = " << rank<< " Memory used: " << max_mem5 << " Mb" << std::endl; 00401 00402 // save 00403 logfile << "\n**r = " << rank<< " Time taken to save = " << max_tsave << " secs" 00404 << std::endl; 00405 logfile << "***r = " << rank<< " Memory used: " << max_mem7 << " Mb" << std::endl; 00406 00407 // cpu times 00408 logfile << "\n CPU TIME:-\n" << " r = " << rank<< " Time taken to load mesh files = " << ctload 00409 << " mins" << std::endl; 00410 00411 logfile << " r = " << rank << " Time taken to copy/move files = " << ctcopymove 00412 << " mins" << std::endl; 00413 00414 logfile << " r = " << rank << " Time taken to merge = " << ctmerge 00415 << " mins" << std::endl; 00416 00417 logfile << " r = " << rank << " Time taken to assign gids = " << ctgid 00418 << " mins" << std::endl; 00419 00420 logfile << " r = " << rank << " Time taken to save mesh = " << ctsave 00421 << " mins" << std::endl; 00422 } 00423 } 00424 #endif 00425 00426 if (rank == 0) { 00427 Timer.GetDateTime(szDateTime); 00428 logfile << "\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX" << std::endl; 00429 logfile << "Ending at : " << szDateTime; 00430 logfile << "Elapsed wall clock time: " << Timer.DiffTime() 00431 << " seconds or " << (Timer.DiffTime()) / 60.0 << " mins\n"; 00432 00433 logfile << "Total CPU time used: " << (double) (clock() - sTime)/(CLOCKS_PER_SEC) << " seconds or " << 00434 (double) (clock() - sTime)/(60*CLOCKS_PER_SEC) 00435 << " mins" << std::endl; 00436 #ifdef USE_MPI 00437 logfile << "Maximum memory used by a processor: " << max_mem7 << " Mb" << std::endl; 00438 #endif 00439 if(procs == 1) 00440 logfile << "Maximum memory used: " << mem7 << " Mb" << std::endl; 00441 logfile << "Total processors used = " << procs << std::endl; 00442 logfile << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n" << std::endl; 00443 00444 } 00445 } 00446 #ifdef USE_MPI 00447 int CoreGen::save_mesh_parallel(const int nrank, const int numprocs) 00448 // ------------------------------------------------------------------------------------------- 00449 // Function: resolve shared entitie and save mesh file in parallel (hdf5 only) 00450 // Input: none 00451 // Output: none 00452 // ------------------------------------------------------------------------------------------- 00453 { 00454 // start saving mesh in parallel 00455 if (nrank == 0) { 00456 logfile << "Saving mesh file in parallel, starting to cleanup sets " << std::endl; 00457 } 00458 00459 // handle sets before saving - delete all unnessary sets - this would save a lot of save time 00460 moab::Tag mattag; 00461 mb->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, mattag ); 00462 moab::Range matsets, this_mat_ents; 00463 mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &mattag, 0, 1, matsets ); 00464 moab::Range::iterator m_it; 00465 moab::EntityHandle old_mat_set, new_mat_set; 00466 00467 for(m_it = matsets.begin(); m_it != matsets.end(); m_it++){ 00468 this_mat_ents.clear(); 00469 old_mat_set = *m_it; 00470 mb->get_entities_by_dimension(old_mat_set, set_DIM, this_mat_ents, true); 00471 int material_id; 00472 mb->tag_get_data(mattag, &old_mat_set, 1, &material_id); 00473 // create a new material set, fill with ents (directly) and set id tag 00474 mb->create_meshset(MESHSET_SET, new_mat_set); 00475 mb->add_entities(new_mat_set, this_mat_ents); 00476 mb->tag_set_data(mattag, &new_mat_set, 1, &material_id); 00477 mb->delete_entities(&old_mat_set, 1); 00478 } 00479 00480 if (nrank == 0) { 00481 logfile << "Done dealing with material sets " << std::endl; 00482 } 00483 00484 moab::Tag nstag; 00485 mb->tag_get_handle( "NEUMANN_SET", 1, MB_TYPE_INTEGER, nstag ); 00486 moab::Range nssets, this_neu_ents; 00487 mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &nstag, 0, 1, nssets ); 00488 moab::Range::iterator n_it; 00489 moab::EntityHandle old_neu_set, new_neu_set; 00490 00491 for(n_it = nssets.begin(); n_it != nssets.end(); n_it++){ 00492 this_neu_ents.clear(); 00493 old_neu_set = *n_it; 00494 mb->get_entities_by_dimension(old_neu_set, (set_DIM-1), this_neu_ents, true); 00495 int neumann_id; 00496 mb->tag_get_data(nstag, &old_mat_set, 1, &neumann_id); 00497 // create a new neumann set, fill with ents (directly) and set id tag 00498 mb->create_meshset(MESHSET_SET, new_neu_set); 00499 mb->add_entities(new_neu_set, this_neu_ents); 00500 mb->tag_set_data(nstag, &new_neu_set, 1, &neumann_id); 00501 mb->delete_entities(&old_neu_set, 1); 00502 } 00503 00504 00505 moab::Tag drtag; 00506 mb->tag_get_handle( "DIRICHLET_SET", 1, MB_TYPE_INTEGER, drtag ); 00507 moab::Range drsets, this_dir_ents; 00508 mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &drtag, 0, 1, drsets ); 00509 moab::Range::iterator d_it; 00510 moab::EntityHandle old_dir_set, new_dir_set; 00511 00512 for(d_it = drsets.begin(); d_it != drsets.end(); d_it++){ 00513 this_dir_ents.clear(); 00514 old_dir_set = *n_it; 00515 mb->get_entities_by_dimension(old_dir_set, 0, this_dir_ents, true); 00516 int dirichlet_id; 00517 mb->tag_get_data(drtag, &old_dir_set, 1, &dirichlet_id); 00518 // create a new neumann set, fill with ents (directly) and set id tag 00519 mb->create_meshset(MESHSET_SET, new_dir_set); 00520 mb->add_entities(new_dir_set, this_dir_ents); 00521 mb->tag_set_data(drtag, &new_dir_set, 1, &dirichlet_id); 00522 mb->delete_entities(&old_dir_set, 1); 00523 } 00524 00525 00526 if (nrank == 0) { 00527 logfile << "Done dealing with all ms ns and ds" << std::endl; 00528 } 00529 00530 00531 // Deleting GEOM_DIMENSION tags and others? 00532 moab::Range all_sets; 00533 moab::EntityHandle curr_set; 00534 mb->get_entities_by_type(0, MBENTITYSET, all_sets); 00535 moab::Range::iterator all_it; 00536 moab::Tag gdtag; 00537 mb->tag_get_handle( "GEOM_DIMENSION", 1, MB_TYPE_INTEGER, gdtag ); 00538 00539 for(all_it = all_sets.begin(); all_it != all_sets.end(); all_it++){ 00540 curr_set = *all_it; 00541 int temp_gid = -1; 00542 mb->tag_get_data(gdtag, &curr_set, 1, &temp_gid); 00543 if(temp_gid != -1){ 00544 mb->delete_entities(&curr_set, 1); 00545 } 00546 } 00547 00548 MPI_Barrier(MPI_COMM_WORLD); 00549 00550 if (nrank == 0) { 00551 logfile << "Done deleting gd sets, now starting to resolve shared ents " << std::endl; 00552 } 00553 00554 00555 // resolve shared sets to create only on MATERIAL_SET 00556 matsets.clear(); 00557 mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &mattag, 0, 1, matsets ); 00558 if(matsets.size() > 0) 00559 pc->resolve_shared_sets( matsets, mattag ); 00560 00561 if (nrank == 0) { 00562 logfile << matsets.size() << "Done resolving material ents " << std::endl; 00563 } 00564 00565 /*// resolve 00566 nssets.clear(); 00567 mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &nstag, 0, 1, nssets ); 00568 if(nssets.size() > 0) 00569 pc->resolve_shared_sets( nssets, nstag ); 00570 */ 00571 // resolve 00572 drsets.clear(); 00573 mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &drtag, 0, 1, drsets ); 00574 if(drsets.size() > 0) 00575 pc->resolve_shared_sets( drsets, drtag ); 00576 00577 // need this barrier before setting pp tag 00578 // MPI_Barrier(MPI_COMM_WORLD); 00579 00580 00581 // Done with deleting recursive sets now create pp tags and save 00582 if (nrank == 0) { 00583 logfile << "setting PARALLEL_PARTITION tag" << std::endl; 00584 logfile << "Saving mesh file in parallel. " << std::endl; 00585 } 00586 moab::Range entities; 00587 moab::EntityHandle meshsetp; 00588 mb->get_entities_by_type(0, MBHEX, entities); 00589 00590 mb->create_meshset(MESHSET_SET, meshsetp); 00591 mb->add_entities(meshsetp, entities); 00592 00593 moab::Tag pp_tag; 00594 00595 mb->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, pp_tag, MB_TAG_SPARSE|MB_TAG_CREAT); 00596 mb->tag_set_data(pp_tag, &meshsetp, 1, &nrank); 00597 00598 00599 //MPI_Barrier(MPI_COMM_WORLD); 00600 00601 00602 // flag specified in input file 00603 if(have_hex27 == true){ 00604 moab::Range entities; 00605 moab::EntityHandle meshset; 00606 mb->get_entities_by_type(0, MBHEX, entities); 00607 mb->create_meshset(MESHSET_SET, meshset); 00608 mb->add_entities(meshset, entities); 00609 // Add nodes along mid- edge, face and region 00610 mb->convert_entities(meshset, true, true, true); 00611 } 00612 /* 00613 moab::Range out_sets; 00614 out_sets.merge(matsets); 00615 out_sets.merge(nssets); 00616 out_sets.merge(drsets); 00617 out_sets.insert(meshsetp); 00618 */ 00619 00620 MPI_Barrier(MPI_COMM_WORLD); 00621 if (nrank == 0) { 00622 logfile << "Before saving mesh file in parallel. " << std::endl; 00623 } 00624 00625 MPI_Barrier(MPI_COMM_WORLD); 00626 00627 moab::ErrorCode rval = mb->write_file(outfile.c_str() , 0,"PARALLEL=WRITE_PART;CPUTIME;"/*DEBUG_IO=2;", out_sets*/); 00628 if(rval != moab::MB_SUCCESS) { 00629 std::cerr<<"Writing output file failed Code:"; 00630 std::string foo = ""; mb->get_last_error(foo); 00631 std::cerr<<"File Error: "<<foo<<std::endl; 00632 return 1; 00633 } 00634 if (nrank == 0) { 00635 logfile << "Done saving mesh file: " << outfile << std::endl; 00636 } 00637 return iBase_SUCCESS; 00638 } 00639 #endif 00640 int CoreGen::save_mesh(int nrank) { 00641 // --------------------------------------------------------------------------- 00642 // Function: save mesh serially from each rank 00643 // Input: none 00644 // Output: none 00645 // --------------------------------------------------------------------------- 00646 00647 // set parallel partition tag 00648 moab::Range entities; 00649 moab::EntityHandle meshset; 00650 mb->get_entities_by_type(0, MBHEX, entities); 00651 00652 mb->create_meshset(MESHSET_SET, meshset); 00653 mb->add_entities(meshset, entities); 00654 00655 moab::Tag pp_tag; 00656 00657 mb->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, pp_tag, MB_TAG_SPARSE|MB_TAG_CREAT); 00658 logfile << "setting PARALLEL_PARTITION tag" << std::endl; 00659 mb->tag_set_data(pp_tag, &meshset, 1, &nrank); 00660 00661 00662 if(have_hex27 == true){ 00663 moab::Range entities; 00664 moab::EntityHandle meshset; 00665 mb->get_entities_by_type(0, MBHEX, entities); 00666 00667 mb->create_meshset(MESHSET_SET, meshset); 00668 mb->add_entities(meshset, entities); 00669 mb->convert_entities(meshset, true, true, true); 00670 } 00671 // export proc- nrank mesh 00672 std::ostringstream os; 00673 std::string fname; 00674 fname = iname; 00675 os << fname << nrank << ".h5m"; 00676 fname = os.str(); 00677 iMesh_save(imesh->instance(), root_set, fname.c_str(), NULL, &err, strlen(fname.c_str()), 0); 00678 ERRORR("Trouble writing output mesh.", err); 00679 logfile << "Saved mesh file: " << fname.c_str() << std::endl; 00680 00681 return iBase_SUCCESS; 00682 } 00683 00684 int CoreGen::save_mesh() { 00685 // --------------------------------------------------------------------------- 00686 // Function: save mesh serially 00687 // Input: none 00688 // Output: none 00689 // --------------------------------------------------------------------------- 00690 // export 00691 if(have_hex27 == true){ 00692 moab::Range entities; 00693 moab::EntityHandle meshset; 00694 mb->get_entities_by_type(0, MBHEX, entities); 00695 00696 mb->create_meshset(MESHSET_SET, meshset); 00697 mb->add_entities(meshset, entities); 00698 mb->convert_entities(meshset, true, true, true); 00699 } 00700 00701 mb->write_mesh(outfile.c_str()); 00702 logfile << "Saved mesh file: " << outfile.c_str() << std::endl; 00703 00704 return iBase_SUCCESS; 00705 } 00706 00707 00708 int CoreGen::save_geometry() { 00709 // --------------------------------------------------------------------------- 00710 // Function: save geometry serially 00711 // Input: none 00712 // Output: none 00713 // --------------------------------------------------------------------------- 00714 /* double dTol = 1e-3; 00715 00716 // getting all entities for merge and imprint 00717 SimpleArray<iBase_EntityHandle> entities_merge, entities_imprint; 00718 iGeom_getEntities(igeom->instance(), root_set, iBase_REGION, 00719 ARRAY_INOUT(entities_merge), &err ); 00720 ERRORR("Trouble writing output geometry.", err); 00721 00722 // merge and imprint before save 00723 logfile << "Merging.." << std::endl; 00724 iGeom_mergeEnts(igeom->instance(), ARRAY_IN(entities_merge), dTol, &err); 00725 ERRORR("Trouble writing output geometry.", err); 00726 00727 iGeom_getEntities( igeom->instance(), root_set, iBase_REGION, ARRAY_INOUT(entities_imprint),&err ); 00728 ERRORR("Trouble writing output geometry.", err); 00729 00730 // logfile << "Imprinting.." << std::endl; 00731 // iGeom_imprintEnts(igeom->instance(), ARRAY_IN(entities_imprint),&err); 00732 // ERRORR("Trouble writing output geometry.", err); 00733 // export 00734 logfile << "Saving geometry file: " << outfile << std::endl; 00735 */ 00736 iGeom_save(igeom->instance(), outfile.c_str(), NULL, &err, 00737 strlen(outfile.c_str()), 0); 00738 ERRORR("Trouble writing output geometry.", err); 00739 logfile << "Saved geometry file: "<< outfile.c_str() <<std::endl; 00740 00741 return iBase_SUCCESS; 00742 } 00743 00744 int CoreGen::distribute_mesh(const int nrank, int numprocs) 00745 // ------------------------------------------------------------------------------------------- 00746 // Function: merge the nodes within a set tolerance in the model 00747 // Input: none 00748 // Output: none 00749 // ------------------------------------------------------------------------------------------- 00750 { 00751 int nback = files.size() - nassys; 00752 if(nrank < ((int) core_alias.size() + nback)){ 00753 if(numprocs > (int) core_alias.size()){ 00754 numprocs = core_alias.size() + nback; 00755 } 00756 00757 rank_load.resize(numprocs); 00758 int extra_procs = numprocs - files.size(); 00759 if(numprocs >= (int) files.size() && numprocs <= (tot_assys + nback)){ 00760 // again fill assm_meshfiles 00761 for(int p=0; p<nassys; p++){ 00762 assm_meshfiles[p]=0; 00763 load_per_assm[p]=0; 00764 } 00765 for(int p=0; p<tot_assys; p++){ 00766 for(int q=0; q<nassys; q++){ 00767 if(strcmp(core_alias[p].c_str(), assm_alias[q].c_str()) ==0) { 00768 assm_meshfiles[q]+=1; 00769 load_per_assm[q]+=1; 00770 } 00771 } 00772 } 00773 if(nrank == 0){ 00774 logfile << " assm_meshfiles" << std::endl; 00775 for(int p=0; p<nassys; p++){ 00776 logfile << assm_meshfiles[p] << " : " << p << std::endl; 00777 } 00778 logfile << " load per assm" << std::endl; 00779 for(int p=0; p<nassys; p++){ 00780 logfile << load_per_assm[p] << " : " << p << std::endl; 00781 } 00782 } 00783 //distribute 00784 for(int i=0; i< (int)files.size(); i++){ 00785 rank_load[i] = i; 00786 if(i< nassys) 00787 times_loaded[i]+=1; 00788 } 00789 00790 double temp = 0; 00791 int assm_load = - 1; 00792 int e= 0; 00793 // compute the rank, mf vector for extra procs 00794 while(e < extra_procs){ 00795 for(int i = 0; i < nassys; i++){ 00796 if (load_per_assm[i] > temp ){ 00797 //if (load_per_assm[i] > temp){ 00798 temp = load_per_assm[i]; 00799 assm_load = i; 00800 } 00801 else if (assm_load == -1){ 00802 logfile << "No assemblies mesh files used in core" << std::endl; 00803 exit(0); 00804 } 00805 } 00806 //assm_meshfiles[assm_load]-=1; 00807 times_loaded[assm_load]+=1; 00808 load_per_assm[assm_load] = (double)assm_meshfiles[assm_load]/(double)times_loaded[assm_load]; 00809 00810 if(nrank == 0) 00811 logfile << assm_load <<": - assembly has a Current Load: " << load_per_assm[assm_load] << std::endl; 00812 int temp_rank = files.size()+ e; 00813 rank_load[temp_rank] = assm_load; 00814 e++; 00815 temp = 0; 00816 } 00817 } 00818 else{ 00819 logfile << "Warning: #procs <= #assys in core, some processor will be idle" << std::endl; 00820 } 00821 if(nrank == 0){ 00822 logfile << "Times loaded 1 " << std::endl; 00823 for(int p=0; p<nassys; p++){ 00824 logfile << times_loaded[p] << " : " << p << std::endl; 00825 } 00826 00827 } 00828 // times_loaded.resize(nassys); 00829 std::vector<std::vector<int> > meshfiles_rank (files.size()); 00830 for(int i=0; i < (int) files.size(); i++){ 00831 for(int j=0; j < (int) rank_load.size(); j++){ 00832 if(rank_load[j]==i){ 00833 meshfiles_rank[i].push_back(j); 00834 // already done above times_loaded[i]+=1; 00835 } 00836 } 00837 } 00838 00839 position_core.resize(numprocs); 00840 00841 for(int i=0; i < (int) files.size(); i++){ 00842 int k = 0; 00843 if(i < (nassys) ){ 00844 for(int j=0; j < (int) assm_location[i].size(); j++){ 00845 if (k >= (int) meshfiles_rank[i].size()){ 00846 k = 0; 00847 } 00848 int p = meshfiles_rank[i][k]; 00849 int q = assm_location[i][j]; 00850 position_core[p].push_back(q); 00851 00852 ++k; 00853 } 00854 } 00855 else{ 00856 // this is background mesh set it -2, no meshfile to copy/move 00857 position_core[i].push_back(-2); 00858 } 00859 } 00860 if(nrank == 0){ 00861 logfile << "Times loaded 1 After" << std::endl; 00862 for(int p=0; p<nassys; p++){ 00863 logfile << times_loaded[p] << " : " << p << std::endl; 00864 } 00865 00866 } 00867 if(nrank == 0){ 00868 logfile << "FINAL scheme 1 assm_meshfiles" << std::endl; 00869 for(int p=0; p<nassys; p++){ 00870 logfile << assm_meshfiles[p] << " : " << p << std::endl; 00871 } 00872 logfile << "FINAL scheme 1 load per assm" << std::endl; 00873 for(int p=0; p<nassys; p++){ 00874 logfile << load_per_assm[p] << " : " << p << std::endl; 00875 } 00876 } 00877 if(nrank == 0){ 00878 logfile << " copy/move task distribution " << std::endl; 00879 for(int i =0; i< numprocs; i++){ 00880 logfile << "rank: " << i << " positions : "; 00881 for(int j=0; j< (int) position_core[i].size(); j++){ 00882 logfile << (int) position_core[i][j] << " "; 00883 } 00884 logfile << "\n" << std::endl; 00885 } 00886 } 00887 } 00888 return 0; 00889 } 00890 00891 00892 int CoreGen::load_meshes_more_procs(const int nrank, int numprocs) 00893 // --------------------------------------------------------------------------- 00894 // Function: loads all the meshes and initializes copymesh object 00895 // Input: none 00896 // Output: none 00897 // --------------------------------------------------------------------------- 00898 { 00899 iMesh_getRootSet(imesh->instance(), &root_set, &err); 00900 ERRORR("Couldn't get the root set", err); 00901 00902 int temp_index = rank_load[nrank]; 00903 00904 iBase_EntitySetHandle orig_set; 00905 iMesh_createEntSet(imesh->instance(), 0, &orig_set, &err); 00906 ERRORR("Couldn't create file set.", err); 00907 00908 // load this file 00909 iMesh_load(imesh->instance(), orig_set, files[temp_index].c_str(), NULL, &err, strlen(files[temp_index].c_str()), 0); 00910 ERRORR("Couldn't read mesh file.", err); 00911 logfile << "Loaded mesh file " << temp_index << " in processor: " << nrank << std::endl; 00912 00913 if(bsameas[temp_index] == 0 && temp_index < nassys){ 00914 if(all_ms_starts[temp_index] != -1 && all_ns_starts[temp_index] !=-1){ 00915 if(!shift_mn_ids(orig_set, temp_index)) 00916 ERRORR("Couldn't shift material and neumann set id's.", 1); 00917 } 00918 } 00919 00920 assys.push_back(orig_set); 00921 assys_index.push_back(temp_index); 00922 00923 // create cm instances for each mesh file 00924 cm.resize(assys.size()); 00925 for (unsigned int i = 0; i < assys.size(); i++) { 00926 ModelEnt *me; 00927 me = NULL; 00928 me = new ModelEnt(mk_core(), iBase_EntitySetHandle(0), /*igeom instance*/0, (moab::EntityHandle)orig_set, 0); 00929 MEntVector assm_set; 00930 assm_set.push_back(me); 00931 cm[i] = (CopyMesh*) mk_core()->construct_meshop("CopyMesh", assm_set); 00932 cm[i]->set_name("copy_move_mesh"); 00933 } 00934 return iBase_SUCCESS; 00935 } 00936 00937 00938 00939 int CoreGen::load_meshes_parallel(const int nrank, int numprocs) 00940 // --------------------------------------------------------------------------- 00941 // Function: loads all the meshes and initializes copymesh object 00942 // Input: none 00943 // Output: none 00944 // --------------------------------------------------------------------------- 00945 { 00946 int nback = files.size() - nassys; 00947 cm.resize(files.size()); 00948 00949 iMesh_getRootSet(imesh->instance(), &root_set, &err); 00950 ERRORR("Couldn't get the root set", err); 00951 if(nrank < ((int) core_alias.size() + nback)){ 00952 if(numprocs > (int) core_alias.size()){ 00953 numprocs = core_alias.size() + nback; 00954 } 00955 00956 #ifdef USE_MPI 00957 if(numprocs > ((int) core_alias.size() + nback)){ 00958 logfile << "Warning: #procs <= #assys in core, some processor will be idle" << std::endl; 00959 } 00960 00961 iBase_EntitySetHandle orig_set; 00962 int temp_index; 00963 int extra_procs = numprocs - files.size(); 00964 std::vector<int> rank_load; 00965 for(int i = 0; i< (int) files.size(); i++){ 00966 temp_index = numprocs*i + nrank; 00967 if(temp_index >= (int) files.size()){ 00968 if (nrank >= (int) files.size() && nrank <= (tot_assys + nback)){ 00969 int temp = 1; 00970 int p = extra_procs; 00971 // compute the rank, mf vector for extra procs 00972 while(p !=0){ 00973 int assm_load = - 1; 00974 for(int i = 0; i < nassys; i++){ 00975 if ((int) assm_meshfiles[i] > temp){ 00976 temp = assm_meshfiles[i]; 00977 assm_load = i; 00978 } 00979 } 00980 if (assm_load == -1){ 00981 continue; 00982 logfile << "Warning: #procs <= #assys in core, some processor will be idle" << std::endl; 00983 } 00984 assm_meshfiles[assm_load]-=1; 00985 rank_load.push_back(assm_load); 00986 --p; 00987 temp = 1; 00988 } 00989 00990 temp_index = nrank - files.size(); 00991 iMesh_createEntSet(imesh->instance(), 0, &orig_set, &err); 00992 ERRORR("Couldn't create file set.", err); 00993 00994 if(!compute_meshtogeom){ 00995 // load this file 00996 iMesh_load(imesh->instance(), orig_set, files[rank_load[temp_index]].c_str(), NULL, &err, strlen(files[rank_load[temp_index]].c_str()), 0); 00997 ERRORR("Couldn't read mesh file.", err); 00998 } 00999 else{ 01000 load_and_compute_meshtogeom(orig_set, files[rank_load[temp_index]]); 01001 } 01002 logfile << "Loaded mesh file " << rank_load[temp_index] << " in processor: " << nrank << std::endl; 01003 01004 ModelEnt *me; 01005 me = NULL; 01006 me = new ModelEnt(mk_core(), iBase_EntitySetHandle(0), /*igeom instance*/0, (moab::EntityHandle)orig_set, 0); 01007 MEntVector assm_set; 01008 assm_set.push_back(me); 01009 cm[i] = (CopyMesh*) mk_core()->construct_meshop("CopyMesh", assm_set); 01010 cm[i]->set_name("copy_move_mesh"); 01011 01012 assys.push_back(orig_set); 01013 assys_index.push_back(rank_load[temp_index]); 01014 break; 01015 } 01016 } 01017 else{ 01018 iMesh_createEntSet(imesh->instance(), 0, &orig_set, &err); 01019 ERRORR("Couldn't create file set.", err); 01020 01021 if(!compute_meshtogeom){ 01022 // load this file 01023 iMesh_load(imesh->instance(), orig_set, files[temp_index].c_str(), NULL, &err, strlen(files[temp_index].c_str()), 0); 01024 ERRORR("Couldn't read mesh file.", err); 01025 } 01026 else{ 01027 load_and_compute_meshtogeom(orig_set, files[temp_index]); 01028 }\ 01029 logfile << "Loaded mesh file " << temp_index << " in processor: " << nrank << std::endl; 01030 01031 ModelEnt *me; 01032 me = NULL; 01033 me = new ModelEnt(mk_core(), iBase_EntitySetHandle(0), /*igeom instance*/0, (moab::EntityHandle)orig_set, 0); 01034 MEntVector assm_set; 01035 assm_set.push_back(me); 01036 cm[i] = (CopyMesh*) mk_core()->construct_meshop("CopyMesh", assm_set); 01037 cm[i]->set_name("copy_move_mesh"); 01038 01039 assys.push_back(orig_set); 01040 assys_index.push_back(temp_index); 01041 } 01042 } 01043 #endif 01044 } 01045 return iBase_SUCCESS; 01046 } 01047 01048 int CoreGen::load_meshes() 01049 // --------------------------------------------------------------------------- 01050 // Function: loads all the meshes and initializes copymesh and merge mesh objects 01051 // Input: none 01052 // Output: none 01053 // --------------------------------------------------------------------------- 01054 { 01055 // create cm instances for each mesh file 01056 cm.resize(files.size()); 01057 iBase_EntitySetHandle orig_set; 01058 // loop reading all mesh files 01059 for (unsigned int i = 0; i < files.size(); i++) { 01060 iMesh_createEntSet(imesh->instance(), 0, &orig_set, &err); 01061 ERRORR("Couldn't create file set.", err); 01062 logfile << "Loading File: " << files[i].c_str() << std::endl; 01063 if(!compute_meshtogeom){ 01064 iMesh_load(imesh->instance(), orig_set, files[i].c_str(), NULL, &err, strlen(files[i].c_str()), 0); 01065 ERRORR("Couldn't read mesh file.", err); 01066 mk_core()->populate_model_ents(0,0,0); 01067 } 01068 else{ 01069 load_and_compute_meshtogeom(orig_set, files[i]); 01070 } 01071 ModelEnt *me; 01072 me = NULL; 01073 me = new ModelEnt(mk_core(), iBase_EntitySetHandle(0), /*igeom instance*/0, (moab::EntityHandle)orig_set, 0); 01074 MEntVector assm_set; 01075 //assm_set.clear(); 01076 assm_set.push_back(me); 01077 cm[i] = (CopyMesh*) mk_core()->construct_meshop("CopyMesh", assm_set); 01078 cm[i]->set_name("copy_move_mesh"); 01079 cm[i]->copy_sets().add_set(orig_set); 01080 01081 //resize this else code will break 01082 //check if we've loaded the same mesh file and need to shift the Material and Neumann set start id's 01083 if(bsameas[i] == 0 && (int) i < nassys){ 01084 if(all_ms_starts[i] != -1 && all_ns_starts[i] !=-1){ 01085 if(!shift_mn_ids(orig_set, i)) 01086 ERRORR("Couldn't shift material and neumann set id's.", 1); 01087 } 01088 } 01089 01090 assys.push_back(orig_set); 01091 assys_index.push_back(i); 01092 } 01093 logfile << "Loaded mesh files." << std::endl; 01094 01095 return iBase_SUCCESS; 01096 } 01097 01098 int CoreGen::load_and_compute_meshtogeom(iBase_EntitySetHandle orig_set, std::string filename) 01099 // --------------------------------------------------------------------------- 01100 // Function: loads all the meshes and initializes copymesh and merge mesh objects 01101 // Input: none 01102 // Output: none 01103 // --------------------------------------------------------------------------- 01104 { 01105 bool isCub = false; 01106 std::string ext; 01107 std::string::size_type idx; 01108 idx = filename.rfind('.'); 01109 if(idx != std::string::npos) 01110 { 01111 ext = filename.substr(idx+1); 01112 } 01113 else 01114 { 01115 // No extension found 01116 } 01117 if(ext == "cub") isCub = true; 01118 01119 if(isCub){ 01120 mk_core()->igeom_instance()->deleteAll(); 01121 01122 iGeom_load(igeom->instance(), filename.c_str(), NULL, &err, 01123 strlen(filename.c_str()), 0); 01124 iMesh_load(imesh->instance(), orig_set, filename.c_str(), NULL, &err, strlen(filename.c_str()), 0); 01125 ERRORR("Couldn't read mesh file.", err); 01126 // populate model eneities here so we can call irel function from MeshKit 01127 01128 mk_core()->irel_pair()->inferAllRelations(); 01129 mk_core()->populate_model_ents(0,0,0); 01130 01131 01132 // now calculate geometric and mesh volume 01133 VerdictWrapper vw(mb); 01134 // Get all the materials in this mesh file 01135 moab::Tag mattag, gidtag; 01136 moab::Range matsets; 01137 mb->tag_get_handle("MATERIAL_SET", 1, MB_TYPE_INTEGER, mattag); 01138 mb->tag_get_handle("GLOBAL_ID", 1, MB_TYPE_INTEGER, gidtag); 01139 mb->get_entities_by_type_and_tag((moab::EntityHandle) orig_set, MBENTITYSET, &mattag, 0, 1, matsets ); 01140 01141 double mtot = 0.0, volume = 0.0; 01142 //loop thru all the materials 01143 moab::Range::iterator set_it; 01144 for (set_it = matsets.begin(); set_it != matsets.end(); set_it++) { 01145 moab::EntityHandle this_set = *set_it; 01146 01147 // get the id for this set 01148 int material_id; 01149 mb->tag_get_data(mattag, &this_set, 1, &material_id); 01150 // add up the geometric volume for these geometric volumes 01151 // get the entity sets there must be one entity set for each volume 01152 std::vector<moab::EntityHandle> geomsets_for_gid; 01153 geomsets_for_gid.clear(); 01154 mb->get_entities_by_type(this_set, MBENTITYSET, geomsets_for_gid); 01155 01156 for(unsigned int volid = 0; volid < geomsets_for_gid.size(); volid++){ 01157 // get the gid of this volume 01158 int my_geom_id = 0; 01159 mb->tag_get_data(gidtag, &geomsets_for_gid[volid], 1, &my_geom_id); 01160 iBase_EntityHandle ent2=NULL; 01161 mk_core()->irel_pair()->getSetEntRelation((iBase_EntitySetHandle) geomsets_for_gid[volid], 1, ent2); 01162 double myvol = 0.0; 01163 if(ent2 != NULL){ 01164 mk_core()->igeom_instance()->measure(&ent2,1,&myvol); 01165 volume+=myvol; 01166 } 01167 } 01168 if(volume == 0){ 01169 logfile << "Trying to compute MESHTOGEOM, but volume obtained is zero." << std::endl; 01170 logfile << "Geometry engine must match file type. Aborting." << std::endl; 01171 // exit(0); don't abort, continue to run. MESHTOGEOM won't be computed rest should still work. 01172 } 01173 01174 std::vector<moab::EntityHandle> elems; 01175 mb->get_entities_by_dimension(this_set, 3, elems, true); 01176 // get all the elements in this material 01177 double mvol_temp = 0.0; 01178 for(unsigned int i=0; i<elems.size();i++){ 01179 mvol_temp = 0.0; 01180 vw.quality_measure(elems[i], moab::MB_VOLUME, mvol_temp); 01181 mtot+=mvol_temp; 01182 } 01183 01184 double meshtogeom = mtot/volume; 01185 meshtogeom_file << material_id << " "<< meshtogeom << std::endl; 01186 logfile << material_id << " " << meshtogeom << std::endl; 01187 moab::Tag mtog_tag; 01188 // now set the meshtogeom tag on this set 01189 mb->tag_get_handle( "MESHTOGEOM", 1, MB_TYPE_DOUBLE, 01190 mtog_tag, MB_TAG_SPARSE|MB_TAG_CREAT ); 01191 mb->tag_set_data(mtog_tag, &(*set_it), 1, &meshtogeom); 01192 01193 elems.clear(); 01194 mtot = 0.0; 01195 volume = 0.0; 01196 } 01197 // 01198 } 01199 else{ 01200 logfile << "Trying to load a file that doesn't have geometry: Aborting" << std::endl; 01201 logfile << "Try loading a .cub file." << std::endl; 01202 // exit(0); don't abort, continue to run. MESHTOGEOM won't be computed rest should still work. 01203 } 01204 return iBase_SUCCESS; 01205 } 01206 01207 int CoreGen::load_geometries() 01208 // --------------------------------------------------------------------------- 01209 // Function: loads all the meshes and initializes copymesh and merge mesh objects 01210 // Input: none 01211 // Output: none 01212 // --------------------------------------------------------------------------- 01213 { 01214 logfile << "\n--Loading geometry files." << std::endl; 01215 01216 iGeom_getRootSet(igeom->instance(), &root_set, &err); 01217 ERRORR("Couldn't get the root set", err); 01218 01219 // create cm instances for each mesh file 01220 cg.resize(files.size()); 01221 01222 iBase_EntitySetHandle orig_set, temp_set, temp_set1; 01223 01224 iGeom_createEntSet(igeom->instance(), 0, &temp_set1, &err); 01225 ERRORR( "Problem creating entity set.",err ); 01226 01227 MEntVector vols_old; 01228 // loop reading all igeom->instance() files 01229 for (unsigned int i = 0; i < files.size(); i++) { 01230 iGeom_createEntSet(igeom->instance(), 0, &orig_set, &err); 01231 ERRORR( "Problem creating entity set.", err); 01232 01233 iGeom_createEntSet(igeom->instance(), 0, &temp_set, &err); 01234 ERRORR( "Problem creating entity set.",err ); 01235 01236 iGeom_load(igeom->instance(), files[i].c_str(), NULL, &err, 01237 strlen(files[i].c_str()), 0); 01238 ERRORR("Couldn't read geometry file.", err); 01239 mk_core()->populate_model_ents(0,-1,-1); 01240 01241 iBase_EntityHandle *entities = NULL; 01242 int entities_ehsize = 0, entities_ehallocated = 0; 01243 01244 iGeom_getEntities(igeom->instance(), root_set, iBase_REGION, &entities, 01245 &entities_ehsize, &entities_ehallocated, &err); 01246 ERRORR( "Problem getting entities." , err); 01247 01248 // add the entity 01249 for (int j = 0; j < entities_ehsize; j++) { 01250 iGeom_addEntToSet(igeom->instance(), entities[j], temp_set, &err); 01251 ERRORR( "Problem adding to set.", err ); 01252 } 01253 01254 iGeom_subtract(igeom->instance(), temp_set, temp_set1, &orig_set, &err); 01255 ERRORR( "Unable to subtract entity sets.", err ); 01256 01257 MEntVector vols; 01258 mk_core()->get_entities_by_dimension(3, vols); 01259 MEntVector vols1;// = vols - vols_old; 01260 01261 std::set_difference(vols.begin(), vols.end(), vols_old.begin(), vols_old.end(), std::inserter(vols1, vols1.begin())); 01262 01263 cg[i] = (CopyGeom*) mk_core()->construct_meshop("CopyGeom", vols1); 01264 cg[i]->set_name("copy_move_geom"); 01265 01266 assys.push_back(orig_set); 01267 01268 assys_index.push_back(i); 01269 01270 // store this set for subtraction with next entity set 01271 temp_set1 = temp_set; 01272 vols_old = vols; 01273 } 01274 logfile << "\n--Loaded geometry files.\n" << std::endl; 01275 01276 return iBase_SUCCESS; 01277 } 01278 01279 int CoreGen::shift_mn_ids(iBase_EntitySetHandle orig_set, int number) 01280 // --------------------------------------------------------------------------- 01281 // Function: loads all the meshes and initializes copymesh and merge mesh objects 01282 // Input: none 01283 // Output: none 01284 // --------------------------------------------------------------------------- 01285 { 01286 logfile << "Swapping MS and NS ids for " << number << std::endl; 01287 // get all the material sets in this assembly 01288 moab::Tag mattag, neutag; 01289 mb->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, mattag ); 01290 mb->tag_get_handle( "NEUMANN_SET", 1, MB_TYPE_INTEGER, neutag ); 01291 01292 int rval = 0; 01293 moab::Range matsets, neusets; 01294 01295 mb->get_entities_by_type_and_tag( (moab::EntityHandle)orig_set, MBENTITYSET, &mattag, 0, 1, matsets ); 01296 mb->get_entities_by_type_and_tag( (moab::EntityHandle)orig_set, MBENTITYSET, &neutag, 0, 1, neusets ); 01297 01298 int i = 0; 01299 moab::Range::iterator set_it; 01300 for (set_it = matsets.begin(); set_it != matsets.end(); set_it++) { 01301 moab::EntityHandle this_set = *set_it; 01302 01303 // get the id for this set 01304 int set_id; 01305 rval = mb->tag_get_data(mattag, &this_set, 1, &set_id); 01306 if(rval != moab::MB_SUCCESS) { 01307 std::cerr<<"getting tag data failed Code:"; 01308 std::string foo = ""; mb->get_last_error(foo); 01309 std::cerr<<"File Error: "<<foo<<std::endl; 01310 return 1; 01311 } 01312 01313 // set the new id for this set 01314 int new_id = all_ms_starts[number] + i; 01315 rval = mb->tag_set_data(mattag, &this_set, 1, &new_id); 01316 if(rval != moab::MB_SUCCESS) { 01317 std::cerr<<"getting tag data failed Code:"; 01318 std::string foo = ""; mb->get_last_error(foo); 01319 std::cerr<<"File Error: "<<foo<<std::endl; 01320 return 1; 01321 } 01322 ++i; 01323 } 01324 01325 int j = 0; 01326 for (set_it = neusets.begin(); set_it != neusets.end(); set_it++) { 01327 moab::EntityHandle this_set = *set_it; 01328 01329 // get the id for this set 01330 int set_id; 01331 rval = mb->tag_get_data(neutag, &this_set, 1, &set_id); 01332 if(rval != moab::MB_SUCCESS) { 01333 std::cerr<<"getting tag data failed Code:"; 01334 std::string foo = ""; mb->get_last_error(foo); 01335 std::cerr<<"File Error: "<<foo<<std::endl; 01336 return 1; 01337 } 01338 01339 // set the new id for this set 01340 int new_id = all_ns_starts[number] + j; 01341 rval = mb->tag_set_data(neutag, &this_set, 1, &new_id); 01342 if(rval != moab::MB_SUCCESS) { 01343 std::cerr<<"getting tag data failed Code:"; 01344 std::string foo = ""; mb->get_last_error(foo); 01345 std::cerr<<"File Error: "<<foo<<std::endl; 01346 return 1; 01347 } 01348 ++j; 01349 } 01350 return 0; 01351 } 01352 01353 int CoreGen::move_verts(iBase_EntitySetHandle set, const double *dx) 01354 // --------------------------------------------------------------------------- 01355 // Function: Change the coordinates for moving the assembly to its first loc. 01356 // Input: none 01357 // Output: none 01358 // --------------------------------------------------------------------------- 01359 { 01360 01361 int verts_ents_alloc = 0, verts_ents_size = 0; 01362 iBase_EntityHandle *verts_ents = NULL; 01363 01364 iMesh_getEntities(imesh->instance(), set, iBase_VERTEX, iMesh_ALL_TOPOLOGIES, 01365 &verts_ents, &verts_ents_alloc, &verts_ents_size, &err); 01366 ERRORR("Failed to get any entities from original set.", iBase_FAILURE); 01367 01368 double *coords = 0; 01369 int coords_alloc = 0, coords_size = 0; 01370 01371 iMesh_getVtxArrCoords(imesh->instance(), verts_ents, verts_ents_size, iBase_INTERLEAVED, 01372 &coords, &coords_alloc, &coords_size, &err); 01373 ERRORR("Failed to get vtx coords from set.", iBase_FAILURE); 01374 01375 for (int i = 0; i < verts_ents_size; i++) { 01376 coords[3 * i] += dx[0]; 01377 coords[3 * i + 1] += dx[1]; 01378 coords[3 * i + 2] += dx[2]; 01379 } 01380 01381 iMesh_setVtxArrCoords(imesh->instance(), verts_ents, verts_ents_size, iBase_INTERLEAVED, 01382 coords, coords_size, &err); 01383 ERRORR("Failed to set vtx coords.", iBase_FAILURE); 01384 01385 return iBase_SUCCESS; 01386 } 01387 01388 int CoreGen::move_geoms(iBase_EntitySetHandle set, const double *dx) 01389 // --------------------------------------------------------------------------- 01390 // Function: Change the coordinates for moving the assembly to its first loc. 01391 // Input: none 01392 // Output: none 01393 // --------------------------------------------------------------------------- 01394 { 01395 int entities_ehsize = 0, entities_ehallocated = 0; 01396 iBase_EntityHandle *entities = NULL; 01397 01398 iGeom_getEntities(igeom->instance(), set, iBase_ALL_TYPES, &entities, &entities_ehsize, 01399 &entities_ehallocated, &err); 01400 ERRORR("Failed to get entities from set.", iBase_FAILURE); 01401 01402 for (int i = 0; i < entities_ehsize; i++) { 01403 iGeom_moveEnt(igeom->instance(), entities[i], dx[0], dx[1], dx[2], &err); 01404 ERRORR("Failed to move geometries.", iBase_FAILURE); 01405 } 01406 return iBase_SUCCESS; 01407 } 01408 01409 int CoreGen::extrude() { 01410 // --------------------------------------------------------------------------- 01411 // Function: extrude 2D surface core 01412 // Input: none 01413 // Output: none 01414 // --------------------------------------------------------------------------- 01415 // extrude if this is a surface mesh 01416 if (set_DIM == 2 && extrude_flag == true) { // if surface geometry and extrude 01417 logfile << "Extruding surface mesh." << std::endl; 01418 //get entities for extrusion 01419 iBase_EntityHandle *ents = NULL; 01420 int ents_alloc = 0, ents_size; 01421 iMesh_getEntities(imesh->instance(), root_set, iBase_FACE, iMesh_ALL_TOPOLOGIES, 01422 &ents, &ents_alloc, &ents_size, &err); 01423 ERRORR("Trouble getting face mesh.", err); 01424 01425 // add entities for extrusion to a set 01426 iBase_EntitySetHandle set; 01427 iMesh_createEntSet(imesh->instance(), false, &set, &err); 01428 ERRORR("Trouble getting face mesh.", err); 01429 01430 iMesh_addEntArrToSet(imesh->instance(), ents, ents_size, set, &err); 01431 ERRORR("Trouble getting face mesh.", err); 01432 01433 ModelEnt me(mk_core(), iBase_EntitySetHandle(0), /*igeom instance*/0, 01434 (moab::EntityHandle)set); 01435 MEntVector selection; 01436 selection.push_back(&me); 01437 01438 // This tag needs to be set to the newly created extrude sets 01439 const char *tag_g1 = "GEOM_DIMENSION"; 01440 iBase_TagHandle gtag; 01441 iMesh_getTagHandle(imesh->instance(), tag_g1, >ag, &err, 14); 01442 ERRORR("Trouble getting geom dimension set.", err); 01443 01444 // This tag needs to be set to the newly created extrude sets 01445 const char *tag_m1 = "MATERIAL_SET"; 01446 iBase_TagHandle mtag; 01447 iMesh_getTagHandle(imesh->instance(), tag_m1, &mtag, &err, 12); 01448 ERRORR("Trouble getting material set.", err); 01449 01450 // This tag needs to be set to the newly created extrude sets 01451 const char *tag_n1 = "NEUMANN_SET"; 01452 iBase_TagHandle ntag; 01453 iMesh_getTagHandle(imesh->instance(), tag_n1, &ntag, &err, 11); 01454 ERRORR("Trouble getting neumann set.", err); 01455 01456 01457 double v[] = { 0, 0, z_height }; 01458 int steps = z_divisions; 01459 em = (ExtrudeMesh*) mk_core()->construct_meshop("ExtrudeMesh", selection); 01460 em->set_transform(Extrude::Translate(v, steps)); 01461 em->copy_faces(true); 01462 01463 01464 SimpleArray<iBase_EntitySetHandle> msets; 01465 iMesh_getEntSetsByTagsRec(imesh->instance(), root_set, &mtag, NULL, 1, 0, 01466 ARRAY_INOUT(msets), &err); 01467 ERRORR("Trouble getting entity set.", err); 01468 01469 for (int i = 0; i < msets.size(); i++) { 01470 em->extrude_sets().add_set((iMesh::EntitySetHandle)msets[i]); 01471 } 01472 01473 // some entity tag types are always copy or expand 01474 em->extrude_sets().add_tag("MATERIAL_SET"); 01475 01476 // run 01477 em->setup_this(); 01478 em->execute_this(); 01479 01480 iMesh_destroyEntSet(imesh->instance(), set, &err); 01481 ERRORR("Error in destroying ent set of faces after extrusion is done.", err); 01482 01483 msets.clear(); 01484 iMesh_getEntSetsByTagsRec(imesh->instance(), root_set, &mtag, NULL, 1, 0, 01485 ARRAY_INOUT(msets), &err); 01486 ERRORR("Trouble getting entity set.", err); 01487 // now delete all the 2D material sets: 01488 for (int i = 0; i < msets.size(); i++) { 01489 int num =0; 01490 iMesh_getNumOfType(imesh->instance(), msets[i], iBase_REGION, &num, &err); 01491 ERRORR("Trouble getting num entities.", err); 01492 if(num == 0){ // material sets with quads 01493 iMesh_destroyEntSet(imesh->instance(), msets[i], &err); 01494 ERRORR("Trouble destroying set.", err); 01495 } 01496 } 01497 01498 01499 // Step 2: get all max. value of neumann sets, then, for newly created NS set a new value and GD =2 tag. 01500 01501 iBase_EntityHandle *ents1 = NULL; 01502 int ents_alloc1 = 0, ents_size1 = 0; 01503 01504 SimpleArray<iBase_EntitySetHandle> nsets; 01505 iMesh_getEntSetsByTagsRec(imesh->instance(), root_set, &ntag, NULL, 01506 1, 0, ARRAY_INOUT(nsets), &err); 01507 ERRORR("Trouble getting entity set.", err); 01508 01509 int max_nset_value = 0; 01510 for (int i = 0; i < nsets.size(); i++) { 01511 int nvalue; 01512 iMesh_getEntSetIntData(imesh->instance(), nsets[i], ntag, &nvalue, &err); 01513 ERRORR("Trouble getting entity set.", err); 01514 if (nvalue > max_nset_value) 01515 max_nset_value = nvalue; 01516 } 01517 01518 for (int i = 0; i < nsets.size(); i++) { 01519 01520 iMesh_getEntities(imesh->instance(), nsets[i], 01521 iBase_FACE, iMesh_ALL_TOPOLOGIES, 01522 &ents1, &ents_alloc1, &ents_size1, &err); 01523 ERRORR("Trouble getting face mesh.", err); 01524 01525 if(ents_size1 > 0) { 01526 // set GEOM_DIMENSION tag = 2 and renumber the neumann set 01527 const int gd = 2; 01528 const int nvalue = max_nset_value + i; 01529 01530 iMesh_setEntSetIntData(imesh->instance(), nsets[i], ntag, nvalue, &err); 01531 ERRORR("Trouble getting entity set.", err); 01532 01533 iMesh_setEntSetIntData(imesh->instance(),nsets[i], gtag, gd, &err); 01534 ERRORR("Trouble getting entity set.", err); 01535 } 01536 ents_alloc1 = 0; 01537 ents_size1 = 0; 01538 *ents1 = NULL; 01539 } 01540 } 01541 return iBase_SUCCESS; 01542 } 01543 01544 int CoreGen::refine_coremodel(){ 01545 01546 #ifdef USE_MPI 01547 std::cout<<"REFINEMENT BLOCK"<<std::endl; 01548 /*********************************************/ 01549 // refine the mesh 01550 /*********************************************/ 01551 int num_levels = deg.size(); 01552 moab::NestedRefine uref(dynamic_cast<Core*>(mb), pc, 0); 01553 std::vector<moab::EntityHandle> lsets; 01554 01555 std::cout<<"Starting hierarchy generation"<<std::endl; 01556 01557 uref.generate_mesh_hierarchy(num_levels, °[0], lsets); 01558 uref.update_special_tags(num_levels, lsets[num_levels]); 01559 std::cout<<"Finished hierarchy generation in "<<uref.timeall.tm_total<<" secs"<<std::endl; 01560 01561 std::cout<<"Time taken for refinement "<<uref.timeall.tm_refine<<" secs"<<std::endl; 01562 std::cout<<"Time taken for resolving shared interface "<<uref.timeall.tm_resolve<<" secs"<<std::endl; 01563 01564 moab::Range all_ents, ents[4]; 01565 for (int l=0; l<num_levels; l++) 01566 { 01567 all_ents.clear(); 01568 ents[0].clear(); ents[1].clear(); ents[2].clear(); ents[3].clear(); 01569 mb->get_entities_by_handle(lsets[l+1], all_ents); //MB_CHK_ERR(error); 01570 01571 for (int k=0; k<4; k++) 01572 ents[k] = all_ents.subset_by_dimension(k); 01573 01574 std::cout<<"Mesh size for level "<<l+1<<" :: nverts = "<<ents[0].size()<<", nedges = "<<ents[1].size()<<", nfaces = "<<ents[2].size()<<", ncells = "<<ents[3].size()<<std::endl; 01575 } 01576 #endif 01577 return 0; 01578 } 01579 01580 int CoreGen::create_neumannset() { 01581 // --------------------------------------------------------------------------- 01582 // Function: create Neumann set on the whole core 01583 // Input: none 01584 // Output: none 01585 // --------------------------------------------------------------------------- 01586 if (nss_flag == true || nsb_flag == true || nst_flag == true || nssall_flag == true) { 01587 logfile << "Creating NeumannSet." << std::endl; 01588 01589 if (extrude_flag == true) 01590 set_DIM = 3; 01591 01592 int err = 0, z_flag = 0, i, ents_alloc = 0, ents_size; 01593 double z1 = 0.0; 01594 iBase_TagHandle ntag1, gtag1, nametag1; 01595 iBase_EntityHandle *ents = NULL; 01596 iBase_EntitySetHandle set = NULL, set_z1 = NULL, set_z2 = NULL; 01597 std::vector<iBase_EntitySetHandle> set_side; 01598 01599 //get entities for skinner 01600 if(set_DIM ==2) { // if surface geometry specified 01601 iMesh_getEntities(imesh->instance(), root_set, 01602 iBase_FACE, iMesh_ALL_TOPOLOGIES, 01603 &ents, &ents_alloc, &ents_size, &err); 01604 } 01605 else { 01606 iMesh_getEntities(imesh->instance(), root_set, 01607 iBase_REGION, iMesh_ALL_TOPOLOGIES, 01608 &ents, &ents_alloc, &ents_size, &err); 01609 } 01610 ERRORR("Trouble getting entities for specifying neumannsets via skinner.", err); 01611 01612 // assign a name to the tag 01613 const char *ch_name1 = "NAME"; 01614 // get tag handle 01615 const char *tag_neumann1 = "NEUMANN_SET"; 01616 const char *global_id1 = "GLOBAL_ID"; 01617 01618 iMesh_getTagHandle(imesh->instance(), tag_neumann1, &ntag1, &err, 12); 01619 ERRORR("Trouble getting handle.", err); 01620 01621 iMesh_getTagHandle(imesh->instance(), global_id1, >ag1, &err, 9); 01622 ERRORR("Trouble getting handle.", err); 01623 01624 iMesh_getTagHandle(imesh->instance(), ch_name1, &nametag1, &err, 4); 01625 ERRORR("Trouble getting handle.", err); 01626 01627 iMesh_createEntSet(imesh->instance(),0, &set, &err); // for all other sides 01628 ERRORR("Trouble creating set handle.", err); 01629 01630 if (set_DIM == 3) { // sets for collecting top and bottom surface 01631 iMesh_createEntSet(imesh->instance(),0, &set_z1, &err); 01632 ERRORR("Trouble creating set handle.", err); 01633 01634 iMesh_createEntSet(imesh->instance(),0, &set_z2, &err); 01635 ERRORR("Trouble creating set handle.", err); 01636 01637 set_side.resize(num_nsside); 01638 for(int i=0; i<num_nsside; i++){ 01639 iMesh_createEntSet(imesh->instance(),0, &set_side[i], &err); 01640 ERRORR("Trouble creating set handle.", err); 01641 } 01642 } 01643 01644 moab::Range tmp_elems; 01645 tmp_elems.insert((EntityHandle*)ents, (moab::EntityHandle*)ents + ents_size); 01646 01647 // get the skin of the entities 01648 moab::Skinner skinner(mb); 01649 moab::Range skin_range; 01650 moab::ErrorCode result; 01651 moab::Range::iterator rit; 01652 01653 result = skinner.find_skin(0, tmp_elems, set_DIM-1, skin_range); 01654 if (MB_SUCCESS != result) return result; 01655 01656 for (rit = skin_range.begin(), i = 0; rit != skin_range.end(); rit++, i++) { 01657 if(set_DIM == 3) { // filter top and bottom 01658 int num_vertex=0, size_vertex =0; 01659 iBase_EntityHandle *vertex = NULL; 01660 iMesh_getEntAdj(imesh->instance(), (iBase_EntityHandle)(*rit), iBase_VERTEX, &vertex, 01661 &num_vertex, &size_vertex, &err); 01662 ERRORR("Trouble getting number of entities after merge.", err); 01663 01664 double *coords = NULL; 01665 int coords_alloc = 0, coords_size=0; 01666 iMesh_getVtxArrCoords(imesh->instance(), vertex, size_vertex, iBase_INTERLEAVED, 01667 &coords, &coords_alloc, &coords_size, &err); 01668 ERRORR("Trouble getting number of entities after merge.", err); 01669 double ztemp = coords[2]; 01670 int flag = 0; 01671 for (int p=1; p<num_vertex; p++) { 01672 double z1 = coords[3*p+2]; 01673 if( fabs(ztemp-z1) >= merge_tol) { 01674 flag = 1; 01675 continue; 01676 } 01677 } 01678 if(flag == 0) { // this is top or bottom surface 01679 if (z_flag == 0) { // store z-coord this is the first call 01680 z_flag = 1; 01681 z1 = ztemp; 01682 } 01683 if( fabs(ztemp-z1) <= merge_tol) { 01684 iMesh_addEntToSet(imesh->instance(), (iBase_EntityHandle)(*rit), set_z1, &err); 01685 ERRORR("Trouble getting number of entities after merge.", err); 01686 } 01687 else { 01688 iMesh_addEntToSet(imesh->instance(), (iBase_EntityHandle)(*rit), set_z2, &err); 01689 ERRORR("Trouble getting number of entities after merge.", err); 01690 } 01691 } 01692 else if (flag == 1) { // add the faces that are not top or bottom surface 01693 // filter the sidesets based on their x and y coords 01694 01695 for(int k=0; k<num_nsside; k++){ 01696 if ( fabs((coords[0])*nsx[k] + (coords[1])*nsy[k] + nsc[k]) <= merge_tol 01697 && fabs((coords[3])*nsx[k] + (coords[4])*nsy[k] + nsc[k]) <= merge_tol 01698 && fabs((coords[6])*nsx[k] + (coords[7])*nsy[k] + nsc[k]) <= merge_tol) { 01699 iMesh_addEntToSet(imesh->instance(), (iBase_EntityHandle)(*rit), (iBase_EntitySetHandle) set_side[k], &err); 01700 ERRORR("Trouble getting number of entities after merge.", err); 01701 continue; 01702 } 01703 else{ // outside the specified 01704 iMesh_addEntToSet(imesh->instance(), (iBase_EntityHandle)(*rit), set, &err); 01705 ERRORR("Trouble getting number of entities after merge.", err); 01706 continue; 01707 } 01708 } 01709 if(num_nsside == 0){ 01710 iMesh_addEntToSet(imesh->instance(), (iBase_EntityHandle)(*rit), set, &err); 01711 ERRORR("Trouble getting number of entities after merge.", err); 01712 } 01713 01714 } 01715 } 01716 else if(set_DIM == 2) { // edges add all for sideset 01717 iMesh_addEntToSet(imesh->instance(), (iBase_EntityHandle)(*rit), set, &err); 01718 ERRORR("Trouble getting number of entities after merge.", err); 01719 } 01720 } 01721 01722 if (set_DIM == 3) { 01723 if (nst_flag == true || nsb_flag == true) { 01724 01725 iMesh_setEntSetIntData( imesh->instance(), set_z1, ntag1, nst_Id, &err); 01726 ERRORR("Trouble getting handle.", err); 01727 01728 iMesh_setEntSetIntData( imesh->instance(), set_z1, gtag1, nst_Id, &err); 01729 ERRORR("Trouble getting handle.", err); 01730 01731 std::string name1 = "core_top_ss"; 01732 iMesh_setEntSetData( imesh->instance(), set_z1, nametag1, name1.c_str(), 11, &err); 01733 ERRORR("Trouble getting handle.", err); 01734 01735 iMesh_setEntSetIntData( imesh->instance(), set_z2, ntag1, nsb_Id, &err); 01736 ERRORR("Trouble getting handle.", err); 01737 01738 iMesh_setEntSetIntData( imesh->instance(), set_z2, gtag1, nsb_Id, &err); 01739 ERRORR("Trouble getting handle.", err); 01740 01741 std::string name2 = "core_bottom_ss"; 01742 iMesh_setEntSetData( imesh->instance(), set_z2, nametag1, name2.c_str(), 14, &err); 01743 ERRORR("Trouble getting handle.", err); 01744 01745 for(int j=0; j<num_nsside; j++){ 01746 iMesh_setEntSetIntData( imesh->instance(), set_side[j], ntag1, nss_Id[j], &err); 01747 ERRORR("Trouble getting handle.", err); 01748 01749 iMesh_setEntSetIntData( imesh->instance(), set_side[j], gtag1, nss_Id[j], &err); 01750 ERRORR("Trouble getting handle.", err); 01751 01752 std::stringstream ss; 01753 ss << j; 01754 std::string name3 = "side" + ss.str(); 01755 iMesh_setEntSetData( imesh->instance(), set_side[j], nametag1, name3.c_str(), 8, &err); 01756 ERRORR("Trouble getting handle.", err); 01757 } 01758 } 01759 } 01760 // same for both 2D and 3D models 01761 if (nssall_flag == true) { 01762 iMesh_setEntSetIntData( imesh->instance(), set, ntag1, nssall_Id, &err); 01763 ERRORR("Trouble getting handle.", err); 01764 01765 iMesh_setEntSetIntData( imesh->instance(), set, gtag1, nssall_Id, &err); 01766 ERRORR("Trouble getting handle.", err); 01767 01768 std::string name = "all_sides"; 01769 iMesh_setEntSetData( imesh->instance(), set, nametag1, name.c_str(), 9, &err); 01770 ERRORR("Trouble getting handle.", err); 01771 } 01772 } 01773 return iBase_SUCCESS; 01774 } 01775 01776 void CoreGen::IOErrorHandler (ErrorStates ECode) const 01777 // --------------------------------------------------------------------------- 01778 // Function: displays error messages related to input file parsing data 01779 // Input: error code 01780 // Output: none 01781 // --------------------------------------------------------------------------- 01782 { 01783 std::cerr << '\n'; 01784 if (ECode == INVALIDINPUT) // invalid input 01785 std::cerr << "Invalid input."; 01786 else if (ECode == ENEGATIVE) // invalid input 01787 std::cerr << "Unexpected negative value."; 01788 else 01789 std::cerr << "Unknown error ...?"; 01790 01791 std::cerr << '\n' << "Error reading input file, line : " << linenumber; 01792 std::cerr << std::endl; 01793 exit (1); 01794 } 01795 01796 int CoreGen::write_minfofile() 01797 // --------------------------------------------------------------------------- 01798 // Function: write the spreadsheet mesh info file based on inputs read from mesh & input file 01799 // Input: none 01800 // Output: none 01801 // --------------------------------------------------------------------------- 01802 { 01803 logfile << "Writing mesh info file indicating elements and pin number" << std::endl; 01804 01805 moab::Tag ntag; 01806 mb->tag_get_handle(NAME_TAG_NAME, NAME_TAG_SIZE, MB_TYPE_OPAQUE, ntag); 01807 moab::Tag mattag; 01808 mb->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, mattag ); 01809 int rval = 0; 01810 moab::Range matsets; 01811 std::vector <EntityHandle> set_ents; 01812 01813 mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &mattag, 0, 1, matsets ); 01814 01815 moab::Range::iterator set_it; 01816 for (set_it = matsets.begin(); set_it != matsets.end(); set_it++) { 01817 moab::EntityHandle this_set = *set_it; 01818 01819 // get the id for this set 01820 int set_id; 01821 rval = mb->tag_get_data(mattag, &this_set, 1, &set_id); 01822 if(rval != moab::MB_SUCCESS) { 01823 std::cerr<<"getting tag data failed Code:"; 01824 std::string foo = ""; mb->get_last_error(foo); 01825 std::cerr<<"File Error: "<<foo<<std::endl; 01826 return 1; 01827 } 01828 char name[NAME_TAG_SIZE]; 01829 rval = mb->tag_get_data(ntag, &this_set, 1, &name); 01830 if(rval != moab::MB_SUCCESS) { 01831 std::cerr<<"getting tag data failed Code:"; 01832 std::string foo = ""; mb->get_last_error(foo); 01833 std::cerr<<"File Error: "<<foo<<std::endl; 01834 return 1; 01835 } 01836 // check for the special block _xp created in AssyGen stage 01837 // now print out elements for each pin on the mesh info file 01838 if(name[0]=='_' && name[1]=='x' && name[2] == 'p'){ 01839 01840 // get the entities in the set, recursively 01841 rval = mb->get_entities_by_dimension(this_set, 3, set_ents, true); 01842 01843 logfile << "Block: " << set_id << " has " 01844 << set_ents.size() << " entities. Name = " << name; 01845 01846 // loop thro' all the elements in all the sets 01847 for (int i=0;i<int(set_ents.size());i++){ 01848 01849 std::vector<EntityHandle> conn; 01850 EntityHandle handle = set_ents[i]; 01851 01852 // get the connectivity of this element 01853 rval = mb->get_connectivity(&handle, 1, conn); 01854 double coords[3]; 01855 double x_sum = 0.0, y_sum = 0.0, z_sum = 0.0; 01856 for (int j = 0; j<int(conn.size()); ++j){ 01857 rval = mb->get_coords(&conn[j], 1, coords); 01858 x_sum+=coords[0]; 01859 y_sum+=coords[1]; 01860 z_sum+=coords[2]; 01861 } 01862 int p = 3; 01863 while(name[p]!='\0'){ 01864 minfo_file << name[p]; 01865 ++p; 01866 } 01867 minfo_file << " \t" << x_sum/conn.size() << " \t" << y_sum/conn.size() << " \t" << z_sum/conn.size() << std::endl; 01868 } 01869 logfile << ". Deleting block " << set_id << std::endl; 01870 mb->delete_entities(&this_set, 1); 01871 set_ents.clear(); 01872 } 01873 01874 } 01875 return 0; 01876 } 01877 01878 01879 int CoreGen::prepareIO(int argc, char *argv[], int nrank, int nprocs, std::string TestDir) 01880 // ----------------------------------------------------------------------------------- 01881 // Function: Obtains file names and opens input/output files and then read/write them 01882 // Input: command line arguments 01883 // Output: none 01884 // ----------------------------------------------------------------------------------- 01885 { 01886 // set rank and total number of processors 01887 rank = nrank; 01888 procs = nprocs; 01889 sTime = clock(); 01890 testdir = TestDir; 01891 if (argc > 1) { 01892 if (argv[1][0] == '-' && argv[1][1] == 'm') { 01893 // set to zero, when run_flag = 1, program runs and does copy, move, merge, extrude, assign gids, save and close 01894 run_flag = 0; 01895 } 01896 } 01897 01898 bool bDone = false; 01899 do { 01900 if (2 == argc) { 01901 if (argv[1][0] == '-' && nrank == 0) { 01902 if (argv[1][1] == 'h') { 01903 logfile << "Usage: coregen [-t -m -h] <coregen input file>" 01904 << std::endl; 01905 logfile << " -t print timing and memory usage info in each step" 01906 << std::endl; 01907 logfile << " -m create makefile only" << std::endl; 01908 logfile << " -h print help" << std::endl; 01909 logfile << "\nInstruction on writing coregen input file can be found at: " 01910 << std::endl; 01911 logfile << " http://press3.mcs.anl.gov/sigma/meshkit/rgg/coregen-input-file-keyword-definitions/" 01912 << std::endl; 01913 exit(0); 01914 } 01915 } 01916 01917 iname = argv[1]; 01918 ifile = iname + ".inp"; 01919 outfile = iname + ".h5m"; 01920 mfile = iname + ".makefile"; 01921 meshtogeomfile = "meshtogeom." + iname; 01922 infofile = iname + "_info.csv"; 01923 minfofile = iname + "_mesh_info.csv"; 01924 logfilename = iname + ".log"; 01925 } else if (3 == argc) { 01926 int i = 1;// will loop through arguments, and process them 01927 for (i = 1; i < argc - 1; i++) { 01928 if (argv[i][0] == '-') { 01929 switch (argv[i][1]) { 01930 case 'm': { 01931 if (nrank == 0) { 01932 logfile << "Creating Make/Info file Only" << std::endl; 01933 } 01934 // only makefile creation specified 01935 iname = argv[2]; 01936 ifile = iname + ".inp"; 01937 outfile = iname + ".h5m"; 01938 mfile = iname + ".makefile"; 01939 meshtogeomfile = "meshtogeom." + iname; 01940 infofile = iname + "_info.csv"; 01941 minfofile = iname + "_mesh_info.csv"; 01942 logfilename = iname + ".log"; 01943 break; 01944 } 01945 case 't': { 01946 mem_tflag = true; 01947 iname = argv[2]; 01948 ifile = iname + ".inp"; 01949 outfile = iname + ".h5m"; 01950 mfile = iname + ".makefile"; 01951 meshtogeomfile = "meshtogeom." + iname; 01952 infofile = iname + "_info.csv"; 01953 minfofile = iname + "_mesh_info.csv"; 01954 logfilename = iname + ".log"; 01955 break; 01956 } 01957 case 'h': { 01958 if (nrank == 0) { 01959 logfile << "Usage: coregen [-t -m -h] <coregen input file>" 01960 << std::endl; 01961 logfile << " -t print timing and memory usage info in each step" 01962 << std::endl; 01963 logfile << " -m create makefile only" 01964 << std::endl; 01965 logfile << " -h print help" << std::endl; 01966 logfile << "\nInstruction on writing coregen input file can also be found at: " 01967 << std::endl; 01968 logfile << " http://press3.mcs.anl.gov/sigma/meshkit/rgg/coregen-input-file-keyword-definitions/" 01969 << std::endl; 01970 exit(0); 01971 break; 01972 } 01973 } 01974 } 01975 } 01976 } 01977 } else { //default case 01978 if (nrank == 0) { 01979 logfile << "Usage: " << argv[0] 01980 << " <input file> WITHOUT EXTENSION" << std::endl; 01981 } 01982 iname = TestDir + "/" + COREGEN_DEFAULT_TEST_FILE; 01983 ifile = iname + ".inp"; 01984 std::string temp = CTEST_FILE_NAME; 01985 outfile = temp + ".h5m"; 01986 mfile = temp + ".makefile"; 01987 meshtogeomfile = "meshtogeom." + temp; 01988 infofile = temp + "_info.csv"; 01989 minfofile = temp + "_mesh_info.csv"; 01990 logfilename = temp + ".log"; 01991 } 01992 01993 // open the file 01994 file_input.open(ifile.c_str(), std::ios::in); 01995 logfile.coss.open(logfilename.c_str(), std::ios::out); 01996 01997 /*********************************************/ 01998 // Print banner on logfile and std out 01999 /*********************************************/ 02000 if (rank == 0) { 02001 banner(); 02002 Timer.GetDateTime(szDateTime); 02003 logfile << "\nStarting out at : " << szDateTime << "\n"; 02004 } 02005 02006 if (!file_input) { 02007 if (nrank == 0) { 02008 logfile << "Default case input file located here: <MeshKit/data>" << std::endl; 02009 logfile << "Usage: coregen [-t -m -h] <coregen input file>" 02010 << std::endl; 02011 logfile << " -t print timing and memory usage info in each step" 02012 << std::endl; 02013 logfile << " -m create makefile only" << std::endl; 02014 logfile << " -h print help" << std::endl; 02015 logfile << "\nInstruction on writing coregen input file can be found at: " 02016 << std::endl; 02017 logfile << " http://press3.mcs.anl.gov/sigma/meshkit/rgg/coregen-input-file-keyword-definitions/" 02018 << std::endl; 02019 logfile << "ERROR - Invalid INPUT FILE specified. Hint: input file name must be specified without extension" << std::endl; 02020 02021 } 02022 file_input.clear(); 02023 exit(1); 02024 } else 02025 bDone = true; // file opened successfully 02026 } while (!bDone); 02027 02028 // open Makefile-rgg 02029 do { 02030 make_file.open(mfile.c_str(), std::ios::out); 02031 if (!make_file) { 02032 if (nrank == 0) { 02033 logfile << "Unable to open makefile for writing" << std::endl; 02034 } 02035 make_file.clear(); 02036 } else 02037 bDone = true; // file opened successfully 02038 } while (!bDone); 02039 02040 if (nrank == 0) { 02041 logfile << "\nEntered input file name: " << ifile << std::endl; 02042 } 02043 02044 // now call the functions to read and write 02045 err = read_inputs_phase1(argc, argv); 02046 ERRORR("Failed to read inputs in phase1.", 1); 02047 02048 err = read_inputs_phase2(argc, argv); 02049 ERRORR("Failed to read inputs in phase2.", 1); 02050 02051 // open meshtogeomfile 02052 if(compute_meshtogeom){ 02053 meshtogeom_file.coss.open(meshtogeomfile.c_str(), std::ios::out); 02054 logfile << "Created meshtogeom file: " << meshtogeomfile << std::endl; 02055 } 02056 // open info file 02057 if(strcmp(info.c_str(),"on") == 0 && nrank == 0){ 02058 do { 02059 info_file.open(infofile.c_str(), std::ios::out); 02060 if (!info_file) { 02061 if (nrank == 0) { 02062 logfile << "Unable to open makefile for writing" << std::endl; 02063 } 02064 info_file.clear(); 02065 } else 02066 bDone = true; // file opened successfully 02067 logfile << "Created core info file: " << infofile << std::endl; 02068 } while (!bDone); 02069 02070 info_file << "assm index" << " \t" << "assm number" << " \t" << "dX" << " \t" << "dY" << " \t" << "dZ" << " \t" << "rank" << std::endl; 02071 } 02072 02073 // open mesh info file 02074 if(strcmp(minfo.c_str(),"on") == 0 && nrank == 0){ 02075 do { 02076 minfo_file.open(minfofile.c_str(), std::ios::out); 02077 if (!info_file) { 02078 if (nrank == 0) { 02079 logfile << "Unable to open minfofile for writing" << std::endl; 02080 } 02081 minfo_file.clear(); 02082 } else 02083 bDone = true; // file opened successfully 02084 logfile << "Created mesh details info file: " << minfofile << std::endl; 02085 } while (!bDone); 02086 minfo_file << "pin_number" << " \t" << "x_centroid" << " \t" << "y_centroid" << " \t" << "z_centroid" << std::endl; 02087 } 02088 if (nrank == 0) { 02089 err = write_makefile(); 02090 ERRORR("Failed to write a makefile.", 1); 02091 } 02092 return 0; 02093 } 02094 02095 02096 int CoreGen::read_inputs_phase1(int argc, char *argv[]) { 02097 // --------------------------------------------------------------------------- 02098 // Function: Reads the dimension and symmetry of the problem 02099 // Input: none 02100 // Output: none 02101 // --------------------------------------------------------------------------- 02102 CParser parse; 02103 for (;;) { 02104 if (!parse.ReadNextLine(file_input, linenumber, input_string, MAXCHARS, 02105 comment)) 02106 ERRORR("Reading input file failed",1); 02107 // logfile << input_string << std::endl; 02108 if (input_string.substr(0, 11) == "problemtype") { 02109 std::istringstream formatString(input_string); 02110 formatString >> card >> prob_type; 02111 if(((strcmp (prob_type.c_str(), "geometry") != 0) 02112 && (strcmp (prob_type.c_str(), "mesh") != 0)) || formatString.fail()) 02113 IOErrorHandler (INVALIDINPUT); 02114 if ((strcmp(prob_type.c_str(), "geometry") == 0)) { 02115 prob_type = "geometry"; 02116 } 02117 } 02118 if (input_string.substr(0, 8) == "geometry" && input_string.substr(0, 12) != "geometrytype") { 02119 std::istringstream formatString(input_string); 02120 formatString >> card >> geometry; 02121 if(((strcmp (geometry.c_str(), "volume") != 0) 02122 && (strcmp (geometry.c_str(), "surface") != 0)) || formatString.fail()) 02123 IOErrorHandler (INVALIDINPUT); 02124 if ((strcmp(geometry.c_str(), "surface") == 0)) { 02125 set_DIM = 2; 02126 } 02127 } 02128 02129 // igeom->instance() engine 02130 if (input_string.substr(0, 10) == "geomengine") { 02131 std::istringstream formatString(input_string); 02132 formatString >> card >> geom_engine; 02133 if(((strcmp (geom_engine.c_str(), "acis") != 0) 02134 && (strcmp (geom_engine.c_str(), "occ") != 0)) || formatString.fail()) 02135 IOErrorHandler (INVALIDINPUT); 02136 } 02137 02138 // symmetry 02139 if (input_string.substr(0, 8) == "symmetry") { 02140 std::istringstream formatString(input_string); 02141 formatString >> card >> symm; 02142 if((symm !=1 && symm !=6 && symm !=12) || formatString.fail()) 02143 IOErrorHandler (INVALIDINPUT); 02144 } 02145 // element type 02146 if (input_string.substr(0, 11) == "elementtype") { 02147 std::istringstream formatString(input_string); 02148 formatString >> card >> etype; 02149 if(etype == "hex27") 02150 have_hex27 = true; 02151 else if (etype == "hex8") 02152 have_hex27 = false; 02153 else 02154 IOErrorHandler (INVALIDINPUT); 02155 } 02156 // UMR parameters 02157 if (input_string.substr(0, 3) == "umr") { 02158 umr_flag = true; 02159 int temp_deg = 0; 02160 std::istringstream formatString(input_string); 02161 formatString >> card >> nDegree; 02162 for(int i=1; i<=nDegree;i++){ 02163 formatString >> temp_deg; 02164 deg.push_back(temp_deg); 02165 } 02166 if(formatString.fail()) 02167 IOErrorHandler (INVALIDINPUT); 02168 } 02169 // merge tolerance 02170 if (input_string.substr(0, 14) == "mergetolerance") { 02171 std::istringstream formatString(input_string); 02172 formatString >> card >> merge_tol; 02173 if(merge_tol < 0 || formatString.fail()) 02174 IOErrorHandler (ENEGATIVE); 02175 } 02176 02177 // save onefile for each proc (multiple) flag 02178 if (input_string.substr(0, 12) == "saveparallel") { 02179 std::istringstream formatString(input_string); 02180 formatString >> card >> savefiles; 02181 if(formatString.fail()) 02182 IOErrorHandler (INVALIDINPUT); 02183 } 02184 // info flag 02185 if (input_string.substr(0, 4) == "info") { 02186 std::istringstream formatString(input_string); 02187 formatString >> card >> info; 02188 if(formatString.fail()) 02189 IOErrorHandler (INVALIDINPUT); 02190 } 02191 // info flag 02192 if (input_string.substr(0, 8) == "meshinfo") { 02193 std::istringstream formatString(input_string); 02194 formatString >> card >> minfo; 02195 if(formatString.fail()) 02196 IOErrorHandler (INVALIDINPUT); 02197 } 02198 // neumannset card 02199 if (input_string.substr(0, 10) == "neumannset") { 02200 std::istringstream formatString(input_string); 02201 std::string nsLoc = "", temp1, temp2, temp3; 02202 double x, y, c; 02203 int nsId = 0; 02204 formatString >> card >> nsLoc >> nsId; 02205 if(nsId < 0 || formatString.fail()) 02206 IOErrorHandler (INVALIDINPUT); 02207 if ((strcmp(nsLoc.c_str(), "top") == 0)) { 02208 nst_flag = true; 02209 nst_Id = nsId; 02210 } else if ((strcmp(nsLoc.c_str(), "bot") == 0)) { 02211 nsb_flag = true; 02212 nsb_Id = nsId; 02213 } else if ((strcmp(nsLoc.c_str(), "sideall") == 0)) { 02214 nssall_flag = true; 02215 nssall_Id = nsId; 02216 } else if ((strcmp(nsLoc.c_str(), "side") == 0)) { 02217 nss_Id.push_back(nsId); 02218 // we are reading the equation of a straight line ax + by + c = 0 02219 formatString >> temp1 >> x >> temp2 >> y >> temp3 >> c; 02220 if(formatString.fail()) 02221 IOErrorHandler (INVALIDINPUT); 02222 nsx.push_back(x); 02223 nsy.push_back(y); 02224 nsc.push_back(c); 02225 02226 ++num_nsside; 02227 nss_flag = true; 02228 } else { 02229 logfile << "Invalid Neumann set specification" << std::endl; 02230 } 02231 } 02232 02233 // breaking condition 02234 if (input_string.substr(0, 3) == "end") { 02235 std::istringstream formatstring(input_string); 02236 break; 02237 } 02238 } 02239 return iBase_SUCCESS; 02240 } 02241 02242 int CoreGen::read_inputs_phase2(int argc, char *argv[]) 02243 // --------------------------------------------------------------------------- 02244 // Function: read all the inputs 02245 // Input: command line arguments 02246 // Output: none 02247 // --------------------------------------------------------------------------- 02248 { 02249 //Rewind the input file 02250 file_input.clear(std::ios_base::goodbit); 02251 file_input.seekg(0L, std::ios::beg); 02252 linenumber = 0; 02253 02254 CParser parse; 02255 for (;;) { 02256 if (!parse.ReadNextLine(file_input, linenumber, input_string, MAXCHARS, 02257 comment)) 02258 ERRORR("Reading input file failed",1); 02259 02260 if (input_string.substr(0, 12) == "geometrytype" ) { 02261 std::istringstream formatString(input_string); 02262 formatString >> card >> geom_type; 02263 if(formatString.fail()) 02264 IOErrorHandler (INVALIDINPUT); 02265 02266 if(geom_type.substr(0,3) == "hex"){ // for all hex type assemblies read the pitch and mesh files names 02267 bool reading_assemblies = false; 02268 do{ 02269 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02270 MAXCHARS, comment)) 02271 ERRORR("Reading input file failed",1); 02272 std::istringstream formatString(input_string); 02273 if (input_string.substr(0, 10) == "assemblies"){ 02274 formatString >> card >> nassys >> pitch; 02275 if(nassys < 0 || formatString.fail()) 02276 IOErrorHandler (INVALIDINPUT); 02277 02278 // reading file and alias names 02279 if(!parse_assembly_names(parse, argc, argv)) 02280 ERRORR("error parsing names of assemblies",1); 02281 reading_assemblies = true; 02282 } 02283 } while (reading_assemblies == false) ; 02284 } 02285 // read lattice info for all assemblies, also assemblies for rect assemblies. 02286 if (geom_type == "hexvertex" && symm == 6) { 02287 02288 // reading lattice 02289 bool reading_lattice = false; 02290 do{ 02291 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02292 MAXCHARS, comment)) 02293 ERRORR("Reading input file failed",1); 02294 std::istringstream formatString1(input_string); 02295 if (input_string.substr(0,7) == "lattice"){ 02296 reading_lattice = true; 02297 formatString1 >> card >> nrings; 02298 if(nrings < 0 || formatString1.fail()) 02299 IOErrorHandler (INVALIDINPUT); 02300 if (nrings % 2 == 0) 02301 tot_assys = (nrings * (nrings)) / 2; 02302 else 02303 tot_assys = ((nrings * (nrings - 1)) / 2) 02304 + (nrings + 1) / 2; 02305 } 02306 } while (reading_lattice == false) ; 02307 02308 // now reading the arrangement 02309 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02310 MAXCHARS, comment)) 02311 ERRORR("Reading input file failed",1); 02312 std::istringstream formatString2(input_string); 02313 for (int i = 1; i <= tot_assys; i++) { 02314 formatString2 >> temp_alias; 02315 if(formatString2.fail()) 02316 IOErrorHandler (INVALIDINPUT); 02317 core_alias.push_back(temp_alias); 02318 } 02319 } 02320 02321 else if (geom_type == "rectangular" && symm == 1) { 02322 02323 bool reading_assemblies = false; 02324 do{ 02325 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02326 MAXCHARS, comment)) 02327 ERRORR("Reading input file failed",1); 02328 std::istringstream formatString(input_string); 02329 if (input_string.substr(0, 10) == "assemblies"){ 02330 formatString >> card >> nassys >> pitchx >> pitchy; 02331 if(nassys < 0 || pitchx < 0 || pitchy< 0 || formatString.fail()) 02332 IOErrorHandler (INVALIDINPUT); 02333 02334 // reading file and alias names 02335 if(!parse_assembly_names(parse, argc, argv)) 02336 ERRORR("error parsing names of assemblies",1); 02337 reading_assemblies = true; 02338 } 02339 } while (reading_assemblies == false) ; 02340 02341 // reading lattice 02342 bool reading_lattice = false; 02343 do{ 02344 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02345 MAXCHARS, comment)) 02346 ERRORR("Reading input file failed",1); 02347 std::istringstream formatString1(input_string); 02348 if (input_string.substr(0,7) == "lattice"){ 02349 reading_lattice = true; 02350 formatString1 >> card >> nringsx >> nringsy; 02351 if(nringsx <= 0 || nringsy <= 0 || formatString1.fail()) 02352 IOErrorHandler (INVALIDINPUT); 02353 tot_assys = nringsx * nringsy; 02354 } 02355 } while (reading_lattice == false) ; 02356 02357 // now reading the arrangement 02358 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02359 MAXCHARS, comment)) 02360 ERRORR("Reading input file failed",1); 02361 std::istringstream formatString2(input_string); 02362 for (int i = 1; i <= tot_assys; i++) { 02363 formatString2 >> temp_alias; 02364 if(formatString2.fail()) 02365 IOErrorHandler (INVALIDINPUT); 02366 core_alias.push_back(temp_alias); 02367 } 02368 } 02369 02370 else if (geom_type == "hexflat" && symm == 6) { 02371 02372 // reading lattice 02373 bool reading_lattice = false; 02374 do{ 02375 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02376 MAXCHARS, comment)) 02377 ERRORR("Reading input file failed",1); 02378 std::istringstream formatString1(input_string); 02379 if (input_string.substr(0,7) == "lattice"){ 02380 reading_lattice = true; 02381 formatString1 >> card >> nrings; 02382 if(nrings < 0 || formatString1.fail()) 02383 IOErrorHandler (INVALIDINPUT); 02384 tot_assys = (nrings * (nrings + 1)) / 2; 02385 } 02386 } while (reading_lattice == false) ; 02387 02388 02389 // now reading the arrangement 02390 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02391 MAXCHARS, comment)) 02392 ERRORR("Reading input file failed",1); 02393 std::istringstream formatString2(input_string); 02394 for (int i = 1; i <= tot_assys; i++) { 02395 formatString2 >> temp_alias; 02396 if(formatString2.fail()) 02397 IOErrorHandler (INVALIDINPUT); 02398 core_alias.push_back(temp_alias); 02399 } 02400 } 02401 02402 else if (geom_type == "hexflat" && symm == 1) { 02403 // reading lattice 02404 bool reading_lattice = false; 02405 do{ 02406 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02407 MAXCHARS, comment)) 02408 ERRORR("Reading input file failed",1); 02409 std::istringstream formatString1(input_string); 02410 if (input_string.substr(0,7) == "lattice"){ 02411 reading_lattice = true; 02412 formatString1 >> card >> nrings; 02413 if(nrings < 0 || formatString1.fail()) 02414 IOErrorHandler (INVALIDINPUT); 02415 tot_assys = 3 * (nrings * (nrings - 1)) + 1; 02416 } 02417 } while (reading_lattice == false) ; 02418 02419 // now reading the arrangement 02420 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02421 MAXCHARS, comment)) 02422 ERRORR("Reading input file failed",1); 02423 std::istringstream formatString2(input_string); 02424 for (int i = 1; i <= tot_assys; i++) { 02425 formatString2 >> temp_alias; 02426 if(formatString2.fail()) 02427 IOErrorHandler (INVALIDINPUT); 02428 core_alias.push_back(temp_alias); 02429 } 02430 } 02431 02432 else if (geom_type == "hexflat" && symm == 12) { 02433 // reading lattice 02434 bool reading_lattice = false; 02435 do{ 02436 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02437 MAXCHARS, comment)) 02438 ERRORR("Reading input file failed",1); 02439 std::istringstream formatString1(input_string); 02440 if (input_string.substr(0,7) == "lattice"){ 02441 reading_lattice = true; 02442 formatString1 >> card >> nrings; 02443 if(nrings < 0 || formatString.fail()) 02444 IOErrorHandler (INVALIDINPUT); 02445 if (nrings % 2 == 0) 02446 tot_assys = (nrings * (nrings + 2)) / 4; 02447 else 02448 tot_assys = ((nrings + 1) * (nrings + 1)) / 4; 02449 } 02450 } while (reading_lattice == false) ; 02451 02452 // now reading the arrangement 02453 if (!parse.ReadNextLine(file_input, linenumber, input_string, 02454 MAXCHARS, comment)) 02455 ERRORR("Reading input file failed",1); 02456 std::istringstream formatString2(input_string); 02457 for (int i = 1; i <= tot_assys; i++) { 02458 formatString2 >> temp_alias; 02459 if(formatString2.fail()) 02460 IOErrorHandler (INVALIDINPUT); 02461 core_alias.push_back(temp_alias); 02462 } 02463 } 02464 02465 else { 02466 ERRORR("Invalid geometry type",1); 02467 } 02468 } 02469 // background mesh 02470 if (input_string.substr(0, 10) == "background") { 02471 std::istringstream formatString(input_string); 02472 formatString >> card >> back_meshfile; 02473 if(formatString.fail()) 02474 IOErrorHandler (INVALIDINPUT); 02475 02476 all_meshfiles.push_back(back_meshfile); 02477 02478 if (iname == COREGEN_DEFAULT_TEST_FILE){ 02479 back_meshfile = testdir + back_meshfile; 02480 } 02481 files.push_back(back_meshfile); 02482 back_mesh = true; 02483 } 02484 // z-height and z-divisions 02485 if (input_string.substr(0, 7) == "extrude") { 02486 extrude_flag = true; 02487 std::istringstream formatString(input_string); 02488 formatString >> card >> z_height >> z_divisions; 02489 if(z_divisions < 0 || formatString.fail()) 02490 IOErrorHandler (INVALIDINPUT); 02491 } 02492 // if keyword MESHTOGEOM is specified set tags on material sets and write a meshtogeom.txt file with #material id and #meshtogeom ratio 02493 if (input_string.substr(0, 10) == "meshtogeom") { 02494 compute_meshtogeom = true; 02495 } 02496 // OutputFileName 02497 if (input_string.substr(0, 14) == "outputfilename") { 02498 std::istringstream formatString(input_string); 02499 formatString >> card >> outfile; 02500 if(formatString.fail()) 02501 IOErrorHandler (INVALIDINPUT); 02502 } 02503 02504 // breaking condition 02505 if (input_string.substr(0, 3) == "end") { 02506 std::istringstream formatstring(input_string); 02507 break; 02508 } 02509 } 02510 // set some variables 02511 assm_meshfiles.resize(nassys); 02512 load_per_assm.resize(nassys); 02513 size_mf.resize(nassys); 02514 times_loaded.resize(nassys); 02515 assm_location.resize(nassys); 02516 bsameas.resize(nassys); 02517 02518 for(int i = 0; i < tot_assys; i++){ 02519 for (int j = 0; j < nassys; j++){ 02520 if (strcmp(core_alias[i].c_str(), assm_alias[j].c_str()) == 0) { 02521 assm_meshfiles[j]+=1; 02522 assm_location[j].push_back(i); 02523 break; 02524 } 02525 } 02526 } 02527 return iBase_SUCCESS; 02528 } 02529 02530 int CoreGen::parse_assembly_names(CParser parse, int argc, char *argv[] ) 02531 // --------------------------------------------------------------------------- 02532 // Function: Reads all the assemblies from CoreGen input file 02533 // Input: none 02534 // Output: none 02535 // --------------------------------------------------------------------------- 02536 { 02537 02538 // reading file and alias names 02539 for (int i = 1; i <= nassys; i++) { 02540 if (!parse.ReadNextLine(file_input, linenumber, 02541 input_string, MAXCHARS, comment, false)) 02542 ERRORR("Reading input file failed",1); 02543 std::istringstream formatString(input_string); 02544 formatString >> meshfile >> mf_alias >> same_as >> reloading_mf >> ms_startid >> ns_startid; 02545 // we don't check for formatting since same_as and parameters after it may not be present. 02546 // variable gets populated correctly in the file 02547 02548 // if meshfile variable is a path then only convert the filename to lower case 02549 unsigned pos = 0; 02550 pos = meshfile.find_last_of("/\\"); 02551 if (pos > 0 && pos < meshfile.length()){ 02552 std::string filename = "", path=""; 02553 path = meshfile.substr(0,pos); 02554 filename = meshfile.substr(pos+1); 02555 std::transform(filename.begin(), filename.end(), filename.begin(), ::tolower); 02556 meshfile = path+"/"+ filename; 02557 } 02558 else{ 02559 std::transform(meshfile.begin(), meshfile.end(), meshfile.begin(), ::tolower); 02560 } 02561 // also convert the alias and reloading_mf name to lower case, since we've changed the actual reading. 02562 std::transform(mf_alias.begin(), mf_alias.end(), mf_alias.begin(), ::tolower); 02563 std::transform(reloading_mf.begin(), reloading_mf.end(), reloading_mf.begin(), ::tolower); 02564 02565 if (same_as == "same_as") 02566 bsameas.push_back(0); 02567 else 02568 bsameas.push_back(1); 02569 if(bsameas[i-1] == 1){ 02570 all_meshfiles.push_back(meshfile); 02571 if (argc == 1){ 02572 meshfile = testdir + "/" + meshfile; 02573 } 02574 files.push_back(meshfile); 02575 assm_alias.push_back(mf_alias); 02576 all_ms_starts.push_back(-1); 02577 all_ns_starts.push_back(-1); 02578 } 02579 else{ 02580 all_meshfiles.push_back(reloading_mf); 02581 if (iname == COREGEN_DEFAULT_TEST_FILE){ 02582 meshfile = testdir + reloading_mf; 02583 } 02584 files.push_back(reloading_mf); 02585 assm_alias.push_back(mf_alias); 02586 all_ms_starts.push_back(ms_startid); 02587 all_ns_starts.push_back(ns_startid); 02588 } 02589 // bsameas = false; 02590 } 02591 return 0; 02592 } 02593 02594 int CoreGen::find_assm(const int i, int &assm_index) 02595 // --------------------------------------------------------------------------- 02596 // Function: find the assembly index (0 to n) for n assemblies for core alias i 02597 // Input: none 02598 // Output: none 02599 // --------------------------------------------------------------------------- 02600 { 02601 int flag = 0; 02602 for (int j = 0; j < nassys; j++) 02603 if (strcmp(core_alias[i].c_str(), assm_alias[j].c_str()) == 0) { 02604 assm_index = j; 02605 flag = 1; 02606 break; 02607 } 02608 if (flag == 0)//nothing found return -1 or no assembly 02609 assm_index = -1; 02610 return iBase_SUCCESS; 02611 } 02612 02613 int CoreGen::write_makefile() 02614 // --------------------------------------------------------------------------- 02615 // Function: write the makefile based on inputs read from input file 02616 // Input: none 02617 // Output: none 02618 // --------------------------------------------------------------------------- 02619 { 02620 std::string name; 02621 std::vector<std::string> f_no_ext, f_sat, f_inp, f_jou, f_injou, f_py; 02622 make_file << "##" << std::endl; 02623 make_file 02624 << "## This makefile is automatically generated by coregen program" 02625 << std::endl; 02626 make_file << "##" << std::endl; 02627 make_file << "## Check your coregen, assygen and cubit location" 02628 << std::endl; 02629 make_file << "##" << std::endl; 02630 make_file << "include ../../common.mk\n" << std::endl; 02631 02632 make_file << "\nCUBIT = trelis -nographics\n" << std::endl; 02633 make_file << "COREGEN = $(MESHKIT_DIR)/bin/coregen\n" << std::endl; 02634 make_file << "ASSYGEN = $(MESHKIT_DIR)/bin/assygen\n" << std::endl; 02635 02636 // remove the ./ if run from the current working TestDirectory 02637 make_file << "MESH_FILES = "; 02638 std::string filename; 02639 for (unsigned int i = 0; i < files.size(); i++) { 02640 if (files[i][0] == '.' && files[i][1] == '/') { 02641 filename = files[i].substr(2, files[i].length()); 02642 } else if (files[i][0] != '.' || files[i][1] != '/') { 02643 int loc1 = files[i].find_last_of("."); 02644 int loc2 = files[i].find_last_of("/"); 02645 filename = files[i].substr(loc2 + 1, loc1); 02646 } else { 02647 filename = files[i]; 02648 } 02649 mk_files.push_back(filename); 02650 make_file << all_meshfiles[i] << " "; 02651 } 02652 02653 // get file names without extension 02654 for (unsigned int i = 0; i < mk_files.size(); i++) { 02655 int loc = mk_files[i].find_first_of("."); 02656 f_no_ext.push_back(mk_files[i].substr(0, loc)); 02657 } 02658 02659 make_file << "\n\nGEOM_FILES = "; 02660 for (unsigned int i = 0; i < mk_files.size(); i++) { 02661 if (geom_engine == "occ") 02662 name = f_no_ext[i] + ".stp"; 02663 else 02664 name = f_no_ext[i] + ".sat"; 02665 f_sat.push_back(name); 02666 make_file << name << " "; 02667 name = ""; 02668 } 02669 02670 make_file << "\n\nJOU_FILES = "; 02671 for (unsigned int i = 0; i < mk_files.size(); i++) { 02672 name = f_no_ext[i] + ".jou"; 02673 f_jou.push_back(name); 02674 make_file << name << " "; 02675 name = ""; 02676 } 02677 02678 make_file << "\n\nPYTHON_FILES = "; 02679 for (unsigned int i = 0; i < mk_files.size(); i++) { 02680 name = f_no_ext[i] + ".py"; 02681 f_py.push_back(name); 02682 make_file << name << " "; 02683 name = ""; 02684 } 02685 02686 make_file << "\n\nINJOU_FILES = "; 02687 for (unsigned int i = 0; i < mk_files.size(); i++) { 02688 name = f_no_ext[i] + ".template.jou"; 02689 f_injou.push_back(name); 02690 make_file << name << " "; 02691 name = ""; 02692 } 02693 02694 make_file << "\n\nASSYGEN_FILES = "; 02695 for (unsigned int i = 0; i < mk_files.size(); i++) { 02696 name = f_no_ext[i] + ".inp"; 02697 f_inp.push_back(name); 02698 make_file << name << " "; 02699 name = ""; 02700 } 02701 02702 make_file << "\n\n" << outfile << " : ${MESH_FILES} " << ifile << std::endl; 02703 make_file << "\t" << "${COREGEN} " << iname << std::endl; 02704 for (unsigned int i = 0; i < mk_files.size(); i++) { 02705 make_file << all_meshfiles[i] << " : " << f_sat[i] << " " << f_jou[i] 02706 << " " << f_injou[i] << std::endl; 02707 make_file << "\t" << "${CUBIT} -batch " << f_jou[i] << "\n" 02708 << std::endl; 02709 02710 make_file << f_sat[i] << " : " << f_py[i] << std::endl; 02711 make_file << "\t" << "${CUBIT} -batch " << f_py[i] << "\n" 02712 << std::endl; 02713 02714 make_file << f_py[i] << " " << f_jou[i] << " " << f_injou[i] << " : " 02715 << f_inp[i] << std::endl; 02716 make_file << "\t" << "${ASSYGEN} " << f_no_ext[i] << "\n" << std::endl; 02717 } 02718 02719 make_file << "make clean :" << "\n\trm *.sat *.cub *.py *.jou" << std::endl; 02720 make_file.close(); 02721 logfile << "Created makefile: " << mfile << std::endl; 02722 return 0; 02723 } 02724 02725 void CoreGen::banner() 02726 // --------------------------------------------------------------------------- 02727 // Function: display the program banner 02728 // Input: none 02729 // Output: none 02730 // --------------------------------------------------------------------------- 02731 { 02732 if(rank == 0 ){ 02733 logfile << '\n'; 02734 logfile 02735 << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX" 02736 << '\n'; 02737 logfile 02738 << "Program to Assemble Nuclear Reactor Assembly Meshes and Form a Core " 02739 << '\n'; 02740 logfile << "\t\t\tArgonne National Laboratory" << '\n'; 02741 logfile << "\t\t\t 2015 " << '\n'; 02742 logfile 02743 << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX" 02744 << '\n'; 02745 } 02746 } 02747 02748 } // namespace MeshKit
1.7.6.1 
