MOAB: Mesh Oriented datABase  (version 5.1.1)
Coupler.cpp
Go to the documentation of this file.
00001 #include "Coupler.hpp"
00002 #include "moab/ParallelComm.hpp"
00003 #include "moab/AdaptiveKDTree.hpp"
00004 #include "ElemUtil.hpp"
00005 #include "moab/CN.hpp"
00006 #include "moab/gs.hpp"
00007 #include "moab/TupleList.hpp"
00008 #include "moab/Error.hpp"
00009 
00010 /* C++ includes */
00011 #include <cstdio>
00012 #include <cassert>
00013 #include <iostream>
00014 #include <algorithm>
00015 #include <sstream>
00016 
00017 #define ERROR(a) {if (MB_SUCCESS != err) std::cerr << a << std::endl;}
00018 #define ERRORR(a,b) {if (MB_SUCCESS != b) {std::cerr << a << std::endl; return b;}}
00019 #define ERRORMPI(a,b) {if (MPI_SUCCESS != b) {std::cerr << a << std::endl; return MB_FAILURE;}}
00020 
00021 #define MASTER_PROC 0
00022 
00023 namespace moab {
00024 
00025 bool debug = false;
00026 int pack_tuples(TupleList *tl, void **ptr);
00027 void unpack_tuples(void *ptr, TupleList **tlp);
00028 
00029 Coupler::Coupler(Interface *impl,
00030                  ParallelComm *pc,
00031                  Range &local_elems,
00032                  int coupler_id,
00033                  bool init_tree,
00034                  int max_ent_dim)
00035   : mbImpl(impl), myPc(pc), myId(coupler_id), numIts(3), max_dim(max_ent_dim), _ntot(0), spherical(false)
00036 {
00037   assert(NULL != impl && (pc || !local_elems.empty()));
00038 
00039   // Keep track of the local points, at least for now
00040   myRange = local_elems;
00041   myTree = NULL;
00042 
00043   // Now initialize the tree
00044   if (init_tree)
00045     initialize_tree();
00046 
00047   // Initialize tuple lists to indicate not initialized
00048   mappedPts = NULL;
00049   targetPts = NULL;
00050   _spectralSource = _spectralTarget = NULL;
00051 }
00052 
00053   /* Destructor
00054    */
00055 Coupler::~Coupler()
00056 {
00057   // This will clear the cache
00058   delete (moab::Element::SpectralHex*)_spectralSource;
00059   delete (moab::Element::SpectralHex*)_spectralTarget;
00060   delete myTree;
00061   delete targetPts;
00062   delete mappedPts;
00063 }
00064 
00065 ErrorCode Coupler::initialize_tree()
00066 {
00067   Range local_ents;
00068 
00069   // Get entities on the local part
00070   ErrorCode result = MB_SUCCESS;
00071   if (myPc)
00072   {
00073     result = myPc->get_part_entities(local_ents, max_dim);
00074     if (local_ents.empty())
00075     {
00076       max_dim--;
00077       result = myPc->get_part_entities(local_ents, max_dim);// go one dimension lower
00078       // right now, this is used for spherical meshes only
00079     }
00080   }
00081   else local_ents = myRange;
00082   if (MB_SUCCESS != result || local_ents.empty()) {
00083     std::cout << "Problems getting source entities" << std::endl;
00084     return result;
00085   }
00086 
00087   // Build the tree for local processor
00088   int max_per_leaf = 6;
00089   for (int i = 0; i < numIts; i++) {
00090     std::ostringstream str;
00091     str << "PLANE_SET=0;"
00092         << "MAX_PER_LEAF=" << max_per_leaf << ";";
00093     if (spherical && !local_ents.empty())
00094     {
00095       // get coordinates of one vertex, and use for radius computation
00096       EntityHandle elem= local_ents[0];
00097       const EntityHandle * conn;
00098       int numn=0;
00099       mbImpl->get_connectivity(elem, conn, numn);
00100       CartVect pos0;
00101       mbImpl -> get_coords(conn, 1, &(pos0[0]) );
00102       double radius = pos0.length();
00103       str<<"SPHERICAL=true;RADIUS="<<radius<<";";
00104     }
00105     FileOptions opts(str.str().c_str());
00106     myTree = new AdaptiveKDTree(mbImpl);
00107     result = myTree->build_tree(local_ents, &localRoot, &opts);
00108     if (MB_SUCCESS != result) {
00109       std::cout << "Problems building tree";
00110       if (numIts != i) {
00111         delete myTree;
00112         max_per_leaf *= 2;
00113         std::cout << "; increasing elements/leaf to "
00114                   << max_per_leaf << std::endl;
00115       }
00116       else {
00117         std::cout << "; exiting" << std::endl;
00118         return result;
00119       }
00120     }
00121     else
00122       break; // Get out of tree building
00123   }
00124 
00125   // Get the bounding box for local tree
00126   if (myPc)
00127     allBoxes.resize(6*myPc->proc_config().proc_size());
00128   else
00129     allBoxes.resize(6);
00130 
00131   unsigned int my_rank = (myPc ? myPc->proc_config().proc_rank() : 0);
00132   BoundBox box;
00133   result = myTree->get_bounding_box(box, &localRoot);
00134   if (MB_SUCCESS != result)
00135     return result;
00136   box.bMin.get(&allBoxes[6*my_rank]);
00137   box.bMax.get(&allBoxes[6*my_rank + 3]);
00138 
00139   // Now communicate to get all boxes
00140   if (myPc) {
00141     int mpi_err;
00142 #if (MPI_VERSION >= 2)
00143     // Use "in place" option
00144     mpi_err = MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL,
00145                             &allBoxes[0], 6, MPI_DOUBLE,
00146                             myPc->proc_config().proc_comm());
00147 #else
00148     {
00149       std::vector<double> allBoxes_tmp(6*myPc->proc_config().proc_size());
00150       mpi_err = MPI_Allgather(&allBoxes[6*my_rank], 6, MPI_DOUBLE,
00151                               &allBoxes_tmp[0], 6, MPI_DOUBLE,
00152                               myPc->proc_config().proc_comm());
00153       allBoxes = allBoxes_tmp;
00154     }
00155 #endif
00156     if (MPI_SUCCESS != mpi_err)
00157       return MB_FAILURE;
00158   }
00159 
00160   /*  std::ostringstream blah;
00161   for(int i = 0; i < allBoxes.size(); i++)
00162   blah << allBoxes[i] << " ";
00163   std::cout << blah.str() << "\n";*/
00164 
00165 #ifdef VERBOSE
00166   double min[3] = {0, 0, 0}, max[3] = {0, 0, 0};
00167   unsigned int dep;
00168   myTree->get_info(localRoot, min, max, dep);
00169   std::cout << "Proc " << my_rank << ": box min/max, tree depth = ("
00170             << min[0] << "," << min[1] << "," << min[2] << "), ("
00171             << max[0] << "," << max[1] << "," << max[2] << "), "
00172             << dep << std::endl;
00173 #endif
00174 
00175   return result;
00176 }
00177 
00178 ErrorCode Coupler::initialize_spectral_elements(EntityHandle rootSource, EntityHandle rootTarget,
00179                                                 bool &specSou, bool &specTar)
00180 {
00181   /*void * _spectralSource;
00182     void * _spectralTarget;*/
00183 
00184   moab::Range spectral_sets;
00185   moab::Tag  sem_tag;
00186   int sem_dims[3];
00187   ErrorCode rval = mbImpl->tag_get_handle("SEM_DIMS", 3, moab::MB_TYPE_INTEGER, sem_tag);
00188   if (moab::MB_SUCCESS != rval) {
00189     std::cout << "Can't find tag, no spectral set\n";
00190     return MB_SUCCESS; // Nothing to do, no spectral elements
00191   }
00192   rval = mbImpl->get_entities_by_type_and_tag(rootSource, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets);
00193   if (moab::MB_SUCCESS != rval || spectral_sets.empty())
00194     std::cout << "Can't get sem set on source\n";
00195   else {
00196     moab::EntityHandle firstSemSet=spectral_sets[0];
00197     rval = mbImpl->tag_get_data(sem_tag, &firstSemSet, 1, (void*)sem_dims);
00198     if (moab::MB_SUCCESS != rval)
00199       return MB_FAILURE;
00200 
00201     if (sem_dims[0]!=sem_dims[1] || sem_dims[0] != sem_dims[2]) {
00202       std::cout << " dimensions are different. bail out\n";
00203       return MB_FAILURE;
00204     }
00205     // Repeat for target sets
00206     spectral_sets.empty();
00207     // Now initialize a source spectral element !
00208     _spectralSource = new moab::Element::SpectralHex(sem_dims[0]);
00209     specSou = true;
00210   }
00211   // Repeat for target source
00212   rval = mbImpl->get_entities_by_type_and_tag(rootTarget, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets);
00213   if (moab::MB_SUCCESS != rval || spectral_sets.empty())
00214     std::cout << "Can't get sem set on target\n";
00215   else {
00216     moab::EntityHandle firstSemSet=spectral_sets[0];
00217     rval = mbImpl->tag_get_data(sem_tag, &firstSemSet, 1, (void*)sem_dims);
00218     if (moab::MB_SUCCESS != rval)
00219       return MB_FAILURE;
00220 
00221     if (sem_dims[0]!=sem_dims[1] || sem_dims[0] != sem_dims[2]) {
00222       std::cout << " dimensions are different. bail out\n";
00223       return MB_FAILURE;
00224     }
00225     // Repeat for target sets
00226     spectral_sets.empty();
00227     // Now initialize a target spectral element !
00228     _spectralTarget = new moab::Element::SpectralHex(sem_dims[0]);
00229     specTar = true;
00230   }
00231 
00232   _ntot = sem_dims[0]*sem_dims[1]*sem_dims[2];
00233   rval = mbImpl->tag_get_handle("SEM_X", _ntot, moab::MB_TYPE_DOUBLE, _xm1Tag);
00234   if (moab::MB_SUCCESS != rval) {
00235      std::cout << "Can't get xm1tag \n";
00236      return MB_FAILURE;
00237   }
00238   rval = mbImpl->tag_get_handle("SEM_Y", _ntot, moab::MB_TYPE_DOUBLE, _ym1Tag);
00239   if (moab::MB_SUCCESS != rval) {
00240      std::cout << "Can't get ym1tag \n";
00241      return MB_FAILURE;
00242   }
00243   rval = mbImpl->tag_get_handle("SEM_Z", _ntot, moab::MB_TYPE_DOUBLE, _zm1Tag);
00244   if (moab::MB_SUCCESS != rval) {
00245      std::cout << "Can't get zm1tag \n";
00246      return MB_FAILURE;
00247   }
00248 
00249   return MB_SUCCESS;
00250 }
00251 
00252 ErrorCode Coupler::locate_points(Range &targ_ents,
00253                                  double rel_eps,
00254                                  double abs_eps,
00255                                  TupleList *tl,
00256                                  bool store_local)
00257 {
00258   // Get locations
00259   std::vector<double> locs(3*targ_ents.size());
00260   Range verts = targ_ents.subset_by_type(MBVERTEX);
00261   ErrorCode rval = mbImpl->get_coords(verts, &locs[0]);
00262   if (MB_SUCCESS != rval)
00263     return rval;
00264   // Now get other ents; reuse verts
00265   unsigned int num_verts = verts.size();
00266   verts = subtract(targ_ents, verts);
00267   // Compute centroids
00268   std::vector<EntityHandle> dum_conn(CN::MAX_NODES_PER_ELEMENT);
00269   std::vector<double> dum_pos(CN::MAX_NODES_PER_ELEMENT);
00270   const EntityHandle *conn;
00271   int num_conn;
00272   double *coords = &locs[num_verts];
00273   // Do this here instead of a function to allow reuse of dum_pos and dum_conn
00274   for (Range::const_iterator rit = verts.begin(); rit != verts.end(); ++rit) {
00275     rval = mbImpl->get_connectivity(*rit, conn, num_conn, false, &dum_conn);
00276     if (MB_SUCCESS != rval)
00277       return rval;
00278     rval = mbImpl->get_coords(conn, num_conn, &dum_pos[0]);
00279     if (MB_SUCCESS != rval)
00280       return rval;
00281     coords[0] = coords[1] = coords[2] = 0.0;
00282     for (int i = 0; i < num_conn; i++) {
00283       coords[0] += dum_pos[3*i];
00284       coords[1] += dum_pos[3*i + 1];
00285       coords[2] += dum_pos[3*i + 2];
00286     }
00287     coords[0] /= num_conn;
00288     coords[1] /= num_conn;
00289     coords[2] /= num_conn;
00290     coords += 3;
00291   }
00292 
00293   if (store_local)
00294     targetEnts = targ_ents;
00295 
00296   return locate_points(&locs[0], targ_ents.size(), rel_eps, abs_eps, tl, store_local);
00297 }
00298 
00299 ErrorCode Coupler::locate_points(double *xyz, unsigned int num_points,
00300                                  double rel_eps,
00301                                  double abs_eps,
00302                                  TupleList *tl,
00303                                  bool store_local)
00304 {
00305   assert(tl || store_local);
00306 
00307   // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing pts to be located
00308   // source_pts: TL(from_proc, tgt_index, src_index): results of source mesh proc point location, ready to send
00309   //             back to tgt procs; src_index of -1 indicates point not located (arguably not useful...)
00310   TupleList target_pts;
00311   target_pts.initialize(2, 0, 0, 3, num_points);
00312   target_pts.enableWriteAccess();
00313 
00314   TupleList source_pts;
00315   mappedPts = new TupleList(0, 0, 1, 3, target_pts.get_max());
00316   mappedPts->enableWriteAccess();
00317 
00318   source_pts.initialize(3, 0, 0, 0, target_pts.get_max());
00319   source_pts.enableWriteAccess();
00320 
00321   mappedPts->set_n(0);
00322   source_pts.set_n(0);
00323   ErrorCode result;
00324 
00325   unsigned int my_rank = (myPc ? myPc->proc_config().proc_rank() : 0);
00326   bool point_located;
00327 
00328   // For each point, find box(es) containing the point,
00329   // appending results to tuple_list;
00330   // keep local points separately, in local_pts, which has pairs
00331   // of <local_index, mapped_index>, where mapped_index is the index
00332   // into the mappedPts tuple list
00333   for (unsigned int i = 0; i < 3*num_points; i += 3) {
00334 
00335     std::vector<int>  procs_to_send_to;
00336     for (unsigned int j = 0; j < (myPc ? myPc->proc_config().proc_size() : 0); j++) {
00337       // Test if point is in proc's box
00338       if ((allBoxes[6*j] <= xyz[i] + abs_eps) && (xyz[i] <= allBoxes[6*j + 3] + abs_eps) &&
00339           (allBoxes[6*j + 1] <= xyz[i + 1] + abs_eps) && (xyz[i + 1] <= allBoxes[6*j + 4] + abs_eps) &&
00340           (allBoxes[6*j + 2] <= xyz[i + 2] + abs_eps) && (xyz[i + 2] <= allBoxes[6*j + 5] + abs_eps)) {
00341         // If in this proc's box, will send to proc to test further
00342         procs_to_send_to.push_back(j);
00343       }
00344     }
00345     if (procs_to_send_to.empty())
00346     {
00347 #ifdef VERBOSE
00348       std::cout << " point index " << i/3 << ": " << xyz[i] << " " << xyz[i+1] << " " << xyz[i+2] << " not found in any box\n";
00349 #endif
00350       // try to find the closest box, and put it in that box, anyway
00351       double min_dist = 1.e+20;
00352       int index = -1;
00353       for (unsigned int j = 0; j < (myPc ? myPc->proc_config().proc_size() : 0); j++)
00354       {
00355         BoundBox box( &allBoxes[6*j]); // form back the box
00356         double distance = box.distance(&xyz[i]); // will compute the distance in 3d, from the box
00357         if (distance < min_dist)
00358         {
00359           index = j;
00360           min_dist = distance;
00361         }
00362       }
00363       if (index ==-1)
00364       {
00365         // need to abort early, nothing we can do
00366         assert("cannot locate any box for some points");
00367         // need a better exit strategy
00368       }
00369 #ifdef VERBOSE
00370       std::cout << " point index " << i/3 << " added to box for proc j:" << index << "\n";
00371 #endif
00372       procs_to_send_to.push_back(index); // will send to just one proc, that has the closest box
00373     }
00374     // we finally decided to populate the tuple list for a list of processors
00375     for (size_t k=0; k<procs_to_send_to.size(); k++)
00376     {
00377       unsigned int j = procs_to_send_to[k];
00378       // Check size of tuple list, grow if we're at max
00379       if (target_pts.get_n() == target_pts.get_max())
00380              target_pts.resize(std::max(10.0, 1.5*target_pts.get_max()));
00381 
00382       target_pts.vi_wr[2*target_pts.get_n()] = j;
00383       target_pts.vi_wr[2*target_pts.get_n() + 1] = i / 3;
00384 
00385       target_pts.vr_wr[3*target_pts.get_n()] = xyz[i];
00386       target_pts.vr_wr[3*target_pts.get_n() + 1] = xyz[i + 1];
00387       target_pts.vr_wr[3*target_pts.get_n() + 2] = xyz[i + 2];
00388       target_pts.inc_n();
00389     }
00390 
00391   } // end for (unsigned int i = 0; ..
00392 
00393   int num_to_me = 0;
00394   for (unsigned int i = 0; i < target_pts.get_n(); i++)
00395     if (target_pts.vi_rd[2*i] == (int)my_rank)
00396       num_to_me++;
00397 #ifdef VERBOSE
00398   printf("rank: %u local points: %u, nb sent target pts: %u mappedPts: %u num to me: %d \n",
00399          my_rank, num_points, target_pts.get_n(), mappedPts->get_n(), num_to_me);
00400 #endif
00401   // Perform scatter/gather, to gather points to source mesh procs
00402   if (myPc) {
00403     (myPc->proc_config().crystal_router())->gs_transfer(1, target_pts, 0);
00404 
00405     num_to_me = 0;
00406     for (unsigned int i = 0; i < target_pts.get_n(); i++) {
00407       if (target_pts.vi_rd[2*i] == (int)my_rank)
00408         num_to_me++;
00409     }
00410 #ifdef VERBOSE
00411     printf("rank: %u after first gs nb received_pts: %u; num_from_me = %d\n",
00412            my_rank, target_pts.get_n(), num_to_me);
00413 #endif
00414     // After scatter/gather:
00415     // target_pts.set_n(# points local proc has to map);
00416     // target_pts.vi_wr[2*i] = proc sending point i
00417     // target_pts.vi_wr[2*i + 1] = index of point i on sending proc
00418     // target_pts.vr_wr[3*i..3*i + 2] = xyz of point i
00419     //
00420     // Mapping builds the tuple list:
00421     // source_pts.set_n(target_pts.get_n())
00422     // source_pts.vi_wr[3*i] = target_pts.vi_wr[2*i] = sending proc
00423     // source_pts.vi_wr[3*i + 1] = index of point i on sending proc
00424     // source_pts.vi_wr[3*i + 2] = index of mapped point (-1 if not mapped)
00425     //
00426     // Also, mapping builds local tuple_list mappedPts:
00427     // mappedPts->set_n( # mapped points );
00428     // mappedPts->vul_wr[i] = local handle of mapped entity
00429     // mappedPts->vr_wr[3*i..3*i + 2] = natural coordinates in mapped entity
00430 
00431     // Test target points against my elements
00432     for (unsigned i = 0; i < target_pts.get_n(); i++) {
00433       result = test_local_box(target_pts.vr_wr + 3*i,
00434                               target_pts.vi_rd[2*i], target_pts.vi_rd[2*i + 1], i,
00435                               point_located, rel_eps, abs_eps, &source_pts);
00436       if (MB_SUCCESS != result)
00437         return result;
00438     }
00439 
00440     // No longer need target_pts
00441     target_pts.reset();
00442 #ifdef VERBOSE
00443     printf("rank: %u nb sent source pts: %u, mappedPts now: %u\n",
00444            my_rank, source_pts.get_n(),  mappedPts->get_n());
00445 #endif
00446       // Send target points back to target procs
00447     (myPc->proc_config().crystal_router())->gs_transfer(1, source_pts, 0);
00448 
00449 #ifdef VERBOSE
00450     printf("rank: %u nb received source pts: %u\n",
00451            my_rank, source_pts.get_n());
00452 #endif
00453   }
00454 
00455   // Store proc/index tuples in targetPts, and/or pass back to application;
00456   // the tuple this gets stored to looks like:
00457   // tl.set_n(# mapped points);
00458   // tl.vi_wr[3*i] = remote proc mapping point
00459   // tl.vi_wr[3*i + 1] = local index of mapped point
00460   // tl.vi_wr[3*i + 2] = remote index of mapped point
00461   //
00462   // Local index is mapped into either myRange, holding the handles of
00463   // local mapped entities, or myXyz, holding locations of mapped pts
00464 
00465   // Store information about located points
00466   TupleList *tl_tmp;
00467   if (!store_local)
00468     tl_tmp = tl;
00469   else {
00470     targetPts = new TupleList();
00471     tl_tmp = targetPts;
00472   }
00473 
00474   tl_tmp->initialize(3, 0, 0, 0, num_points);
00475   tl_tmp->set_n(num_points); // Automatically sets tl to write_enabled
00476   // Initialize so we know afterwards how many pts weren't located
00477   std::fill(tl_tmp->vi_wr, tl_tmp->vi_wr + 3*num_points, -1);
00478 
00479   unsigned int local_pts = 0;
00480   for (unsigned int i = 0; i < source_pts.get_n(); i++) {
00481     if (-1 != source_pts.vi_rd[3*i + 2]) { // Why bother sending message saying "i don't have the point" if it gets discarded?
00482       int tgt_index = 3*source_pts.vi_rd[3*i + 1];
00483       // Prefer always entities that are local, from the source_pts
00484       // if a local entity was already found to contain the target point, skip
00485       // tl_tmp->vi_wr[tgt_index] is -1 initially, but it could already be set with
00486       // a remote processor
00487       if (tl_tmp->vi_wr[tgt_index] != (int)my_rank) {
00488         tl_tmp->vi_wr[tgt_index]     = source_pts.vi_rd[3*i];
00489         tl_tmp->vi_wr[tgt_index + 1] = source_pts.vi_rd[3*i + 1];
00490         tl_tmp->vi_wr[tgt_index + 2] = source_pts.vi_rd[3*i + 2];
00491       }
00492     }
00493   }
00494 
00495   // Count missing points
00496   unsigned int missing_pts = 0;
00497   for (unsigned int i = 0; i < num_points; i++) {
00498     if (tl_tmp->vi_rd[3*i + 1] == -1) {
00499       missing_pts++;
00500 #ifdef VERBOSE
00501       printf("missing point at index i:  %d -> %15.10f %15.10f %15.10f\n", i, xyz[3*i], xyz[3*i + 1], xyz[3*i + 2]);
00502 #endif
00503     }
00504     else if (tl_tmp->vi_rd[3*i] == (int)my_rank)
00505       local_pts++;
00506   }
00507 #ifdef VERBOSE
00508   printf("rank: %u point location: wanted %u got %u locally, %u remote, missing %u\n",
00509          my_rank, num_points, local_pts, num_points - missing_pts - local_pts, missing_pts);
00510 #endif
00511   assert(0 == missing_pts); // Will likely break on curved geometries
00512 
00513   // No longer need source_pts
00514   source_pts.reset();
00515 
00516   // Copy into tl if passed in and storing locally
00517   if (tl && store_local) {
00518     tl->initialize(3, 0, 0, 0, num_points);
00519     tl->enableWriteAccess();
00520     memcpy(tl->vi_wr, tl_tmp->vi_rd, 3 * tl_tmp->get_n() * sizeof(int));
00521     tl->set_n(tl_tmp->get_n());
00522     tl->disableWriteAccess();
00523   }
00524 
00525   tl_tmp->disableWriteAccess();
00526 
00527   // Done
00528   return MB_SUCCESS;
00529 }
00530 
00531 ErrorCode Coupler::test_local_box(double *xyz,
00532                                   int from_proc, int remote_index, int /*index*/,
00533                                   bool &point_located,
00534                                   double rel_eps, double abs_eps,
00535                                   TupleList *tl)
00536 {
00537   std::vector<EntityHandle> entities;
00538   std::vector<CartVect> nat_coords;
00539   bool canWrite = false;
00540   if (tl) {
00541     canWrite = tl->get_writeEnabled();
00542     if (!canWrite) {
00543       tl->enableWriteAccess();
00544       canWrite = true;
00545     }
00546   }
00547 
00548   if (rel_eps && !abs_eps) {
00549     // Relative epsilon given, translate to absolute epsilon using box dimensions
00550     BoundBox box;
00551     myTree->get_bounding_box(box, &localRoot);
00552     abs_eps = rel_eps * box.diagonal_length();
00553   }
00554 
00555   ErrorCode result = nat_param(xyz, entities, nat_coords, abs_eps);
00556   if (MB_SUCCESS != result)
00557     return result;
00558 
00559   // If we didn't find any ents and we're looking locally, nothing more to do
00560   if (entities.empty()) {
00561     if (tl->get_n() == tl->get_max())
00562       tl->resize(std::max(10.0, 1.5 * tl->get_max()));
00563 
00564     tl->vi_wr[3 * tl->get_n()] = from_proc;
00565     tl->vi_wr[3 * tl->get_n() + 1] = remote_index;
00566     tl->vi_wr[3 * tl->get_n() + 2] = -1;
00567     tl->inc_n();
00568 
00569     point_located = false;
00570     return MB_SUCCESS;
00571   }
00572 
00573   // Grow if we know we'll exceed size
00574   if (mappedPts->get_n() + entities.size() >= mappedPts->get_max())
00575     mappedPts->resize(std::max(10.0, 1.5 * mappedPts->get_max()));
00576 
00577   std::vector<EntityHandle>::iterator eit = entities.begin();
00578   std::vector<CartVect>::iterator ncit = nat_coords.begin();
00579 
00580   mappedPts->enableWriteAccess();
00581   for (; eit != entities.end(); ++eit, ++ncit) {
00582     // Store in tuple mappedPts
00583     mappedPts->vr_wr[3*mappedPts->get_n()] = (*ncit)[0];
00584     mappedPts->vr_wr[3*mappedPts->get_n() + 1] = (*ncit)[1];
00585     mappedPts->vr_wr[3*mappedPts->get_n() + 2] = (*ncit)[2];
00586     mappedPts->vul_wr[mappedPts->get_n()] = *eit;
00587     mappedPts->inc_n();
00588 
00589     // Also store local point, mapped point indices
00590     if (tl->get_n() == tl->get_max())
00591       tl->resize(std::max(10.0, 1.5*tl->get_max()));
00592 
00593     // Store in tuple source_pts
00594     tl->vi_wr[3*tl->get_n()] = from_proc;
00595     tl->vi_wr[3*tl->get_n() + 1] = remote_index;
00596     tl->vi_wr[3*tl->get_n() + 2] = mappedPts->get_n()-1;
00597     tl->inc_n();
00598   }
00599 
00600   point_located = true;
00601 
00602   if (tl && !canWrite)
00603     tl->disableWriteAccess();
00604 
00605   return MB_SUCCESS;
00606 }
00607 
00608 ErrorCode Coupler::interpolate( Coupler::Method method,
00609                                 const std::string &interp_tag,
00610                                 double *interp_vals,
00611                                 TupleList *tl,
00612                                 bool normalize)
00613 {
00614   Tag tag;
00615   ErrorCode result ;
00616   if (_spectralSource) {
00617     result = mbImpl->tag_get_handle(interp_tag.c_str(), _ntot, MB_TYPE_DOUBLE, tag);MB_CHK_SET_ERR(result, "Failed to get handle for interpolation tag \"" << interp_tag << "\"");
00618   }
00619   else {
00620     result = mbImpl->tag_get_handle(interp_tag.c_str(), 1, MB_TYPE_DOUBLE, tag);MB_CHK_SET_ERR(result, "Failed to get handle for interpolation tag \"" << interp_tag << "\"");
00621   }
00622 
00623   return interpolate(method, tag, interp_vals, tl, normalize);
00624 }
00625 
00626 ErrorCode Coupler::interpolate(Coupler::Method *methods,
00627                                Tag *tags,
00628                                int *points_per_method,
00629                                int num_methods,
00630                                double *interp_vals,
00631                                TupleList *tl,
00632                                bool /* normalize */)
00633 {
00634   //if (!((LINEAR_FE == method) || (CONSTANT == method)))
00635   // return MB_FAILURE;
00636 
00637   // remote pts first
00638   TupleList *tl_tmp = (tl ? tl : targetPts);
00639 
00640   ErrorCode result = MB_SUCCESS;
00641 
00642   unsigned int pts_total = 0;
00643   for (int i = 0; i < num_methods; i++)
00644     pts_total += points_per_method[i];
00645 
00646   // If tl was passed in non-NULL, just have those points, otherwise have targetPts plus
00647   // locally mapped pts
00648   if (pts_total != tl_tmp->get_n())
00649     return MB_FAILURE;
00650 
00651   TupleList tinterp;
00652   tinterp.initialize(5, 0, 0, 1, tl_tmp->get_n());
00653   int t = 0;
00654   tinterp.enableWriteAccess();
00655   for (int i = 0; i < num_methods; i++) {
00656     for (int j = 0; j < points_per_method[i]; j++) {
00657       tinterp.vi_wr[5*t] = tl_tmp->vi_rd[3*t];
00658       tinterp.vi_wr[5*t + 1] = tl_tmp->vi_rd[3*t + 1];
00659       tinterp.vi_wr[5*t + 2] = tl_tmp->vi_rd[3*t + 2];
00660       tinterp.vi_wr[5*t + 3] = methods[i];
00661       tinterp.vi_wr[5*t + 4] = i;
00662       tinterp.vr_wr[t] = 0.0;
00663       tinterp.inc_n();
00664       t++;
00665     }
00666   }
00667 
00668   // Scatter/gather interpolation points
00669   if (myPc) {
00670     (myPc->proc_config().crystal_router())->gs_transfer(1, tinterp, 0);
00671 
00672     // Perform interpolation on local source mesh; put results into
00673     // tinterp.vr_wr[i]
00674 
00675     mappedPts->enableWriteAccess();
00676     for (unsigned int i = 0; i < tinterp.get_n(); i++) {
00677       int mindex = tinterp.vi_rd[5*i + 2];
00678       Method method = (Method)tinterp.vi_rd[5*i + 3];
00679       Tag tag = tags[tinterp.vi_rd[5*i + 4]];
00680 
00681       result = MB_FAILURE;
00682       if (LINEAR_FE == method || QUADRATIC_FE == method || SPHERICAL==method) {
00683         result = interp_field(mappedPts->vul_rd[mindex],
00684                               CartVect(mappedPts->vr_wr + 3*mindex),
00685                               tag, tinterp.vr_wr[i]);
00686       }
00687       else if (CONSTANT == method) {
00688         result = constant_interp(mappedPts->vul_rd[mindex],
00689                                  tag, tinterp.vr_wr[i]);
00690       }
00691 
00692       if (MB_SUCCESS != result)
00693         return result;
00694     }
00695 
00696     // Scatter/gather interpolation data
00697     myPc->proc_config().crystal_router()->gs_transfer(1, tinterp, 0);
00698   }
00699 
00700   // Copy the interpolated field as a unit
00701   for (unsigned int i = 0; i < tinterp.get_n(); i++)
00702     interp_vals[tinterp.vi_rd[5*i + 1]] = tinterp.vr_rd[i];
00703 
00704   // Done
00705   return MB_SUCCESS;
00706 }
00707 
00708 ErrorCode Coupler::nat_param(double xyz[3],
00709                              std::vector<EntityHandle> &entities,
00710                              std::vector<CartVect> &nat_coords,
00711                              double epsilon)
00712 {
00713   if (!myTree)
00714     return MB_FAILURE;
00715 
00716   AdaptiveKDTreeIter treeiter;
00717   ErrorCode result = myTree->get_tree_iterator(localRoot, treeiter);
00718   if (MB_SUCCESS != result) {
00719     std::cout << "Problems getting iterator" << std::endl;
00720     return result;
00721   }
00722 
00723   EntityHandle closest_leaf;
00724   if (epsilon) {
00725     std::vector<double> dists;
00726     std::vector<EntityHandle> leaves;
00727 
00728     // Two tolerances
00729     result = myTree->distance_search(xyz, epsilon, leaves,
00730                                      /*iter_tol*/ epsilon,
00731                                      /*inside_tol*/ 10*epsilon,
00732                                      &dists, NULL, &localRoot);
00733     if (leaves.empty())
00734       // Not found returns success here, with empty list, just like case with no epsilon
00735       return MB_SUCCESS;
00736     // Get closest leaf
00737     double min_dist = *dists.begin();
00738     closest_leaf = *leaves.begin();
00739     std::vector<EntityHandle>::iterator vit = leaves.begin() + 1;
00740     std::vector<double>::iterator dit = dists.begin() + 1;
00741     for (; vit != leaves.end() && min_dist; ++vit, ++dit) {
00742       if (*dit < min_dist) {
00743         min_dist = *dit;
00744         closest_leaf = *vit;
00745       }
00746     }
00747   }
00748   else {
00749     result = myTree->point_search(xyz, treeiter, 1.0e-10, 1.0e-6, NULL, &localRoot);
00750     if (MB_ENTITY_NOT_FOUND == result) // Point is outside of myTree's bounding box
00751       return MB_SUCCESS;
00752     else if (MB_SUCCESS != result) {
00753       std::cout << "Problems getting leaf \n";
00754       return result;
00755     }
00756     closest_leaf = treeiter.handle();
00757   }
00758 
00759   // Find natural coordinates of point in element(s) in that leaf
00760   CartVect tmp_nat_coords;
00761   Range range_leaf;
00762   result = mbImpl->get_entities_by_dimension(closest_leaf, max_dim, range_leaf, false);
00763   if (MB_SUCCESS != result)
00764     std::cout << "Problem getting leaf in a range" << std::endl;
00765 
00766   // Loop over the range_leaf
00767   for (Range::iterator iter = range_leaf.begin(); iter != range_leaf.end(); ++iter) {
00768     // Test to find out in which entity the point is
00769     // Get the EntityType and create the appropriate Element::Map subtype
00770     // If spectral, do not need coordinates, just the GL points
00771     EntityType etype = mbImpl->type_from_handle(*iter);
00772     if (NULL != this->_spectralSource && MBHEX == etype) {
00773       EntityHandle eh = *iter;
00774       const double * xval;
00775       const double * yval;
00776       const double * zval;
00777       ErrorCode rval = mbImpl->tag_get_by_ptr(_xm1Tag, &eh, 1, (const void**)&xval);
00778       if (moab::MB_SUCCESS != rval) {
00779         std::cout << "Can't get xm1 values \n";
00780         return MB_FAILURE;
00781       }
00782       rval = mbImpl->tag_get_by_ptr(_ym1Tag, &eh, 1, (const void**)&yval);
00783       if (moab::MB_SUCCESS != rval) {
00784         std::cout << "Can't get ym1 values \n";
00785         return MB_FAILURE;
00786       }
00787       rval = mbImpl->tag_get_by_ptr(_zm1Tag, &eh, 1, (const void**)&zval);
00788       if (moab::MB_SUCCESS != rval) {
00789         std::cout << "Can't get zm1 values \n";
00790         return MB_FAILURE;
00791       }
00792       Element::SpectralHex* spcHex = (Element::SpectralHex*)_spectralSource;
00793 
00794       spcHex->set_gl_points((double*)xval, (double*)yval, (double*)zval);
00795       try {
00796         tmp_nat_coords = spcHex->ievaluate(CartVect(xyz), epsilon ); // introduce
00797         bool inside = spcHex->inside_nat_space(CartVect(tmp_nat_coords), epsilon);
00798         if (!inside) {
00799 #ifdef VERBOSE
00800           std::cout << "point " << xyz[0] << " " << xyz[1] << " " << xyz[2] <<
00801               " is not converging inside hex " << mbImpl->id_from_handle(eh) << "\n";
00802 #endif
00803           continue; // It is possible that the point is outside, so it will not converge
00804         }
00805       }
00806       catch (Element::Map::EvaluationError&) {
00807         continue;
00808       }
00809 
00810 
00811     }
00812     else {
00813       const EntityHandle *connect;
00814       int num_connect;
00815 
00816       // Get connectivity
00817       result = mbImpl->get_connectivity(*iter, connect, num_connect, true);
00818       if (MB_SUCCESS != result)
00819         return result;
00820 
00821       // Get coordinates of the vertices
00822       std::vector<CartVect> coords_vert(num_connect);
00823       result = mbImpl->get_coords(connect, num_connect, &(coords_vert[0][0]));
00824       if (MB_SUCCESS != result) {
00825         std::cout << "Problems getting coordinates of vertices\n";
00826         return result;
00827       }
00828       CartVect pos(xyz);
00829       if (MBHEX == etype) {
00830         if (8 == num_connect) {
00831           Element::LinearHex hexmap(coords_vert);
00832           if (!hexmap.inside_box(pos, epsilon))
00833             continue;
00834           try {
00835             tmp_nat_coords = hexmap.ievaluate(pos, epsilon);
00836             bool inside = hexmap.inside_nat_space(tmp_nat_coords, epsilon);
00837             if (!inside)
00838               continue;
00839           }
00840           catch (Element::Map::EvaluationError&) {
00841             continue;
00842           }
00843         }
00844         else if (27 == num_connect) {
00845           Element::QuadraticHex hexmap(coords_vert);
00846          if (!hexmap.inside_box(pos, epsilon))
00847            continue;
00848          try {
00849            tmp_nat_coords = hexmap.ievaluate(pos, epsilon);
00850            bool inside = hexmap.inside_nat_space(tmp_nat_coords, epsilon);
00851            if (!inside)
00852              continue;
00853          }
00854          catch (Element::Map::EvaluationError&) {
00855            continue;
00856          }
00857         }
00858         else // TODO this case not treated yet, no interpolation
00859           continue;
00860       }
00861       else if (MBTET == etype) {
00862         Element::LinearTet tetmap(coords_vert);
00863         // This is just a linear solve; unless degenerate, will not except
00864         tmp_nat_coords = tetmap.ievaluate(pos);
00865         bool inside = tetmap.inside_nat_space(tmp_nat_coords, epsilon);
00866         if (!inside)
00867           continue;
00868       }
00869       else if (MBQUAD == etype && spherical) {
00870         Element::SphericalQuad sphermap(coords_vert);
00871         /* skip box test, because it can filter out good elements with high curvature
00872          * if (!sphermap.inside_box(pos, epsilon))
00873           continue;*/
00874         try {
00875           tmp_nat_coords = sphermap.ievaluate(pos, epsilon);
00876           bool inside = sphermap.inside_nat_space(tmp_nat_coords, epsilon);
00877           if (!inside)
00878             continue;
00879         }
00880         catch (Element::Map::EvaluationError&) {
00881           continue;
00882         }
00883 
00884       }
00885       else if (MBTRI == etype && spherical) {
00886         Element::SphericalTri sphermap(coords_vert);
00887         /* skip box test, because it can filter out good elements with high curvature
00888          * if (!sphermap.inside_box(pos, epsilon))
00889             continue;*/
00890         try {
00891           tmp_nat_coords = sphermap.ievaluate(pos, epsilon);
00892           bool inside = sphermap.inside_nat_space(tmp_nat_coords, epsilon);
00893           if (!inside)
00894             continue;
00895         }
00896         catch (Element::Map::EvaluationError&) {
00897           continue;
00898         }
00899       }
00900 
00901       else if (MBQUAD == etype) {
00902         Element::LinearQuad quadmap(coords_vert);
00903         if (!quadmap.inside_box(pos, epsilon))
00904           continue;
00905         try {
00906           tmp_nat_coords = quadmap.ievaluate(pos, epsilon);
00907           bool inside = quadmap.inside_nat_space(tmp_nat_coords, epsilon);
00908           if (!inside)
00909             continue;
00910         }
00911         catch (Element::Map::EvaluationError&) {
00912           continue;
00913         }
00914         if (!quadmap.inside_nat_space(tmp_nat_coords, epsilon))
00915           continue;
00916       }
00917       /*
00918       else if (etype == MBTRI){
00919         Element::LinearTri trimap(coords_vert);
00920         if (!trimap.inside_box( pos, epsilon))
00921           continue;
00922         try {
00923           tmp_nat_coords = trimap.ievaluate(pos, epsilon);
00924           bool inside = trimap.inside_nat_space(tmp_nat_coords, epsilon);
00925           if (!inside) continue;
00926         }
00927         catch (Element::Map::EvaluationError) {
00928           continue;
00929         }
00930         if (!trimap.inside_nat_space(tmp_nat_coords, epsilon))
00931           continue;
00932       }
00933       */
00934       else if (etype == MBEDGE){
00935         Element::LinearEdge edgemap(coords_vert);
00936         try {
00937           tmp_nat_coords = edgemap.ievaluate(CartVect(xyz), epsilon);
00938         }
00939         catch (Element::Map::EvaluationError) {
00940           continue;
00941         }
00942         if (!edgemap.inside_nat_space(tmp_nat_coords, epsilon))
00943           continue;
00944       }
00945       else {
00946         std::cout << "Entity not Hex/Tet/Quad/Tri/Edge. Please verify." << std::endl;
00947         continue;
00948       }
00949     }
00950     // If we get here then we've found the coordinates.
00951     // Save them and the entity and return success.
00952     entities.push_back(*iter);
00953     nat_coords.push_back(tmp_nat_coords);
00954     return MB_SUCCESS;
00955   }
00956 
00957   // Didn't find any elements containing the point
00958   return MB_SUCCESS;
00959 }
00960 
00961 ErrorCode Coupler::interp_field(EntityHandle elem,
00962                                 CartVect nat_coord,
00963                                 Tag tag,
00964                                 double &field)
00965 {
00966   if (_spectralSource) {
00967     // Get tag values at the GL points for some field (Tag)
00968     const double * vx;
00969     ErrorCode rval = mbImpl-> tag_get_by_ptr(tag, &elem, 1, (const void**)&vx);
00970     if (moab::MB_SUCCESS != rval) {
00971       std::cout << "Can't get field values for the tag \n";
00972       return MB_FAILURE;
00973     }
00974     Element::SpectralHex * spcHex = (Element::SpectralHex*)_spectralSource;
00975     field = spcHex->evaluate_scalar_field(nat_coord, vx);
00976   }
00977   else {
00978     double vfields[27]; // Will work for linear hex, quadratic hex or Tets
00979     moab::Element::Map *elemMap = NULL;
00980     int num_verts = 0;
00981     // Get the EntityType
00982     // Get the tag values at the vertices
00983     const EntityHandle *connect;
00984     int num_connect;
00985     ErrorCode result = mbImpl->get_connectivity(elem, connect, num_connect);
00986     if (MB_SUCCESS != result)
00987       return result;
00988     EntityType etype = mbImpl->type_from_handle(elem);
00989     if (MBHEX == etype) {
00990       if (8 == num_connect) {
00991         elemMap = new moab::Element::LinearHex();
00992         num_verts = 8;
00993       }
00994       else { /* (MBHEX == etype && 27 == num_connect) */
00995         elemMap = new moab::Element::QuadraticHex();
00996         num_verts = 27;
00997       }
00998     }
00999     else if (MBTET == etype) {
01000       elemMap = new moab::Element::LinearTet();
01001       num_verts = 4;
01002     }
01003     else if (MBQUAD == etype) {
01004       elemMap = new moab::Element::LinearQuad();
01005       num_verts = 4;
01006     }
01007     else if (MBTRI == etype) {
01008       elemMap = new moab::Element::LinearTri();
01009       num_verts = 3;
01010     }
01011     else
01012       return MB_FAILURE;
01013 
01014     result = mbImpl->tag_get_data(tag, connect, std::min(num_verts, num_connect), vfields);
01015     if (MB_SUCCESS != result) {
01016       delete elemMap;
01017       return result;
01018     }
01019 
01020     // Function for the interpolation
01021     field = 0;
01022 
01023     // Check the number of vertices
01024     assert(num_connect >= num_verts);
01025 
01026     // Calculate the field
01027     try {
01028       field = elemMap->evaluate_scalar_field(nat_coord, vfields);
01029     }
01030     catch (moab::Element::Map::EvaluationError&) {
01031       delete elemMap;
01032       return MB_FAILURE;
01033     }
01034 
01035     delete elemMap;
01036   }
01037 
01038   return MB_SUCCESS;
01039 }
01040 
01041 // Simplest "interpolation" for element-based source fields. Set the value of the field
01042 // at the target point to that of the field in the source element it lies in.
01043 ErrorCode Coupler::constant_interp(EntityHandle elem,
01044                                    Tag tag,
01045                                    double &field)
01046 {
01047   double tempField;
01048 
01049   // Get the tag values at the vertices
01050   ErrorCode result = mbImpl->tag_get_data(tag, &elem, 1, &tempField);
01051   if (MB_SUCCESS != result)
01052     return result;
01053 
01054   field = tempField;
01055 
01056   return MB_SUCCESS;
01057 }
01058 
01059 // Normalize a field over the entire mesh represented by the root_set.
01060 ErrorCode Coupler::normalize_mesh(EntityHandle          root_set,
01061                                   const char            *norm_tag,
01062                                   Coupler::IntegType    integ_type,
01063                                   int                   num_integ_pts)
01064 {
01065   ErrorCode err;
01066 
01067   // SLAVE START ****************************************************************
01068   // Search for entities based on tag_handles and tag_values
01069   std::vector< std::vector<EntityHandle> > entity_sets;
01070   std::vector< std::vector<EntityHandle> > entity_groups;
01071 
01072   // put the root_set into entity_sets
01073   std::vector<EntityHandle> ent_set;
01074   ent_set.push_back(root_set);
01075   entity_sets.push_back(ent_set);
01076 
01077   // get all entities from root_set and put into entity_groups
01078   std::vector<EntityHandle> entities;
01079   err = mbImpl->get_entities_by_handle(root_set, entities, true);
01080   ERRORR("Failed to get entities in root_set.", err);
01081 
01082   entity_groups.push_back(entities);
01083 
01084   // Call do_normalization() to continue common normalization processing
01085   err = do_normalization(norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts);
01086   ERRORR("Failure in do_normalization().", err);
01087   // SLAVE END   ****************************************************************
01088 
01089   return err;
01090 }
01091 
01092 // Normalize a field over the subset of entities identified by the tags and values passed
01093 ErrorCode Coupler::normalize_subset(EntityHandle        root_set,
01094                                     const char          *norm_tag,
01095                                     const char          **tag_names,
01096                                     int                 num_tags,
01097                                     const char          **tag_values,
01098                                     Coupler::IntegType  integ_type,
01099                                     int                 num_integ_pts)
01100 {
01101   moab::ErrorCode err;
01102   std::vector<Tag> tag_handles;
01103 
01104   // Lookup tag handles from tag names
01105   for (int t = 0; t < num_tags; t++) {
01106     // get tag handle & size
01107     Tag th;
01108     err = mbImpl->tag_get_handle(tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY);
01109     ERRORR("Failed to get tag handle.", err);
01110     tag_handles.push_back(th);
01111   }
01112 
01113   return normalize_subset(root_set,
01114                           norm_tag,
01115                           &tag_handles[0],
01116                           num_tags,
01117                           tag_values,
01118                           integ_type,
01119                           num_integ_pts);
01120 }
01121 
01122 ErrorCode Coupler::normalize_subset(EntityHandle       root_set,
01123                                     const char         *norm_tag,
01124                                     Tag                *tag_handles,
01125                                     int                num_tags,
01126                                     const char         **tag_values,
01127                                     Coupler::IntegType integ_type,
01128                                     int                num_integ_pts)
01129 {
01130   ErrorCode err;
01131 
01132   // SLAVE START ****************************************************************
01133   // Search for entities based on tag_handles and tag_values
01134   std::vector< std::vector<EntityHandle> > entity_sets;
01135   std::vector< std::vector<EntityHandle> > entity_groups;
01136 
01137   err = get_matching_entities(root_set, tag_handles, tag_values, num_tags,
01138                               &entity_sets, &entity_groups);
01139   ERRORR("Failed to get matching entities.", err);
01140 
01141   // Call do_normalization() to continue common normalization processing
01142   err = do_normalization(norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts);
01143   ERRORR("Failure in do_normalization().", err);
01144   // SLAVE END   ****************************************************************
01145 
01146   return err;
01147 }
01148 
01149 ErrorCode Coupler::do_normalization(const char                               *norm_tag,
01150                                     std::vector<std::vector<EntityHandle> >  &entity_sets,
01151                                     std::vector<std::vector<EntityHandle> >  &entity_groups,
01152                                     Coupler::IntegType                       integ_type,
01153                                     int                                      num_integ_pts)
01154 {
01155   // SLAVE START ****************************************************************
01156   ErrorCode err;
01157   int ierr = 0;
01158 
01159   // Setup data for parallel computing
01160   int nprocs, rank;
01161   ierr = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
01162   ERRORMPI("Getting number of procs failed.", ierr);
01163   ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
01164   ERRORMPI("Getting rank failed.", ierr);
01165 
01166   // Get the integrated field value for each group(vector) of entities.
01167   // If no entities are in a group then a zero will be put in the list
01168   // of return values.
01169   unsigned int num_ent_grps = entity_groups.size();
01170   std::vector<double> integ_vals(num_ent_grps);
01171 
01172   err = get_group_integ_vals(entity_groups, integ_vals, norm_tag, num_integ_pts, integ_type);
01173   ERRORR("Failed to get integrated field values for groups in mesh.", err);
01174   // SLAVE END   ****************************************************************
01175 
01176   // SLAVE/MASTER START #########################################################
01177   // Send list of integrated values back to master proc. The ordering of the
01178   // values will match the ordering of the entity groups (i.e. vector of vectors)
01179   // sent from master to slaves earlier.  The values for each entity group will
01180   // be summed during the transfer.
01181   std::vector<double> sum_integ_vals(num_ent_grps);
01182 
01183   if (nprocs > 1) {
01184     // If parallel then send the values back to the master.
01185     ierr = MPI_Reduce(&integ_vals[0], &sum_integ_vals[0], num_ent_grps, MPI_DOUBLE,
01186                      MPI_SUM, MASTER_PROC, myPc->proc_config().proc_comm());
01187     ERRORMPI("Transfer and reduction of integrated values failed.", ierr);
01188   }
01189   else {
01190     // Otherwise just copy the vector
01191     sum_integ_vals = integ_vals;
01192   }
01193   // SLAVE/MASTER END   #########################################################
01194 
01195   // MASTER START ***************************************************************
01196   // Calculate the normalization factor for each group by taking the
01197   // inverse of each integrated field value. Put the normalization factor
01198   // for each group back into the list in the same order.
01199   for (unsigned int i = 0; i < num_ent_grps; i++) {
01200     double val = sum_integ_vals[i];
01201     if (fabs(val) > 1e-8) sum_integ_vals[i] = 1.0/val;
01202     else {
01203       sum_integ_vals[i] = 0.0; /* VSM: not sure what we should do here ? */
01204       /* commenting out error below since if integral(value)=0.0, then normalization
01205          is probably unnecessary to start with ? */
01206       /* ERRORR("Integrating an invalid field -- integral("<<norm_tag<<") = "<<val<<".", err); */
01207     }
01208   }
01209   // MASTER END   ***************************************************************
01210 
01211   // MASTER/SLAVE START #########################################################
01212   if (nprocs > 1) {
01213     // If parallel then broadcast the normalization factors to the procs.
01214     ierr = MPI_Bcast(&sum_integ_vals[0], num_ent_grps, MPI_DOUBLE, MASTER_PROC,
01215                     myPc->proc_config().proc_comm());
01216     ERRORMPI("Broadcast of normalization factors failed.", ierr);
01217   }
01218   // MASTER/SLAVE END   #########################################################
01219 
01220   // SLAVE START ****************************************************************
01221   // Save the normalization factors to a new tag with name of norm_tag's value
01222   // and the string "_normF" appended.  This new tag will be created on the entity
01223   // set that contains all of the entities from a group.
01224 
01225   err = apply_group_norm_factor(entity_sets, sum_integ_vals, norm_tag, integ_type);
01226   ERRORR("Failed to set the normalization factor for groups in mesh.", err);
01227   // SLAVE END   ****************************************************************
01228 
01229   return err;
01230 }
01231 
01232 // Functions supporting the subset normalization function
01233 
01234 // Retrieve groups of entities matching tags and values if present
01235 ErrorCode Coupler::get_matching_entities(EntityHandle                            root_set,
01236                                          const char                              **tag_names,
01237                                          const char                              **tag_values,
01238                                          int                                     num_tags,
01239                                          std::vector<std::vector<EntityHandle> > *entity_sets,
01240                                          std::vector<std::vector<EntityHandle> > *entity_groups)
01241 {
01242   ErrorCode err;
01243   std::vector<Tag> tag_handles;
01244 
01245   for (int t = 0; t < num_tags; t++) {
01246     // Get tag handle & size
01247     Tag th;
01248     err = mbImpl->tag_get_handle(tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY);
01249     ERRORR("Failed to get tag handle.", err);
01250     tag_handles.push_back(th);
01251   }
01252 
01253   return get_matching_entities(root_set, &tag_handles[0], tag_values, num_tags,
01254                                entity_sets, entity_groups);
01255 }
01256 
01257 // Retrieve groups of entities matching tags and values if present
01258 ErrorCode Coupler::get_matching_entities(EntityHandle                            root_set,
01259                                          Tag                                     *tag_handles,
01260                                          const char                              **tag_values,
01261                                          int                                     num_tags,
01262                                          std::vector<std::vector<EntityHandle> > *entity_sets,
01263                                          std::vector<std::vector<EntityHandle> > *entity_groups)
01264 {
01265   // SLAVE START ****************************************************************
01266 
01267   // Setup data for parallel computing
01268   ErrorCode err;
01269   int ierr = 0;
01270   int nprocs, rank;
01271   ierr = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
01272   ERRORMPI("Getting number of procs failed.", ierr);
01273   ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
01274   ERRORMPI("Getting rank failed.", ierr);
01275 
01276   Range ent_sets;
01277   err = mbImpl->get_entities_by_type_and_tag(root_set, moab::MBENTITYSET, tag_handles,
01278                                              (const void* const *)tag_values,
01279                                              num_tags, ent_sets, Interface::INTERSECT, 0);
01280   ERRORR("Core::get_entities_by_type_and_tag failed.", err);
01281 
01282   TupleList *tag_list = NULL;
01283   err = create_tuples(ent_sets, tag_handles, num_tags, &tag_list);
01284   ERRORR("Failed to create tuples from entity sets.", err);
01285 
01286   // Free up range
01287   ent_sets.clear();
01288   // SLAVE END   ****************************************************************
01289 
01290   // If we are running in a multi-proc session then send tuple list back to master
01291   // proc for consolidation. Otherwise just copy the pointer to the tuple_list.
01292   TupleList *cons_tuples;
01293   if (nprocs > 1) {
01294     // SLAVE/MASTER START #########################################################
01295 
01296     // Pack the tuple_list in a buffer.
01297     uint *tuple_buf;
01298     int tuple_buf_len;
01299     tuple_buf_len = pack_tuples(tag_list, (void**)&tuple_buf);
01300 
01301     // Free tag_list here as its not used again if nprocs > 1
01302     tag_list->reset();
01303 
01304     // Send back the buffer sizes to the master proc
01305     int *recv_cnts = (int*)malloc(nprocs * sizeof(int));
01306     int *offsets   = (int*)malloc(nprocs * sizeof(int));
01307     uint *all_tuples_buf = NULL;
01308 
01309     MPI_Gather(&tuple_buf_len, 1, MPI_INT, recv_cnts, 1, MPI_INT, MASTER_PROC,
01310                myPc->proc_config().proc_comm());
01311     ERRORMPI("Gathering buffer sizes failed.", err);
01312 
01313     // Allocate a buffer large enough for all the data
01314     if (rank == MASTER_PROC) {
01315       int all_tuples_len = recv_cnts[0];
01316       offsets[0] = 0;
01317       for (int i = 1; i < nprocs; i++) {
01318         offsets[i] = offsets[i - 1] + recv_cnts[i - 1];
01319         all_tuples_len += recv_cnts[i];
01320       }
01321 
01322       all_tuples_buf = (uint*)malloc(all_tuples_len * sizeof(uint));
01323     }
01324 
01325     // Send all buffers to the master proc for consolidation
01326     MPI_Gatherv((void*)tuple_buf, tuple_buf_len, MPI_INT,
01327                 (void*)all_tuples_buf, recv_cnts, offsets, MPI_INT, MASTER_PROC,
01328                 myPc->proc_config().proc_comm());
01329     ERRORMPI("Gathering tuple_lists failed.", err);
01330     free(tuple_buf); // malloc'd in pack_tuples
01331 
01332     if (rank == MASTER_PROC) {
01333       // Unpack the tuple_list from the buffer.
01334       TupleList **tl_array = (TupleList**)malloc(nprocs * sizeof(TupleList*));
01335       for (int i = 0; i < nprocs; i++)
01336         unpack_tuples((void*) &all_tuples_buf[offsets[i]], &tl_array[i]);
01337 
01338       // Free all_tuples_buf here as it is only allocated on the MASTER_PROC
01339       free(all_tuples_buf);
01340       // SLAVE/MASTER END   #########################################################
01341 
01342       // MASTER START ***************************************************************
01343       // Consolidate all tuple_lists into one tuple_list with no duplicates.
01344       err = consolidate_tuples(tl_array, nprocs, &cons_tuples);
01345       ERRORR("Failed to consolidate tuples.", err);
01346 
01347       for (int i = 0; i < nprocs; i++)
01348         tl_array[i]->reset();
01349       free(tl_array);
01350       // MASTER END   ***************************************************************
01351     }
01352 
01353     // Free offsets and recv_cnts as they are allocated on all procs
01354     free(offsets);
01355     free(recv_cnts);
01356 
01357     // MASTER/SLAVE START #########################################################
01358     // Broadcast condensed tuple list back to all procs.
01359     uint *ctl_buf;
01360     int ctl_buf_sz;
01361     if (rank == MASTER_PROC)
01362       ctl_buf_sz = pack_tuples(cons_tuples, (void**)&ctl_buf);
01363 
01364     // Send buffer size
01365     ierr = MPI_Bcast(&ctl_buf_sz, 1, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm());
01366     ERRORMPI("Broadcasting tuple_list size failed.", ierr);
01367 
01368     // Allocate a buffer in the other procs
01369     if (rank != MASTER_PROC)
01370       ctl_buf = (uint*)malloc(ctl_buf_sz * sizeof(uint));
01371 
01372     ierr = MPI_Bcast((void*)ctl_buf, ctl_buf_sz, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm());
01373     ERRORMPI("Broadcasting tuple_list failed.", ierr);
01374 
01375     if (rank != MASTER_PROC)
01376       unpack_tuples(ctl_buf, &cons_tuples);
01377     free(ctl_buf);
01378     // MASTER/SLAVE END   #########################################################
01379   }
01380   else
01381     cons_tuples = tag_list;
01382 
01383   // SLAVE START ****************************************************************
01384   // Loop over the tuple list getting the entities with the tags in the tuple_list entry
01385   uint mi, ml, mul, mr;
01386   cons_tuples->getTupleSize(mi, ml, mul, mr);
01387 
01388   for (unsigned int i = 0; i < cons_tuples->get_n(); i++) {
01389     // Get Entity Sets that match the tags and values.
01390 
01391     // Convert the data in the tuple_list to an array of pointers to the data
01392     // in the tuple_list as that is what the iMesh API call is expecting.
01393     int **vals = (int**)malloc(mi * sizeof(int*));
01394     for (unsigned int j = 0; j < mi; j++)
01395       vals[j] = (int *)&(cons_tuples->vi_rd[(i*mi) + j]);
01396 
01397     // Get entities recursively based on type and tag data
01398     err = mbImpl->get_entities_by_type_and_tag(root_set, moab::MBENTITYSET, tag_handles,
01399                                                (const void* const *)vals,
01400                                                mi, ent_sets, Interface::INTERSECT, 0);
01401     ERRORR("Core::get_entities_by_type_and_tag failed.", err);
01402     if (debug) std::cout << "ent_sets_size=" << ent_sets.size() << std::endl;
01403 
01404     // Free up the array of pointers
01405     free(vals);
01406 
01407     // Loop over the entity sets and then free the memory for ent_sets.
01408     std::vector<EntityHandle> ent_set_hdls;
01409     std::vector<EntityHandle> ent_hdls;
01410     for (unsigned int j = 0; j < ent_sets.size(); j++) {
01411       // Save the entity set
01412       ent_set_hdls.push_back(ent_sets[j]);
01413 
01414       // Get all entities for the entity set
01415       Range ents;
01416 
01417       /* VSM: do we need to filter out entity sets ? */
01418       err = mbImpl->get_entities_by_handle(ent_sets[j], ents, false);
01419       ERRORR("Core::get_entities_by_handle failed.", err);
01420       if (debug) std::cout << "ents_size=" << ents.size() << std::endl;
01421 
01422       // Save all of the entities from the entity set and free the memory for ents.
01423       for (unsigned int k = 0; k < ents.size(); k++) {
01424         ent_hdls.push_back(ents[k]);
01425       }
01426       ents.clear();
01427       if (debug) std::cout << "ent_hdls.size=" << ent_hdls.size() << std::endl;
01428     }
01429 
01430     // Free the entity set list for next tuple iteration.
01431     ent_sets.clear();
01432 
01433     // Push ent_set_hdls onto entity_sets, ent_hdls onto entity_groups
01434     // and clear both ent_set_hdls and ent_hdls.
01435     entity_sets->push_back(ent_set_hdls);
01436     ent_set_hdls.clear();
01437     entity_groups->push_back(ent_hdls);
01438     ent_hdls.clear();
01439     if (debug) std::cout << "entity_sets->size=" << entity_sets->size()
01440                          << ", entity_groups->size=" << entity_groups->size() << std::endl;
01441   }
01442 
01443   cons_tuples->reset();
01444   // SLAVE END   ****************************************************************
01445 
01446   return err;
01447 }
01448 
01449 // Return a tuple_list containing  tag values for each Entity Set
01450 // The tuple_list will have a column for each tag and a row for each
01451 // Entity Set. It is assumed all of the tags are integer tags.
01452 ErrorCode Coupler::create_tuples(Range        &ent_sets,
01453                                  const char   **tag_names,
01454                                  unsigned int num_tags,
01455                                  TupleList    **tuple_list)
01456 {
01457   ErrorCode err;
01458   std::vector<Tag> tag_handles;
01459 
01460   for (unsigned int t = 0; t < num_tags; t++) {
01461     // Get tag handle & size
01462     Tag th;
01463     err = mbImpl->tag_get_handle(tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY);
01464     ERRORR("Failed to get tag handle.", err);
01465     tag_handles.push_back(th);
01466   }
01467 
01468   return create_tuples(ent_sets, &tag_handles[0], num_tags, tuple_list);
01469 }
01470 
01471 // Return a tuple_list containing  tag values for each Entity Set
01472 // The tuple_list will have a column for each tag and a row for each
01473 // Entity Set.  It is assumed all of the tags are integer tags.
01474 ErrorCode Coupler::create_tuples(Range        &ent_sets,
01475                                  Tag          *tag_handles,
01476                                  unsigned int num_tags,
01477                                  TupleList    **tuples)
01478 {
01479   // ASSUMPTION: All tags are of type integer.  This may need to be expanded in future.
01480   ErrorCode err;
01481 
01482   // Allocate a tuple_list for the number of entity sets passed in
01483   TupleList *tag_tuples = new TupleList(num_tags, 0, 0, 0, (int)ent_sets.size());
01484   //tag_tuples->initialize(num_tags, 0, 0, 0, num_sets);
01485   uint mi, ml, mul, mr;
01486   tag_tuples->getTupleSize(mi, ml, mul, mr);
01487   tag_tuples->enableWriteAccess();
01488 
01489   if (mi == 0)
01490     ERRORR("Failed to initialize tuple_list.", MB_FAILURE);
01491 
01492   // Loop over the filtered entity sets retrieving each matching tag value one by one.
01493   int val;
01494   for (unsigned int i = 0; i < ent_sets.size(); i++) {
01495     for (unsigned int j = 0; j < num_tags; j++) {
01496       EntityHandle set_handle = ent_sets[i];
01497       err = mbImpl->tag_get_data(tag_handles[j], &set_handle, 1, &val);
01498       ERRORR("Failed to get integer tag data.", err);
01499       tag_tuples->vi_wr[i*mi + j] = val;
01500     }
01501 
01502     // If we get here there was no error so increment n in the tuple_list
01503     tag_tuples->inc_n();
01504   }
01505   tag_tuples->disableWriteAccess();
01506   *tuples = tag_tuples;
01507 
01508   return MB_SUCCESS;
01509 }
01510 
01511 // Consolidate tuple_lists into one list with no duplicates
01512 ErrorCode Coupler::consolidate_tuples(TupleList     **all_tuples,
01513                                       unsigned int  num_tuples,
01514                                       TupleList     **unique_tuples)
01515 {
01516   int total_rcv_tuples = 0;
01517   int offset = 0, copysz = 0;
01518   unsigned num_tags = 0;
01519 
01520   uint ml, mul, mr;
01521   uint *mi = (uint *)malloc(sizeof(uint) * num_tuples);
01522 
01523   for(unsigned int i = 0; i < num_tuples; i++) {
01524     all_tuples[i]->getTupleSize(mi[i], ml, mul, mr);
01525   }
01526 
01527   for (unsigned int i = 0; i < num_tuples; i++) {
01528     if (all_tuples[i] != NULL) {
01529       total_rcv_tuples += all_tuples[i]->get_n();
01530       num_tags = mi[i];
01531     }
01532   }
01533   const unsigned int_size = sizeof(sint);
01534   const unsigned int_width = num_tags * int_size;
01535 
01536   // Get the total size of all of the tuple_lists in all_tuples.
01537   for (unsigned int i = 0; i < num_tuples; i++) {
01538     if (all_tuples[i] != NULL)
01539       total_rcv_tuples += all_tuples[i]->get_n();
01540   }
01541 
01542   // Copy the tuple_lists into a single tuple_list.
01543   TupleList *all_tuples_list = new TupleList(num_tags, 0, 0, 0, total_rcv_tuples);
01544   all_tuples_list->enableWriteAccess();
01545   //all_tuples_list->initialize(num_tags, 0, 0, 0, total_rcv_tuples);
01546   for (unsigned int i = 0; i < num_tuples; i++) {
01547     if (all_tuples[i] != NULL) {
01548       copysz = all_tuples[i]->get_n() * int_width;
01549       memcpy(all_tuples_list->vi_wr+offset, all_tuples[i]->vi_rd, copysz);
01550       offset = offset + (all_tuples[i]->get_n() * mi[i]);
01551       all_tuples_list->set_n(all_tuples_list->get_n() + all_tuples[i]->get_n());
01552     }
01553   }
01554 
01555   // Sort the new tuple_list. Use a radix type sort, starting with the last (or least significant)
01556   // tag column in the vi array and working towards the first (or most significant) tag column.
01557   TupleList::buffer sort_buffer;
01558   sort_buffer.buffer_init(2 * total_rcv_tuples * int_width);
01559   for (int i = num_tags - 1; i >= 0; i--) {
01560     all_tuples_list->sort(i, &sort_buffer);
01561   }
01562 
01563   // Cycle through the sorted list eliminating duplicates.
01564   // Keep counters to the current end of the tuple_list (w/out dups) and the last tuple examined.
01565   unsigned int end_idx = 0, last_idx = 1;
01566   while (last_idx < all_tuples_list->get_n()) {
01567     if (memcmp(all_tuples_list->vi_rd + (end_idx*num_tags), all_tuples_list->vi_rd + (last_idx*num_tags), int_width) == 0) {
01568       // Values equal - skip
01569       last_idx += 1;
01570     }
01571     else {
01572       // Values different - copy
01573       // Move up the end index
01574       end_idx += 1;
01575       memcpy(all_tuples_list->vi_wr + (end_idx*num_tags), all_tuples_list->vi_rd + (last_idx*num_tags), int_width);
01576       last_idx += 1;
01577     }
01578   }
01579   // Update the count in all_tuples_list
01580   all_tuples_list->set_n(end_idx + 1);
01581 
01582   // Resize the tuple_list
01583   all_tuples_list->resize(all_tuples_list->get_n());
01584 
01585   // Set the output parameter
01586   *unique_tuples = all_tuples_list;
01587 
01588   return MB_SUCCESS;
01589 }
01590 
01591 // Calculate integrated field values for groups of entities
01592 ErrorCode Coupler::get_group_integ_vals(std::vector<std::vector<EntityHandle> > &groups,
01593                                         std::vector<double> &integ_vals,
01594                                         const char *norm_tag,
01595                                         int /*num_integ_vals*/,
01596                                         Coupler::IntegType integ_type)
01597 {
01598   ErrorCode err;
01599 
01600   std::vector<std::vector<EntityHandle> >::iterator iter_i;
01601   std::vector<EntityHandle>::iterator iter_j;
01602   double grp_intrgr_val, intgr_val;
01603 
01604   // Get the tag handle for norm_tag
01605   Tag norm_hdl;
01606   err = mbImpl->tag_get_handle(norm_tag, 1, moab::MB_TYPE_DOUBLE, norm_hdl, moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT);
01607   ERRORR("Failed to get norm_tag handle.", err);
01608 
01609   // Check size of integ_vals vector
01610   if (integ_vals.size() != groups.size())
01611     integ_vals.resize(groups.size());
01612 
01613   // Loop over the groups(vectors) of entities
01614   unsigned int i;
01615   for (i = 0, iter_i = groups.begin(); iter_i != groups.end(); i++, ++iter_i) {
01616     grp_intrgr_val = 0;
01617 
01618     // Loop over the all the entities in the group, integrating
01619     // the field_fn over the entity in iter_j
01620     for (iter_j = (*iter_i).begin(); iter_j != (*iter_i).end(); ++iter_j) {
01621       EntityHandle ehandle = (*iter_j);
01622 
01623       // Check that the entity in iter_j is of the same dimension as the
01624       // integ_type we are performing
01625       EntityType j_type;
01626       j_type = mbImpl->type_from_handle(ehandle);
01627       ERRORR("Failed to get entity type.", err);
01628       // Skip any entities in the group that are not of the type being considered
01629       if ((integ_type == VOLUME) && (j_type < MBTET || j_type >= MBENTITYSET))
01630         continue;
01631 
01632       intgr_val = 0;
01633 
01634       // Retrieve the vertices from the element
01635       const EntityHandle* verts = NULL;
01636       int connectivity_size = 0;
01637 
01638       err = mbImpl->get_connectivity(ehandle, verts, connectivity_size, false);
01639       ERRORR("Failed to get vertices from entity.", err);
01640 
01641       // Get the vertex coordinates and the field values at the vertices.
01642       double *coords = (double*) malloc(sizeof(double) * (3*connectivity_size));
01643       /* TODO: VSM: check if this works for lower dimensions also without problems */
01644       /* if (3 == geom_dim) */
01645       err = mbImpl->get_coords(verts, connectivity_size, coords);
01646       ERRORR("Failed to get vertex coordinates.", err);
01647 
01648       /* allocate the field data array */
01649       double *vfield = (double*) malloc(sizeof(double) * (connectivity_size));
01650       err = mbImpl->tag_get_data(norm_hdl, verts, connectivity_size, vfield);
01651       if (MB_SUCCESS != err) {
01652         free(coords);
01653       }
01654       ERRORR("Failed to get vertex coordinates.", err);
01655 
01656       // Get coordinates of all corner vertices (in normal order) and
01657       // put in array of CartVec.
01658       std::vector<CartVect> vertices(connectivity_size);
01659 
01660       // Put the vertices into a CartVect vector
01661       double *x = coords;
01662       for (int j = 0; j < connectivity_size; j++, x += 3) {
01663         vertices[j] = CartVect(x);
01664       }
01665       free(coords);
01666 
01667       moab::Element::Map *elemMap;
01668       if (j_type == MBHEX) {
01669         if (connectivity_size == 8)
01670           elemMap = new moab::Element::LinearHex(vertices);
01671         else
01672           elemMap = new moab::Element::QuadraticHex(vertices);
01673       }
01674       else if (j_type == MBTET) {
01675         elemMap = new moab::Element::LinearTet(vertices);
01676       }
01677       else if (j_type == MBQUAD) {
01678         elemMap = new moab::Element::LinearQuad(vertices);
01679       }
01680       /*
01681       else if (j_type == MBTRI) {
01682         elemMap = new moab::Element::LinearTri(vertices);
01683       }
01684       */
01685       else if (j_type == MBEDGE) {
01686         elemMap = new moab::Element::LinearEdge(vertices);
01687       }
01688       else
01689         ERRORR("Unknown topology type.", MB_UNSUPPORTED_OPERATION);
01690 
01691       // Set the vertices in the Map and perform the integration
01692       try {
01693         /* VSM: Do we need this call ?? */
01694         // elemMap->set_vertices(vertices);
01695 
01696         // Perform the actual integration over the element
01697         intgr_val = elemMap->integrate_scalar_field(vfield);
01698 
01699         // Combine the result with those of the group
01700         grp_intrgr_val += intgr_val;
01701       }
01702       catch (moab::Element::Map::ArgError&) {
01703         std::cerr << "Failed to set vertices on Element::Map." << std::endl;
01704       }
01705       catch (moab::Element::Map::EvaluationError&) {
01706         std::cerr << "Failed to get inverse evaluation of coordinate on Element::Map." << std::endl;
01707       }
01708 
01709       delete(elemMap);
01710       free(vfield);
01711     }
01712 
01713     // Set the group integrated value in the vector
01714     integ_vals[i] = grp_intrgr_val;
01715   }
01716 
01717   return err;
01718 }
01719 
01720 // Apply a normalization factor to group of entities
01721 ErrorCode Coupler::apply_group_norm_factor(std::vector<std::vector<EntityHandle> > &entity_sets,
01722                                            std::vector<double>                     &norm_factors,
01723                                            const char                              *norm_tag,
01724                                            Coupler::IntegType                      /*integ_type*/)
01725 {
01726   ErrorCode err;
01727 
01728   // Construct the new tag for the normalization factor from the norm_tag name
01729   // and "_normf".
01730   int norm_tag_len = strlen(norm_tag);
01731   const char* normf_appd = "_normf";
01732   int normf_appd_len = strlen(normf_appd);
01733 
01734   char* normf_tag = (char*)malloc(norm_tag_len + normf_appd_len + 1);
01735   char* tmp_ptr = normf_tag;
01736 
01737   memcpy(tmp_ptr, norm_tag, norm_tag_len);
01738   tmp_ptr += norm_tag_len;
01739   memcpy(tmp_ptr, normf_appd, normf_appd_len);
01740   tmp_ptr += normf_appd_len;
01741   *tmp_ptr = '\0';
01742 
01743   Tag normf_hdl;
01744   // Check to see if the tag exists.  If not then create it and get the handle.
01745   err = mbImpl->tag_get_handle(normf_tag, 1, moab::MB_TYPE_DOUBLE, normf_hdl, moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT);
01746   ERRORR("Failed to create normalization factor tag.", err);
01747   if (normf_hdl == NULL) {
01748     std::string msg("Failed to create normalization factor tag named '");
01749     msg += std::string(normf_tag) + std::string("'");
01750     ERRORR(msg.c_str(), MB_FAILURE);
01751   }
01752   free(normf_tag);
01753 
01754   std::vector< std::vector<EntityHandle> >::iterator iter_i;
01755   std::vector<EntityHandle>::iterator iter_j;
01756   std::vector<double>::iterator iter_f;
01757   double grp_norm_factor = 0.0;
01758 
01759   // Loop over the entity sets
01760   for (iter_i = entity_sets.begin(), iter_f = norm_factors.begin();
01761        (iter_i != entity_sets.end()) && (iter_f != norm_factors.end());
01762        ++iter_i, ++iter_f) {
01763     grp_norm_factor = *iter_f;
01764 
01765     // Loop over the all the entity sets in iter_i and set the
01766     // new normf_tag with the norm factor value on each
01767     for (iter_j = (*iter_i).begin(); iter_j != (*iter_i).end(); ++iter_j) {
01768       EntityHandle entset = *iter_j;
01769 
01770       std::cout << "Coupler: applying normalization for entity set="
01771                 << entset << ",  normalization_factor=" << grp_norm_factor << std::endl;
01772 
01773       err = mbImpl->tag_set_data(normf_hdl, &entset, 1, &grp_norm_factor);
01774       ERRORR("Failed to set normalization factor on entity set.", err);
01775     }
01776   }
01777 
01778   return MB_SUCCESS;
01779 }
01780 
01781 #define UINT_PER_X(X) ((sizeof(X) + sizeof(uint) - 1) / sizeof(uint))
01782 #define UINT_PER_REAL UINT_PER_X(realType)
01783 #define UINT_PER_LONG UINT_PER_X(slong)
01784 #define UINT_PER_UNSIGNED UINT_PER_X(unsigned)
01785 
01786 // Function for packing tuple_list.  Returns number of uints copied into buffer.
01787 int pack_tuples(TupleList* tl, void** ptr)
01788 {
01789   uint mi, ml, mul, mr;
01790   tl->getTupleSize(mi, ml, mul, mr);
01791 
01792   uint n = tl->get_n();
01793 
01794   int sz_buf = 1 + 4*UINT_PER_UNSIGNED +
01795                tl->get_n() * (mi +
01796                               ml*UINT_PER_LONG +
01797                               mul*UINT_PER_LONG +
01798                               mr*UINT_PER_REAL);
01799 
01800   uint *buf = (uint*) malloc(sz_buf*sizeof(uint));
01801   *ptr = (void*) buf;
01802 
01803   // Copy n
01804   memcpy(buf, &n,   sizeof(uint)),     buf += 1;
01805   // Copy mi
01806   memcpy(buf, &mi,  sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
01807   // Copy ml
01808   memcpy(buf, &ml,  sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
01809   // Copy mul
01810   memcpy(buf, &mul, sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
01811   // Copy mr
01812   memcpy(buf, &mr,  sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
01813   // Copy vi_wr
01814   memcpy(buf, tl->vi_rd,  tl->get_n()*mi*sizeof(sint)),   buf += tl->get_n()*mi;
01815   // Copy vl_wr
01816   memcpy(buf, tl->vl_rd,  tl->get_n()*ml*sizeof(slong)),  buf += tl->get_n()*ml*UINT_PER_LONG;
01817   // Copy vul_wr
01818   memcpy(buf, tl->vul_rd, tl->get_n()*mul*sizeof(ulong)), buf += tl->get_n()*mul*UINT_PER_LONG;
01819   // Copy vr_wr
01820   memcpy(buf, tl->vr_rd,  tl->get_n()*mr*sizeof(realType)),   buf += tl->get_n()*mr*UINT_PER_REAL;
01821 
01822   return sz_buf;
01823 }
01824 
01825 // Function for packing tuple_list
01826 void unpack_tuples(void* ptr, TupleList** tlp)
01827 {
01828   TupleList *tl = new TupleList();
01829   *tlp = tl;
01830 
01831   uint nt;
01832   unsigned mit, mlt, mult, mrt;
01833   uint *buf = (uint*)ptr;
01834 
01835   // Get n
01836   memcpy(&nt,   buf, sizeof(uint)),     buf += 1;
01837   // Get mi
01838   memcpy(&mit,  buf, sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
01839   // Get ml
01840   memcpy(&mlt,  buf, sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
01841   // Get mul
01842   memcpy(&mult, buf, sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
01843   // Get mr
01844   memcpy(&mrt,  buf, sizeof(unsigned)), buf += UINT_PER_UNSIGNED;
01845 
01846   // Initialize tl
01847   tl->initialize(mit, mlt, mult, mrt, nt);
01848   tl->enableWriteAccess();
01849   tl->set_n(nt);
01850 
01851   uint mi, ml, mul, mr;
01852   tl->getTupleSize(mi, ml, mul, mr);
01853 
01854   // Get vi_rd
01855   memcpy(tl->vi_wr,  buf, tl->get_n()*mi*sizeof(sint)),   buf += tl->get_n()*mi;
01856   // Get vl_rd
01857   memcpy(tl->vl_wr,  buf, tl->get_n()*ml*sizeof(slong)),  buf += tl->get_n()*ml*UINT_PER_LONG;
01858   // Get vul_rd
01859   memcpy(tl->vul_wr, buf, tl->get_n()*mul*sizeof(ulong)), buf += tl->get_n()*mul*UINT_PER_LONG;
01860   // Get vr_rd
01861   memcpy(tl->vr_wr,  buf, tl->get_n()*mr*sizeof(realType)),   buf += tl->get_n()*mr*UINT_PER_REAL;
01862 
01863   tl->disableWriteAccess();
01864   return;
01865 }
01866 
01867 ErrorCode Coupler::get_gl_points_on_elements(Range &targ_elems, std::vector<double> &vpos, int &numPointsOfInterest)
01868 {
01869   numPointsOfInterest = targ_elems.size() * _ntot;
01870   vpos.resize(3 * numPointsOfInterest);
01871   int ielem = 0;
01872   for (Range::iterator eit = targ_elems.begin(); eit != targ_elems.end(); ++eit, ielem += _ntot * 3) {
01873     EntityHandle eh = *eit;
01874     const double * xval;
01875     const double * yval;
01876     const double * zval;
01877     ErrorCode rval = mbImpl-> tag_get_by_ptr(_xm1Tag, &eh, 1, (const void**)&xval);
01878     if (moab::MB_SUCCESS != rval) {
01879       std::cout << "Can't get xm1 values \n";
01880       return MB_FAILURE;
01881     }
01882     rval = mbImpl-> tag_get_by_ptr(_ym1Tag, &eh, 1, (const void**)&yval);
01883     if (moab::MB_SUCCESS != rval) {
01884       std::cout << "Can't get ym1 values \n";
01885       return MB_FAILURE;
01886     }
01887     rval = mbImpl-> tag_get_by_ptr(_zm1Tag, &eh, 1, (const void**)&zval);
01888     if (moab::MB_SUCCESS != rval) {
01889       std::cout << "Can't get zm1 values \n";
01890       return MB_FAILURE;
01891     }
01892     // Now, in a stride, populate vpos
01893     for (int i = 0; i < _ntot; i++) {
01894       vpos[ielem + 3*i    ] = xval[i];
01895       vpos[ielem + 3*i + 1] = yval[i];
01896       vpos[ielem + 3*i + 2] = zval[i];
01897     }
01898   }
01899 
01900   return MB_SUCCESS;
01901 }
01902 
01903 } // namespace_moab
01904 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines