moab
|
00001 #include "moab/ParallelComm.hpp" 00002 #include "MBParallelConventions.h" 00003 #include "moab/Core.hpp" 00004 #include "Coupler.hpp" 00005 #include "iMesh_extensions.h" 00006 #include "moab_mpi.h" 00007 #include "ElemUtil.hpp" 00008 #include <iostream> 00009 #include <iomanip> 00010 #include <sstream> 00011 #include <assert.h> 00012 #include "MBiMesh.hpp" 00013 00014 using namespace moab; 00015 00016 #define STRINGIFY_(A) #A 00017 #define STRINGIFY(A) STRINGIFY_(A) 00018 #ifdef MESHDIR 00019 std::string TestDir( STRINGIFY(MESHDIR) ); 00020 #else 00021 std::string TestDir("."); 00022 #endif 00023 00024 #define RRA(a) if (MB_SUCCESS != result) {\ 00025 std::string tmp_str; mbImpl->get_last_error(tmp_str);\ 00026 tmp_str.append("\n"); tmp_str.append(a);\ 00027 dynamic_cast<Core*>(mbImpl)->get_error_handler()->set_last_error(tmp_str.c_str()); \ 00028 return result;} 00029 00030 #define PRINT_LAST_ERROR \ 00031 if (MB_SUCCESS != result) {\ 00032 std::string tmp_str;\ 00033 std::cout << "Failure; message:" << std::endl;\ 00034 mbImpl->get_last_error(tmp_str);\ 00035 std::cout << tmp_str << std::endl;\ 00036 MPI_Abort(MPI_COMM_WORLD, result); \ 00037 return result;\ 00038 } 00039 00040 #define PRINT_LAST_ERR \ 00041 if (iBase_SUCCESS != err) {\ 00042 std::string tmp_str;\ 00043 std::cout << "Failure; message:" << std::endl;\ 00044 mbImpl->get_last_error(tmp_str);\ 00045 std::cout << tmp_str << std::endl;\ 00046 MPI_Abort(MPI_COMM_WORLD, result); \ 00047 return result;\ 00048 } 00049 00050 void print_usage(); 00051 00052 ErrorCode get_file_options(int argc, char **argv, 00053 std::vector<std::string> &meshFiles, 00054 Coupler::Method &method, 00055 std::string &interpTag, 00056 std::string &gNormTag, 00057 std::string &ssNormTag, 00058 std::vector<const char *> &ssTagNames, 00059 std::vector<const char *> &ssTagValues, 00060 std::string &readOpts, 00061 std::string &outFile, 00062 std::string &writeOpts, 00063 std::string &dbgFile, 00064 bool &help, 00065 double & epsilon); 00066 00067 // ErrorCode get_file_options(int argc, char **argv, 00068 // std::vector<const char *> &filenames, 00069 // std::string &tag_name, 00070 // std::string &out_fname, 00071 // std::string &opts); 00072 00073 ErrorCode report_iface_ents(Interface *mbImpl, 00074 std::vector<ParallelComm *> &pcs, 00075 bool print_results); 00076 00077 ErrorCode test_interpolation(Interface *mbImpl, 00078 Coupler::Method method, 00079 std::string &interpTag, 00080 std::string &gNormTag, 00081 std::string &ssNormTag, 00082 std::vector<const char *> &ssTagNames, 00083 std::vector<const char *> &ssTagValues, 00084 iBase_EntitySetHandle *roots, 00085 std::vector<ParallelComm *> &pcs, 00086 double &instant_time, 00087 double &pointloc_time, 00088 double &interp_time, 00089 double &gnorm_time, 00090 double &ssnorm_time, 00091 double & toler); 00092 00093 void reduceMax(double &v){ 00094 double buf; 00095 00096 MPI_Barrier(MPI_COMM_WORLD); 00097 MPI_Allreduce( &v, &buf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); 00098 00099 v=buf; 00100 } 00101 00102 00103 int main(int argc, char **argv) 00104 { 00105 // need to init MPI first, to tell how many procs and rank 00106 int err = MPI_Init(&argc, &argv); 00107 00108 std::vector<const char *> ssTagNames, ssTagValues; 00109 std::vector<std::string> meshFiles; 00110 std::string interpTag, gNormTag, ssNormTag, readOpts, outFile, writeOpts, dbgFile; 00111 Coupler::Method method = Coupler::CONSTANT; 00112 00113 ErrorCode result = MB_SUCCESS; 00114 bool help = false; 00115 double toler = 5.e-10; 00116 result = get_file_options(argc, argv, meshFiles, method, interpTag, 00117 gNormTag, ssNormTag, ssTagNames, ssTagValues, 00118 readOpts, outFile, writeOpts, dbgFile, help, toler); 00119 00120 if (result != MB_SUCCESS || help) { 00121 print_usage(); 00122 err = MPI_Finalize(); 00123 return 1; 00124 } 00125 00126 int nprocs, rank; 00127 err = MPI_Comm_size(MPI_COMM_WORLD, &nprocs); 00128 err = MPI_Comm_rank(MPI_COMM_WORLD, &rank); 00129 00130 // redirect stdout and stderr if dbgFile is not null 00131 if (!dbgFile.empty()) { 00132 std::stringstream dfname; 00133 dfname << dbgFile << rank << ".txt"; 00134 if (!std::freopen(dfname.str().c_str(), "a", stdout)) return false; 00135 if (!std::freopen(dfname.str().c_str(), "a", stderr)) return false; 00136 } 00137 00138 // create MOAB instance based on that 00139 Interface *mbImpl = new Core(); 00140 if (NULL == mbImpl) return 1; 00141 00142 iMesh_Instance iMeshInst = reinterpret_cast<iMesh_Instance>( new MBiMesh(mbImpl) ); 00143 00144 // read in mesh(es) 00145 std::vector<ParallelComm *> pcs(meshFiles.size()); 00146 00147 // Create root sets for each mesh using the iMesh API. Then pass these 00148 // to the load_file functions to be populated. 00149 iBase_EntitySetHandle *roots = (iBase_EntitySetHandle *) malloc(sizeof(iBase_EntitySetHandle) * meshFiles.size()); 00150 00151 for (unsigned int i = 0; i < meshFiles.size(); i++) { 00152 pcs[i] = new ParallelComm(mbImpl, MPI_COMM_WORLD); 00153 int index = pcs[i]->get_id(); 00154 std::string newReadopts; 00155 std::ostringstream extraOpt; 00156 extraOpt << ";PARALLEL_COMM=" << index; 00157 newReadopts = readOpts+extraOpt.str(); 00158 00159 iMesh_createEntSet(iMeshInst, 0, &(roots[i]), &err); 00160 result = mbImpl->load_file( meshFiles[i].c_str(), (EntityHandle *)&roots[i], newReadopts.c_str() ); 00161 //result = rps[i]->load_file(meshFiles[i].c_str(), (EntityHandle *)&roots[i], FileOptions(readOpts.c_str())); 00162 PRINT_LAST_ERROR; 00163 } 00164 00165 result = report_iface_ents(mbImpl, pcs, true); 00166 PRINT_LAST_ERROR; 00167 00168 double instant_time=0.0, pointloc_time=0.0, interp_time=0.0, gnorm_time=0.0, ssnorm_time=0.0; 00169 // test interpolation and global normalization and subset normalization 00170 00171 result = test_interpolation(mbImpl, method, interpTag, gNormTag, ssNormTag, 00172 ssTagNames, ssTagValues, roots, pcs, 00173 instant_time, pointloc_time, interp_time, 00174 gnorm_time, ssnorm_time, toler); 00175 PRINT_LAST_ERROR; 00176 00177 00178 reduceMax(instant_time); 00179 reduceMax(pointloc_time); 00180 reduceMax(interp_time); 00181 00182 if(0==rank) 00183 printf("\nMax time : %g %g %g (inst loc interp -- %d procs )\n", instant_time, 00184 pointloc_time, interp_time, nprocs); 00185 00186 // output mesh 00187 if (!outFile.empty()) { 00188 Range partSets; 00189 // only save the target mesh 00190 partSets.insert((EntityHandle)roots[1]); 00191 std::string newwriteOpts; 00192 std::ostringstream extraOpt; 00193 extraOpt << ";PARALLEL_COMM=" << 1; 00194 newwriteOpts = writeOpts+extraOpt.str(); 00195 result = mbImpl->write_file(outFile.c_str(), NULL, newwriteOpts.c_str(), partSets); 00196 PRINT_LAST_ERROR; 00197 std::cout << "Wrote " << outFile << std::endl; 00198 std::cout << "mbcoupler_test complete." << std::endl; 00199 } 00200 00201 for (unsigned int i = 0; i < meshFiles.size(); i++) { 00202 delete pcs[i]; 00203 } 00204 00205 delete mbImpl; 00206 //may be leaking iMeshInst, don't care since it's end of program. Remove above deletes? 00207 00208 err = MPI_Finalize(); 00209 00210 return 0; 00211 } 00212 00213 ErrorCode report_iface_ents(Interface *mbImpl, 00214 std::vector<ParallelComm *> &pcs, 00215 const bool print_results) 00216 { 00217 Range iface_ents[6]; 00218 ErrorCode result = MB_SUCCESS, tmp_result; 00219 00220 // now figure out which vertices are shared 00221 for (unsigned int p = 0; p < pcs.size(); p++) { 00222 for (int i = 0; i < 4; i++) { 00223 tmp_result = pcs[p]->get_iface_entities(-1, i, iface_ents[i]); 00224 00225 if (MB_SUCCESS != tmp_result) { 00226 std::cerr << "get_iface_entities returned error on proc " 00227 << pcs[p]->proc_config().proc_rank() << "; message: " << std::endl; 00228 std::string last_error; 00229 result = mbImpl->get_last_error(last_error); 00230 if (last_error.empty()) std::cerr << "(none)" << std::endl; 00231 else std::cerr << last_error << std::endl; 00232 result = tmp_result; 00233 } 00234 if (0 != i) iface_ents[4].merge(iface_ents[i]); 00235 } 00236 } 00237 00238 // report # iface entities 00239 result = mbImpl->get_adjacencies(iface_ents[4], 0, false, iface_ents[5], 00240 Interface::UNION); 00241 00242 int rank; 00243 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 00244 00245 if (print_results || iface_ents[0].size() != iface_ents[5].size()) { 00246 std::cerr << "Proc " << rank << " iface entities: " << std::endl; 00247 for (int i = 0; i < 4; i++) 00248 std::cerr << " " << iface_ents[i].size() << " " 00249 << i << "d iface entities." << std::endl; 00250 std::cerr << " (" << iface_ents[5].size() 00251 << " verts adj to other iface ents)" << std::endl; 00252 } 00253 00254 return result; 00255 } 00256 00257 // Print usage 00258 void print_usage() { 00259 std::cerr << "Usage: "; 00260 std::cerr << "mbcoupler_test -meshes <source_mesh> <target_mesh> -itag <interp_tag> [-gnorm <gnorm_tag>] [-ssnorm <ssnorm_tag> <ssnorm_selection>] [-ropts <roptions>] [-outfile <out_file> [-wopts <woptions>]] [-dbgout [<dbg_file>]]" << std::endl; 00261 std::cerr << " -meshes" << std::endl; 00262 std::cerr << " Read in mesh files <source_mesh> and <target_mesh>." << std::endl; 00263 std::cerr << " -itag" << std::endl; 00264 std::cerr << " Interpolate tag <interp_tag> from source mesh to target mesh." << std::endl; 00265 std::cerr << " -gnorm" << std::endl; 00266 std::cerr << " Normalize the value of tag <gnorm_tag> over then entire mesh and save to" << std::endl; 00267 std::cerr << " tag \"<gnorm_tag>_normf\" on the mesh set. Do this for all meshes." << std::endl; 00268 std::cerr << " -ssnorm" << std::endl; 00269 std::cerr << " Normalize the value of tag <ssnorm_tag> over subsets of a mesh and save to" << std::endl; 00270 std::cerr << " tag \"<ssnorm_tag>_normf\" on the Entity Set for each subset. Subsets are selected" << std::endl; 00271 std::cerr << " using criteria in <ssnorm_selection>. Do this for all meshes." << std::endl; 00272 std::cerr << " -ropts" << std::endl; 00273 std::cerr << " Read in the mesh files using options in <roptions>." << std::endl; 00274 std::cerr << " -outfile" << std::endl; 00275 std::cerr << " Write out target mesh to <out_file>." << std::endl; 00276 std::cerr << " -wopts" << std::endl; 00277 std::cerr << " Write out mesh files using options in <woptions>." << std::endl; 00278 std::cerr << " -dbgout" << std::endl; 00279 std::cerr << " Write stdout and stderr streams to the file \'<dbg_file>.txt\'." << std::endl; 00280 std::cerr << " -eps" << std::endl; 00281 std::cerr << " epsilon" << std::endl; 00282 std::cerr << " -meth <method> (0=CONSTANT, 1=LINEAR_FE, 2=QUADRATIC_FE, 3=SPECTRAL)" << std::endl; 00283 } 00284 00285 // Check first character for a '-'. 00286 // Return true if one is found. False otherwise. 00287 bool check_for_flag(const char *str) { 00288 if (str[0] == '-') 00289 return true; 00290 else 00291 return false; 00292 } 00293 00294 // New get_file_options() function with added possibilities for mbcoupler_test. 00295 ErrorCode get_file_options(int argc, char **argv, 00296 std::vector<std::string> &meshFiles, 00297 Coupler::Method &method, 00298 std::string &interpTag, 00299 std::string &gNormTag, 00300 std::string &ssNormTag, 00301 std::vector<const char *> &ssTagNames, 00302 std::vector<const char *> &ssTagValues, 00303 std::string &readOpts, 00304 std::string &outFile, 00305 std::string &writeOpts, 00306 std::string &dbgFile, 00307 bool &help, 00308 double & epsilon) 00309 { 00310 // Initialize some of the outputs to null values indicating not present 00311 // in the argument list. 00312 gNormTag = ""; 00313 ssNormTag = ""; 00314 readOpts = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=3.0.1;CPUTIME"; 00315 outFile = ""; 00316 writeOpts = "PARALLEL=WRITE_PART;CPUTIME"; 00317 dbgFile = ""; 00318 std::string defaultDbgFile = argv[0]; // The executable name will be the default debug output file. 00319 00320 // These will indicate if we've gotten our required parameters at the end of parsing. 00321 bool haveMeshes = false; 00322 bool haveInterpTag = false; 00323 00324 // Loop over the values in argv pulling out an parsing each one 00325 int npos = 1; 00326 00327 if (argc > 1 && argv[1] == std::string("-h")) { 00328 help = true; 00329 return MB_SUCCESS; 00330 } 00331 00332 while (npos < argc) { 00333 if (argv[npos] == std::string("-meshes")) { 00334 // Parse out the mesh filenames 00335 npos++; 00336 int numFiles = 2; 00337 meshFiles.resize(numFiles); 00338 for (int i = 0; i < numFiles; i++) { 00339 if ((npos < argc) && (!check_for_flag(argv[npos]))) 00340 meshFiles[i] = argv[npos++]; 00341 else { 00342 std::cerr << " ERROR - missing correct number of mesh filenames" << std::endl; 00343 return MB_FAILURE; 00344 } 00345 } 00346 00347 haveMeshes = true; 00348 } 00349 else if (argv[npos] == std::string("-itag")) { 00350 // Parse out the interpolation tag 00351 npos++; 00352 if ((npos < argc) && (!check_for_flag(argv[npos]))) 00353 interpTag = argv[npos++]; 00354 else { 00355 std::cerr << " ERROR - missing <interp_tag>" << std::endl; 00356 return MB_FAILURE; 00357 } 00358 00359 haveInterpTag = true; 00360 } 00361 else if (argv[npos] == std::string("-meth")) { 00362 // Parse out the interpolation tag 00363 npos++; 00364 if (argv[npos][0] == '0') method = Coupler::CONSTANT; 00365 else if (argv[npos][0] == '1') method = Coupler::LINEAR_FE; 00366 else if (argv[npos][0] == '2') method = Coupler::QUADRATIC_FE; 00367 else if (argv[npos][0] == '3') method = Coupler::SPECTRAL; 00368 else { 00369 std::cerr << " ERROR - unrecognized method number " << method << std::endl; 00370 return MB_FAILURE; 00371 } 00372 npos++; 00373 } 00374 else if (argv[npos] == std::string("-eps")) { 00375 // Parse out the tolerance 00376 npos++; 00377 if ((npos < argc) && (!check_for_flag(argv[npos]))) 00378 epsilon = atof(argv[npos++]); 00379 else { 00380 std::cerr << " ERROR - missing <epsilon>" << std::endl; 00381 return MB_FAILURE; 00382 } 00383 } 00384 else if (argv[npos] == std::string("-gnorm")) { 00385 // Parse out the global normalization tag 00386 npos++; 00387 if ((npos < argc) && (!check_for_flag(argv[npos]))) 00388 gNormTag = argv[npos++]; 00389 else { 00390 std::cerr << " ERROR - missing <gnorm_tag>" << std::endl; 00391 return MB_FAILURE; 00392 } 00393 } 00394 else if (argv[npos] == std::string("-ssnorm")) { 00395 // Parse out the subset normalization tag and selection criteria 00396 npos++; 00397 if ((npos < argc) && (!check_for_flag(argv[npos]))) 00398 ssNormTag = argv[npos++]; 00399 else { 00400 std::cerr << " ERROR - missing <ssnorm_tag>" << std::endl; 00401 return MB_FAILURE; 00402 } 00403 00404 if ((npos < argc) && (!check_for_flag(argv[npos]))) { 00405 char* opts = argv[npos++]; 00406 char sep1[1] = {';'}; 00407 char sep2[1] = {'='}; 00408 bool end_vals_seen = false; 00409 std::vector<char *> tmpTagOpts; 00410 00411 // first get the options 00412 for (char* i = strtok(opts, sep1); i; i = strtok(0, sep1)) { 00413 tmpTagOpts.push_back(i); 00414 } 00415 00416 // parse out the name and val or just name. 00417 for (unsigned int j = 0; j < tmpTagOpts.size(); j++) { 00418 char* e = strtok(tmpTagOpts[j], sep2); 00419 ssTagNames.push_back(e); 00420 e = strtok(0, sep2); 00421 if (e != NULL) { 00422 // We have a value 00423 if (end_vals_seen) { 00424 // ERROR we should not have a value after none are seen 00425 std::cerr << " ERROR - new value seen after end of values in <ssnorm_selection>" << std::endl; 00426 return MB_FAILURE; 00427 } 00428 // Otherwise get the value string from e and convert it to an int 00429 int *valp = new int; 00430 *valp = atoi(e); 00431 ssTagValues.push_back((const char *) valp); 00432 } 00433 else { 00434 // Otherwise there is no '=' so push a null on the list 00435 end_vals_seen = true; 00436 ssTagValues.push_back((const char *) 0); 00437 } 00438 } 00439 } 00440 else { 00441 std::cerr << " ERROR - missing <ssnorm_selection>" << std::endl; 00442 return MB_FAILURE; 00443 } 00444 } 00445 else if (argv[npos] == std::string("-ropts")) { 00446 // Parse out the mesh file read options 00447 npos++; 00448 if ((npos < argc) && (!check_for_flag(argv[npos]))) 00449 readOpts = argv[npos++]; 00450 else { 00451 std::cerr << " ERROR - missing <roptions>" << std::endl; 00452 return MB_FAILURE; 00453 } 00454 } 00455 else if (argv[npos] == std::string("-outfile")) { 00456 // Parse out the output file name 00457 npos++; 00458 if ((npos < argc) && (!check_for_flag(argv[npos]))) 00459 outFile = argv[npos++]; 00460 else { 00461 std::cerr << " ERROR - missing <out_file>" << std::endl; 00462 return MB_FAILURE; 00463 } 00464 } 00465 else if (argv[npos] == std::string("-wopts")) { 00466 // Parse out the output file write options 00467 npos++; 00468 if ((npos < argc) && (!check_for_flag(argv[npos]))) 00469 writeOpts = argv[npos++]; 00470 else { 00471 std::cerr << " ERROR - missing <woptions>" << std::endl; 00472 return MB_FAILURE; 00473 } 00474 } 00475 else if (argv[npos] == std::string("-dbgout")) { 00476 // Parse out the debug output file name. 00477 // If no name then use the default. 00478 npos++; 00479 if ((npos < argc) && (!check_for_flag(argv[npos]))) 00480 dbgFile = argv[npos++]; 00481 else 00482 dbgFile = defaultDbgFile; 00483 } 00484 else { 00485 // Unrecognized parameter. Skip it and move along. 00486 std::cerr << " ERROR - Unrecognized parameter:" << argv[npos] << std::endl; 00487 std::cerr << " Skipping..." << std::endl; 00488 npos++; 00489 } 00490 } 00491 00492 if (!haveMeshes) { 00493 meshFiles.resize(2); 00494 meshFiles[0] = std::string(TestDir + "/64bricks_1khex.h5m"); 00495 meshFiles[1] = std::string(TestDir + "/64bricks_12ktet.h5m"); 00496 std::cout << "Mesh files not entered; using default files " 00497 << meshFiles[0] << " and " << meshFiles[1] << std::endl; 00498 } 00499 00500 if (!haveInterpTag) { 00501 interpTag = "vertex_field"; 00502 std::cout << "Interpolation field name not given, using default of " << interpTag << std::endl; 00503 } 00504 00505 #ifdef HDF5_FILE 00506 if (1 == argc) { 00507 std::cout << "No arguments given; using output file dum.h5m." << std::endl; 00508 outFile = "dum.h5m"; 00509 } 00510 #endif 00511 00512 return MB_SUCCESS; 00513 } 00514 00515 00516 // End new get_file_options() 00517 00518 // ErrorCode get_file_options(int argc, char **argv, 00519 // std::vector<const char *> &filenames, 00520 // std::string &interp_tag, 00521 // std::string &out_fname, 00522 // std::string &opts) 00523 // { 00524 // int npos = 1; 00525 // int nfiles = atoi(argv[npos++]); 00526 00527 // // get mesh filenames 00528 // filenames.resize(nfiles); 00529 // for (int i = 0; i < nfiles; i++) filenames[i] = argv[npos++]; 00530 00531 // if (npos < argc) interp_tag = argv[npos++]; 00532 00533 // if (npos < argc) out_fname = argv[npos++]; 00534 00535 // // get partition information 00536 // const char *tag_name = "GEOM_DIMENSION"; 00537 // int tag_val = -1; 00538 // int distrib = 1; 00539 // int with_ghosts = 0; 00540 // if (npos < argc) tag_name = argv[npos++]; 00541 // if (npos < argc) tag_val = strtol(argv[npos++], NULL, 0); 00542 // if (npos < argc) distrib = strtol(argv[npos++], NULL, 0); 00543 // if (npos < argc) with_ghosts = strtol(argv[npos++], NULL, 0); 00544 00545 // if (-1 == tag_val && !strcmp(tag_name, "GEOM_DIMENSION")) 00546 // tag_val = 3; 00547 00548 // std::ostringstream options; 00549 // options << "PARALLEL=READ_DELETE;PARTITION=" << tag_name; 00550 00551 // if (-1 != tag_val) 00552 // options << ";PARTITION_VAL=" << tag_val; 00553 00554 // if (1 == distrib) 00555 // options << ";PARTITION_DISTRIBUTE"; 00556 00557 // options << ";PARALLEL_RESOLVE_SHARED_ENTS"; 00558 00559 // if (1 == with_ghosts) 00560 // options << ";PARALLEL_GHOSTS=3.0.1"; 00561 00562 // options << ";CPUTIME"; 00563 00564 // opts = options.str(); 00565 00566 // return MB_SUCCESS; 00567 // } 00568 00569 ErrorCode test_interpolation(Interface *mbImpl, 00570 Coupler::Method method, 00571 std::string &interpTag, 00572 std::string &gNormTag, 00573 std::string &ssNormTag, 00574 std::vector<const char *> &ssTagNames, 00575 std::vector<const char *> &ssTagValues, 00576 iBase_EntitySetHandle *roots, 00577 std::vector<ParallelComm *> &pcs, 00578 double &instant_time, 00579 double &pointloc_time, 00580 double &interp_time, 00581 double &gnorm_time, 00582 double &ssnorm_time, 00583 double & toler) 00584 { 00585 assert(method >= Coupler::CONSTANT && method <= Coupler::SPECTRAL); 00586 00587 // source is 1st mesh, target is 2nd 00588 Range src_elems, targ_elems, targ_verts; 00589 ErrorCode result = pcs[0]->get_part_entities(src_elems, 3); 00590 PRINT_LAST_ERROR; 00591 00592 double start_time = MPI_Wtime(); 00593 00594 // instantiate a coupler, which also initializes the tree 00595 Coupler mbc(mbImpl, pcs[0], src_elems, 0); 00596 00597 // initialize spectral elements, if they exist 00598 bool specSou=false, specTar = false; 00599 // result = mbc.initialize_spectral_elements((EntityHandle)roots[0], (EntityHandle)roots[1], specSou, specTar); 00600 00601 00602 instant_time = MPI_Wtime(); 00603 00604 // get points from the target mesh to interpolate 00605 // we have to treat differently the case when the target is a spectral mesh 00606 // in that case, the points of interest are the GL points, not the vertex nodes 00607 std::vector<double> vpos; // this will have the positions we are interested in 00608 int numPointsOfInterest = 0; 00609 if (!specTar) 00610 {// usual case 00611 Range tmp_verts; 00612 00613 // first get all vertices adj to partition entities in target mesh 00614 result = pcs[1]->get_part_entities(targ_elems, 3); 00615 PRINT_LAST_ERROR; 00616 if (Coupler::CONSTANT == method) 00617 targ_verts = targ_elems; 00618 else 00619 result = mbImpl->get_adjacencies(targ_elems, 0, false, targ_verts, 00620 Interface::UNION); 00621 PRINT_LAST_ERROR; 00622 00623 // then get non-owned verts and subtract 00624 result = pcs[1]->get_pstatus_entities(0, PSTATUS_NOT_OWNED, tmp_verts); 00625 PRINT_LAST_ERROR; 00626 targ_verts = subtract( targ_verts, tmp_verts); 00627 // get position of these entities; these are the target points 00628 numPointsOfInterest = (int)targ_verts.size(); 00629 vpos.resize(3*targ_verts.size()); 00630 result = mbImpl->get_coords(targ_verts, &vpos[0]); 00631 PRINT_LAST_ERROR; 00632 } 00633 else 00634 { 00635 // in this case, the target mesh is spectral, we want values 00636 // interpolated on the GL positions; for each element, get the GL points, and construct CartVect!!! 00637 result = pcs[1]->get_part_entities(targ_elems, 3); 00638 PRINT_LAST_ERROR; 00639 result = mbc.get_gl_points_on_elements(targ_elems, vpos, numPointsOfInterest); 00640 PRINT_LAST_ERROR; 00641 } 00642 00643 00644 00645 00646 // locate those points in the source mesh 00647 std::cout<<"rank "<< pcs[0]->proc_config().proc_rank() << " points of interest: " << numPointsOfInterest << "\n"; 00648 result = mbc.locate_points(&vpos[0], numPointsOfInterest, 0, toler); 00649 PRINT_LAST_ERROR; 00650 00651 pointloc_time = MPI_Wtime(); 00652 00653 // now interpolate tag onto target points 00654 std::vector<double> field(numPointsOfInterest); 00655 00656 result = mbc.interpolate(method, interpTag, &field[0]); 00657 PRINT_LAST_ERROR; 00658 00659 interp_time = MPI_Wtime(); 00660 00661 // do global normalization if specified 00662 if (!gNormTag.empty()) { 00663 int err; 00664 00665 // Normalize the source mesh 00666 err = mbc.normalize_mesh(roots[0], gNormTag.c_str(), Coupler::VOLUME, 4); 00667 PRINT_LAST_ERR; 00668 00669 // Normalize the target mesh 00670 err = mbc.normalize_mesh(roots[1], gNormTag.c_str(), Coupler::VOLUME, 4); 00671 PRINT_LAST_ERR; 00672 } 00673 00674 gnorm_time = MPI_Wtime(); 00675 00676 // do subset normalization if specified 00677 00678 if (!ssNormTag.empty()) { 00679 int err; 00680 00681 err = mbc.normalize_subset(roots[0], 00682 ssNormTag.c_str(), 00683 &ssTagNames[0], 00684 ssTagNames.size(), 00685 &ssTagValues[0], 00686 Coupler::VOLUME, 00687 4); 00688 PRINT_LAST_ERR; 00689 00690 err = mbc.normalize_subset(roots[1], 00691 ssNormTag.c_str(), 00692 &ssTagNames[0], 00693 ssTagNames.size(), 00694 &ssTagValues[0], 00695 Coupler::VOLUME, 00696 4); 00697 PRINT_LAST_ERR; 00698 } 00699 00700 ssnorm_time = MPI_Wtime(); 00701 00702 ssnorm_time -= gnorm_time; 00703 gnorm_time -= interp_time; 00704 interp_time -= pointloc_time; 00705 pointloc_time -= instant_time; 00706 instant_time -= start_time; 00707 00708 // set field values as tag on target vertices 00709 if (specSou) 00710 { 00711 // create a new tag for the values on the target 00712 Tag tag; 00713 std::string newtag = interpTag +"_TAR"; 00714 result = mbImpl->tag_get_handle(newtag.c_str(), 1, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT|MB_TAG_DENSE); 00715 result = mbImpl->tag_set_data(tag, targ_verts, &field[0]); 00716 PRINT_LAST_ERROR; 00717 } 00718 else 00719 { 00720 if (!specTar) 00721 { 00722 // use original tag 00723 Tag tag; 00724 result = mbImpl->tag_get_handle(interpTag.c_str(), 1, MB_TYPE_DOUBLE, tag); PRINT_LAST_ERROR; 00725 result = mbImpl->tag_set_data(tag, targ_verts, &field[0]); 00726 PRINT_LAST_ERROR; 00727 } 00728 else 00729 { 00730 // we have the field values computed at each GL points, on each element 00731 // in the target mesh 00732 // we need to create a new tag, on elements, of size _ntot, to hold 00733 // all those values. 00734 // so it turns out we need ntot. maybe we can compute it from the 00735 // number of values computed, divided by number of elements 00736 int ntot = numPointsOfInterest/targ_elems.size(); 00737 Tag tag; 00738 std::string newtag = interpTag +"_TAR"; 00739 result = mbImpl->tag_get_handle(newtag.c_str(), ntot, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT|MB_TAG_DENSE); 00740 result = mbImpl->tag_set_data(tag, targ_elems, &field[0]); 00741 PRINT_LAST_ERROR; 00742 00743 } 00744 } 00745 // done 00746 return MB_SUCCESS; 00747 }