moab
mbcoupler_test.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines