Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
Coupler.cpp
Go to the documentation of this file.
1 #include "Coupler.hpp"
2 #include "moab/ParallelComm.hpp"
4 #include "ElemUtil.hpp"
5 #include "moab/CN.hpp"
6 #include "moab/gs.hpp"
7 #include "moab/TupleList.hpp"
8 #include "moab/Error.hpp"
9 
10 /* C++ includes */
11 #include <cstdio>
12 #include <cassert>
13 #include <iostream>
14 #include <algorithm>
15 #include <sstream>
16 
17 #define ERROR( a ) \
18  { \
19  if( MB_SUCCESS != err ) std::cerr << ( a ) << std::endl; \
20  }
21 #define ERRORR( a, b ) \
22  { \
23  if( MB_SUCCESS != ( b ) ) \
24  { \
25  std::cerr << ( a ) << std::endl; \
26  return b; \
27  } \
28  }
29 #define ERRORMPI( a, b ) \
30  { \
31  if( MPI_SUCCESS != ( b ) ) \
32  { \
33  std::cerr << ( a ) << std::endl; \
34  return MB_FAILURE; \
35  } \
36  }
37 
38 #define MASTER_PROC 0
39 
40 namespace moab
41 {
42 
43 bool debug = false;
44 int pack_tuples( TupleList* tl, void** ptr );
45 void unpack_tuples( void* ptr, TupleList** tlp );
46 
48  ParallelComm* pc,
49  Range& local_elems,
50  int coupler_id,
51  bool init_tree,
52  int max_ent_dim )
53  : mbImpl( impl ), myPc( pc ), myId( coupler_id ), numIts( 3 ), max_dim( max_ent_dim ), _ntot( 0 ),
54  spherical( false )
55 {
56  assert( NULL != impl && ( pc || !local_elems.empty() ) );
57 
58  // Keep track of the local points, at least for now
59  myRange = local_elems;
60  myTree = NULL;
61 
62  // Now initialize the tree
63  if( init_tree ) initialize_tree();
64 
65  // Initialize tuple lists to indicate not initialized
66  mappedPts = NULL;
67  targetPts = NULL;
69 }
70 
71 /* Destructor
72  */
74 {
75  // This will clear the cache
78  delete myTree;
79  delete targetPts;
80  delete mappedPts;
81 }
82 
84 {
85  Range local_ents;
86 
87  // Get entities on the local part
88  ErrorCode result = MB_SUCCESS;
89  if( myPc )
90  {
91  result = myPc->get_part_entities( local_ents, max_dim );
92  if( local_ents.empty() )
93  {
94  max_dim--;
95  result = myPc->get_part_entities( local_ents, max_dim ); // go one dimension lower
96  // right now, this is used for spherical meshes only
97  }
98  }
99  else
100  local_ents = myRange;
101  if( MB_SUCCESS != result || local_ents.empty() )
102  {
103  std::cout << "Problems getting source entities" << std::endl;
104  return result;
105  }
106 
107  // Build the tree for local processor
108  int max_per_leaf = 6;
109  for( int i = 0; i < numIts; i++ )
110  {
111  std::ostringstream str;
112  str << "PLANE_SET=0;"
113  << "MAX_PER_LEAF=" << max_per_leaf << ";";
114  if( spherical && !local_ents.empty() )
115  {
116  // get coordinates of one vertex, and use for radius computation
117  EntityHandle elem = local_ents[0];
118  const EntityHandle* conn;
119  int numn = 0;
120  mbImpl->get_connectivity( elem, conn, numn );
121  CartVect pos0;
122  mbImpl->get_coords( conn, 1, &( pos0[0] ) );
123  double radius = pos0.length();
124  str << "SPHERICAL=true;RADIUS=" << radius << ";";
125  }
126  FileOptions opts( str.str().c_str() );
127  myTree = new AdaptiveKDTree( mbImpl );
128  result = myTree->build_tree( local_ents, &localRoot, &opts );
129  if( MB_SUCCESS != result )
130  {
131  std::cout << "Problems building tree";
132  if( numIts != i )
133  {
134  delete myTree;
135  max_per_leaf *= 2;
136  std::cout << "; increasing elements/leaf to " << max_per_leaf << std::endl;
137  }
138  else
139  {
140  std::cout << "; exiting" << std::endl;
141  return result;
142  }
143  }
144  else
145  break; // Get out of tree building
146  }
147 
148  // Get the bounding box for local tree
149  if( myPc )
150  allBoxes.resize( 6 * myPc->proc_config().proc_size() );
151  else
152  allBoxes.resize( 6 );
153 
154  unsigned int my_rank = ( myPc ? myPc->proc_config().proc_rank() : 0 );
155  BoundBox box;
156  result = myTree->get_bounding_box( box, &localRoot );
157  if( MB_SUCCESS != result ) return result;
158  box.bMin.get( &allBoxes[6 * my_rank] );
159  box.bMax.get( &allBoxes[6 * my_rank + 3] );
160 
161  // Now communicate to get all boxes
162  if( myPc )
163  {
164  int mpi_err;
165 #if( MPI_VERSION >= 2 )
166  // Use "in place" option
167  mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 6, MPI_DOUBLE,
168  myPc->proc_config().proc_comm() );
169 #else
170  {
171  std::vector< double > allBoxes_tmp( 6 * myPc->proc_config().proc_size() );
172  mpi_err = MPI_Allgather( &allBoxes[6 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE,
173  myPc->proc_config().proc_comm() );
174  allBoxes = allBoxes_tmp;
175  }
176 #endif
177  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
178  }
179 
180  /* std::ostringstream blah;
181  for(int i = 0; i < allBoxes.size(); i++)
182  blah << allBoxes[i] << " ";
183  std::cout << blah.str() << "\n";*/
184 
185 #ifdef VERBOSE
186  double min[3] = { 0, 0, 0 }, max[3] = { 0, 0, 0 };
187  unsigned int dep;
188  myTree->get_info( localRoot, min, max, dep );
189  std::cout << "Proc " << my_rank << ": box min/max, tree depth = (" << min[0] << "," << min[1] << "," << min[2]
190  << "), (" << max[0] << "," << max[1] << "," << max[2] << "), " << dep << std::endl;
191 #endif
192 
193  return result;
194 }
195 
197  EntityHandle rootTarget,
198  bool& specSou,
199  bool& specTar )
200 {
201  /*void * _spectralSource;
202  void * _spectralTarget;*/
203 
204  moab::Range spectral_sets;
205  moab::Tag sem_tag;
206  int sem_dims[3];
207  ErrorCode rval = mbImpl->tag_get_handle( "SEM_DIMS", 3, moab::MB_TYPE_INTEGER, sem_tag );
208  if( moab::MB_SUCCESS != rval )
209  {
210  std::cout << "Can't find tag, no spectral set\n";
211  return MB_SUCCESS; // Nothing to do, no spectral elements
212  }
213  rval = mbImpl->get_entities_by_type_and_tag( rootSource, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets );
214  if( moab::MB_SUCCESS != rval || spectral_sets.empty() )
215  std::cout << "Can't get sem set on source\n";
216  else
217  {
218  moab::EntityHandle firstSemSet = spectral_sets[0];
219  rval = mbImpl->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims );
220  if( moab::MB_SUCCESS != rval ) return MB_FAILURE;
221 
222  if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] )
223  {
224  std::cout << " dimensions are different. bail out\n";
225  return MB_FAILURE;
226  }
227  // Repeat for target sets
228  spectral_sets.empty();
229  // Now initialize a source spectral element !
230  _spectralSource = new moab::Element::SpectralHex( sem_dims[0] );
231  specSou = true;
232  }
233  // Repeat for target source
234  rval = mbImpl->get_entities_by_type_and_tag( rootTarget, moab::MBENTITYSET, &sem_tag, NULL, 1, spectral_sets );
235  if( moab::MB_SUCCESS != rval || spectral_sets.empty() )
236  std::cout << "Can't get sem set on target\n";
237  else
238  {
239  moab::EntityHandle firstSemSet = spectral_sets[0];
240  rval = mbImpl->tag_get_data( sem_tag, &firstSemSet, 1, (void*)sem_dims );
241  if( moab::MB_SUCCESS != rval ) return MB_FAILURE;
242 
243  if( sem_dims[0] != sem_dims[1] || sem_dims[0] != sem_dims[2] )
244  {
245  std::cout << " dimensions are different. bail out\n";
246  return MB_FAILURE;
247  }
248  // Repeat for target sets
249  spectral_sets.empty();
250  // Now initialize a target spectral element !
251  _spectralTarget = new moab::Element::SpectralHex( sem_dims[0] );
252  specTar = true;
253  }
254 
255  _ntot = sem_dims[0] * sem_dims[1] * sem_dims[2];
257  if( moab::MB_SUCCESS != rval )
258  {
259  std::cout << "Can't get xm1tag \n";
260  return MB_FAILURE;
261  }
263  if( moab::MB_SUCCESS != rval )
264  {
265  std::cout << "Can't get ym1tag \n";
266  return MB_FAILURE;
267  }
269  if( moab::MB_SUCCESS != rval )
270  {
271  std::cout << "Can't get zm1tag \n";
272  return MB_FAILURE;
273  }
274 
275  return MB_SUCCESS;
276 }
277 
278 ErrorCode Coupler::locate_points( Range& targ_ents, double rel_eps, double abs_eps, TupleList* tl, bool store_local )
279 {
280  // Get locations
281  std::vector< double > locs( 3 * targ_ents.size() );
282  Range verts = targ_ents.subset_by_type( MBVERTEX );
283  ErrorCode rval = mbImpl->get_coords( verts, &locs[0] );
284  if( MB_SUCCESS != rval ) return rval;
285  // Now get other ents; reuse verts
286  unsigned int num_verts = verts.size();
287  verts = subtract( targ_ents, verts );
288  // Compute centroids
289  std::vector< EntityHandle > dum_conn( CN::MAX_NODES_PER_ELEMENT );
290  std::vector< double > dum_pos( CN::MAX_NODES_PER_ELEMENT );
291  const EntityHandle* conn;
292  int num_conn;
293  double* coords = &locs[num_verts];
294  // Do this here instead of a function to allow reuse of dum_pos and dum_conn
295  for( Range::const_iterator rit = verts.begin(); rit != verts.end(); ++rit )
296  {
297  rval = mbImpl->get_connectivity( *rit, conn, num_conn, false, &dum_conn );
298  if( MB_SUCCESS != rval ) return rval;
299  rval = mbImpl->get_coords( conn, num_conn, &dum_pos[0] );
300  if( MB_SUCCESS != rval ) return rval;
301  coords[0] = coords[1] = coords[2] = 0.0;
302  for( int i = 0; i < num_conn; i++ )
303  {
304  coords[0] += dum_pos[3 * i];
305  coords[1] += dum_pos[3 * i + 1];
306  coords[2] += dum_pos[3 * i + 2];
307  }
308  coords[0] /= num_conn;
309  coords[1] /= num_conn;
310  coords[2] /= num_conn;
311  coords += 3;
312  }
313 
314  if( store_local ) targetEnts = targ_ents;
315 
316  return locate_points( &locs[0], targ_ents.size(), rel_eps, abs_eps, tl, store_local );
317 }
318 
320  unsigned int num_points,
321  double rel_eps,
322  double abs_eps,
323  TupleList* tl,
324  bool store_local )
325 {
326  assert( tl || store_local );
327 
328  // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing
329  // pts to be located source_pts: TL(from_proc, tgt_index, src_index): results of source mesh
330  // proc point location, ready to send
331  // back to tgt procs; src_index of -1 indicates point not located (arguably not
332  // useful...)
333  TupleList target_pts;
334  target_pts.initialize( 2, 0, 0, 3, num_points );
335  target_pts.enableWriteAccess();
336 
337  TupleList source_pts;
338  mappedPts = new TupleList( 0, 0, 1, 3, target_pts.get_max() );
340 
341  source_pts.initialize( 3, 0, 0, 0, target_pts.get_max() );
342  source_pts.enableWriteAccess();
343 
344  mappedPts->set_n( 0 );
345  source_pts.set_n( 0 );
346  ErrorCode result;
347 
348  unsigned int my_rank = ( myPc ? myPc->proc_config().proc_rank() : 0 );
349  bool point_located;
350 
351  // For each point, find box(es) containing the point,
352  // appending results to tuple_list;
353  // keep local points separately, in local_pts, which has pairs
354  // of <local_index, mapped_index>, where mapped_index is the index
355  // into the mappedPts tuple list
356  for( unsigned int i = 0; i < 3 * num_points; i += 3 )
357  {
358 
359  std::vector< int > procs_to_send_to;
360  for( unsigned int j = 0; j < ( myPc ? myPc->proc_config().proc_size() : 0 ); j++ )
361  {
362  // Test if point is in proc's box
363  if( ( allBoxes[6 * j] <= xyz[i] + abs_eps ) && ( xyz[i] <= allBoxes[6 * j + 3] + abs_eps ) &&
364  ( allBoxes[6 * j + 1] <= xyz[i + 1] + abs_eps ) && ( xyz[i + 1] <= allBoxes[6 * j + 4] + abs_eps ) &&
365  ( allBoxes[6 * j + 2] <= xyz[i + 2] + abs_eps ) && ( xyz[i + 2] <= allBoxes[6 * j + 5] + abs_eps ) )
366  {
367  // If in this proc's box, will send to proc to test further
368  procs_to_send_to.push_back( j );
369  }
370  }
371  if( procs_to_send_to.empty() )
372  {
373 #ifdef VERBOSE
374  std::cout << " point index " << i / 3 << ": " << xyz[i] << " " << xyz[i + 1] << " " << xyz[i + 2]
375  << " not found in any box\n";
376 #endif
377  // try to find the closest box, and put it in that box, anyway
378  double min_dist = 1.e+20;
379  int index = -1;
380  for( unsigned int j = 0; j < ( myPc ? myPc->proc_config().proc_size() : 0 ); j++ )
381  {
382  BoundBox box( &allBoxes[6 * j] ); // form back the box
383  double distance = box.distance( &xyz[i] ); // will compute the distance in 3d, from the box
384  if( distance < min_dist )
385  {
386  index = j;
387  min_dist = distance;
388  }
389  }
390  if( index == -1 )
391  {
392  // need to abort early, nothing we can do
393  assert( "cannot locate any box for some points" );
394  // need a better exit strategy
395  }
396 #ifdef VERBOSE
397  std::cout << " point index " << i / 3 << " added to box for proc j:" << index << "\n";
398 #endif
399  procs_to_send_to.push_back( index ); // will send to just one proc, that has the closest box
400  }
401  // we finally decided to populate the tuple list for a list of processors
402  for( size_t k = 0; k < procs_to_send_to.size(); k++ )
403  {
404  unsigned int j = procs_to_send_to[k];
405  // Check size of tuple list, grow if we're at max
406  if( target_pts.get_n() == target_pts.get_max() )
407  target_pts.resize( std::max( 10.0, 1.5 * target_pts.get_max() ) );
408 
409  target_pts.vi_wr[2 * target_pts.get_n()] = j;
410  target_pts.vi_wr[2 * target_pts.get_n() + 1] = i / 3;
411 
412  target_pts.vr_wr[3 * target_pts.get_n()] = xyz[i];
413  target_pts.vr_wr[3 * target_pts.get_n() + 1] = xyz[i + 1];
414  target_pts.vr_wr[3 * target_pts.get_n() + 2] = xyz[i + 2];
415  target_pts.inc_n();
416  }
417 
418  } // end for (unsigned int i = 0; ..
419 
420 #ifdef VERBOSE
421  int num_to_me = 0;
422  for( unsigned int i = 0; i < target_pts.get_n(); i++ )
423  if( target_pts.vi_rd[2 * i] == (int)my_rank ) num_to_me++;
424  printf( "rank: %u local points: %u, nb sent target pts: %u mappedPts: %u num to me: %d \n", my_rank, num_points,
425  target_pts.get_n(), mappedPts->get_n(), num_to_me );
426 #endif
427  // Perform scatter/gather, to gather points to source mesh procs
428  if( myPc )
429  {
430  ( myPc->proc_config().crystal_router() )->gs_transfer( 1, target_pts, 0 );
431 
432 #ifdef VERBOSE
433  num_to_me = 0;
434  for( unsigned int i = 0; i < target_pts.get_n(); i++ )
435  {
436  if( target_pts.vi_rd[2 * i] == (int)my_rank ) num_to_me++;
437  }
438  printf( "rank: %u after first gs nb received_pts: %u; num_from_me = %d\n", my_rank, target_pts.get_n(),
439  num_to_me );
440 #endif
441  // After scatter/gather:
442  // target_pts.set_n(# points local proc has to map);
443  // target_pts.vi_wr[2*i] = proc sending point i
444  // target_pts.vi_wr[2*i + 1] = index of point i on sending proc
445  // target_pts.vr_wr[3*i..3*i + 2] = xyz of point i
446  //
447  // Mapping builds the tuple list:
448  // source_pts.set_n(target_pts.get_n())
449  // source_pts.vi_wr[3*i] = target_pts.vi_wr[2*i] = sending proc
450  // source_pts.vi_wr[3*i + 1] = index of point i on sending proc
451  // source_pts.vi_wr[3*i + 2] = index of mapped point (-1 if not mapped)
452  //
453  // Also, mapping builds local tuple_list mappedPts:
454  // mappedPts->set_n( # mapped points );
455  // mappedPts->vul_wr[i] = local handle of mapped entity
456  // mappedPts->vr_wr[3*i..3*i + 2] = natural coordinates in mapped entity
457 
458  // Test target points against my elements
459  for( unsigned i = 0; i < target_pts.get_n(); i++ )
460  {
461  result = test_local_box( target_pts.vr_wr + 3 * i, target_pts.vi_rd[2 * i], target_pts.vi_rd[2 * i + 1], i,
462  point_located, rel_eps, abs_eps, &source_pts );
463  if( MB_SUCCESS != result ) return result;
464  }
465 
466  // No longer need target_pts
467  target_pts.reset();
468 #ifdef VERBOSE
469  printf( "rank: %u nb sent source pts: %u, mappedPts now: %u\n", my_rank, source_pts.get_n(),
470  mappedPts->get_n() );
471 #endif
472  // Send target points back to target procs
473  ( myPc->proc_config().crystal_router() )->gs_transfer( 1, source_pts, 0 );
474 
475 #ifdef VERBOSE
476  printf( "rank: %u nb received source pts: %u\n", my_rank, source_pts.get_n() );
477 #endif
478  }
479 
480  // Store proc/index tuples in targetPts, and/or pass back to application;
481  // the tuple this gets stored to looks like:
482  // tl.set_n(# mapped points);
483  // tl.vi_wr[3*i] = remote proc mapping point
484  // tl.vi_wr[3*i + 1] = local index of mapped point
485  // tl.vi_wr[3*i + 2] = remote index of mapped point
486  //
487  // Local index is mapped into either myRange, holding the handles of
488  // local mapped entities, or myXyz, holding locations of mapped pts
489 
490  // Store information about located points
491  TupleList* tl_tmp;
492  if( !store_local )
493  tl_tmp = tl;
494  else
495  {
496  targetPts = new TupleList();
497  tl_tmp = targetPts;
498  }
499 
500  tl_tmp->initialize( 3, 0, 0, 0, num_points );
501  tl_tmp->set_n( num_points ); // Automatically sets tl to write_enabled
502  // Initialize so we know afterwards how many pts weren't located
503  std::fill( tl_tmp->vi_wr, tl_tmp->vi_wr + 3 * num_points, -1 );
504 
505  unsigned int local_pts = 0;
506  for( unsigned int i = 0; i < source_pts.get_n(); i++ )
507  {
508  if( -1 != source_pts.vi_rd[3 * i + 2] )
509  { // Why bother sending message saying "i don't have the point" if it gets discarded?
510  int tgt_index = 3 * source_pts.vi_rd[3 * i + 1];
511  // Prefer always entities that are local, from the source_pts
512  // if a local entity was already found to contain the target point, skip
513  // tl_tmp->vi_wr[tgt_index] is -1 initially, but it could already be set with
514  // a remote processor
515  if( tl_tmp->vi_wr[tgt_index] != (int)my_rank )
516  {
517  tl_tmp->vi_wr[tgt_index] = source_pts.vi_rd[3 * i];
518  tl_tmp->vi_wr[tgt_index + 1] = source_pts.vi_rd[3 * i + 1];
519  tl_tmp->vi_wr[tgt_index + 2] = source_pts.vi_rd[3 * i + 2];
520  }
521  }
522  }
523 
524  // Count missing points
525  unsigned int missing_pts = 0;
526  for( unsigned int i = 0; i < num_points; i++ )
527  {
528  if( tl_tmp->vi_rd[3 * i + 1] == -1 )
529  {
530  missing_pts++;
531 #ifdef VERBOSE
532  printf( "missing point at index i: %d -> %15.10f %15.10f %15.10f\n", i, xyz[3 * i], xyz[3 * i + 1],
533  xyz[3 * i + 2] );
534 #endif
535  }
536  else if( tl_tmp->vi_rd[3 * i] == (int)my_rank )
537  local_pts++;
538  }
539  printf( "rank: %u point location: wanted %u got %u locally, %u remote, missing %u\n", my_rank, num_points,
540  local_pts, num_points - missing_pts - local_pts, missing_pts );
541  assert( 0 == missing_pts ); // Will likely break on curved geometries
542 
543  // No longer need source_pts
544  source_pts.reset();
545 
546  // Copy into tl if passed in and storing locally
547  if( tl && store_local )
548  {
549  tl->initialize( 3, 0, 0, 0, num_points );
550  tl->enableWriteAccess();
551  memcpy( tl->vi_wr, tl_tmp->vi_rd, 3 * tl_tmp->get_n() * sizeof( int ) );
552  tl->set_n( tl_tmp->get_n() );
553  tl->disableWriteAccess();
554  }
555 
556  tl_tmp->disableWriteAccess();
557 
558  // Done
559  return MB_SUCCESS;
560 }
561 
563  int from_proc,
564  int remote_index,
565  int /*index*/,
566  bool& point_located,
567  double rel_eps,
568  double abs_eps,
569  TupleList* tl )
570 {
571  std::vector< EntityHandle > entities;
572  std::vector< CartVect > nat_coords;
573  bool canWrite = false;
574  if( tl )
575  {
576  canWrite = tl->get_writeEnabled();
577  if( !canWrite )
578  {
579  tl->enableWriteAccess();
580  canWrite = true;
581  }
582  }
583 
584  if( rel_eps && !abs_eps )
585  {
586  // Relative epsilon given, translate to absolute epsilon using box dimensions
587  BoundBox box;
589  abs_eps = rel_eps * box.diagonal_length();
590  }
591 
592  ErrorCode result = nat_param( xyz, entities, nat_coords, abs_eps );
593  if( MB_SUCCESS != result ) return result;
594 
595  // If we didn't find any ents and we're looking locally, nothing more to do
596  if( entities.empty() )
597  {
598  if( tl->get_n() == tl->get_max() ) tl->resize( std::max( 10.0, 1.5 * tl->get_max() ) );
599 
600  tl->vi_wr[3 * tl->get_n()] = from_proc;
601  tl->vi_wr[3 * tl->get_n() + 1] = remote_index;
602  tl->vi_wr[3 * tl->get_n() + 2] = -1;
603  tl->inc_n();
604 
605  point_located = false;
606  return MB_SUCCESS;
607  }
608 
609  // Grow if we know we'll exceed size
610  if( mappedPts->get_n() + entities.size() >= mappedPts->get_max() )
611  mappedPts->resize( std::max( 10.0, 1.5 * mappedPts->get_max() ) );
612 
613  std::vector< EntityHandle >::iterator eit = entities.begin();
614  std::vector< CartVect >::iterator ncit = nat_coords.begin();
615 
617  for( ; eit != entities.end(); ++eit, ++ncit )
618  {
619  // Store in tuple mappedPts
620  mappedPts->vr_wr[3 * mappedPts->get_n()] = ( *ncit )[0];
621  mappedPts->vr_wr[3 * mappedPts->get_n() + 1] = ( *ncit )[1];
622  mappedPts->vr_wr[3 * mappedPts->get_n() + 2] = ( *ncit )[2];
623  mappedPts->vul_wr[mappedPts->get_n()] = *eit;
624  mappedPts->inc_n();
625 
626  // Also store local point, mapped point indices
627  if( tl->get_n() == tl->get_max() ) tl->resize( std::max( 10.0, 1.5 * tl->get_max() ) );
628 
629  // Store in tuple source_pts
630  tl->vi_wr[3 * tl->get_n()] = from_proc;
631  tl->vi_wr[3 * tl->get_n() + 1] = remote_index;
632  tl->vi_wr[3 * tl->get_n() + 2] = mappedPts->get_n() - 1;
633  tl->inc_n();
634  }
635 
636  point_located = true;
637 
638  if( tl && !canWrite ) tl->disableWriteAccess();
639 
640  return MB_SUCCESS;
641 }
642 
644  const std::string& interp_tag,
645  double* interp_vals,
646  TupleList* tl,
647  bool normalize )
648 {
649  Tag tag;
650  ErrorCode result;
651  if( _spectralSource )
652  {
653  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 << "\"" );
654  }
655  else
656  {
657  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 << "\"" );
658  }
659 
660  return interpolate( method, tag, interp_vals, tl, normalize );
661 }
662 
664  Tag* tags,
665  int* points_per_method,
666  int num_methods,
667  double* interp_vals,
668  TupleList* tl,
669  bool /* normalize */ )
670 {
671  // if (!((LINEAR_FE == method) || (CONSTANT == method)))
672  // return MB_FAILURE;
673 
674  // remote pts first
675  TupleList* tl_tmp = ( tl ? tl : targetPts );
676 
677  ErrorCode result = MB_SUCCESS;
678 
679  unsigned int pts_total = 0;
680  for( int i = 0; i < num_methods; i++ )
681  pts_total += points_per_method[i];
682 
683  // If tl was passed in non-NULL, just have those points, otherwise have targetPts plus
684  // locally mapped pts
685  if( pts_total != tl_tmp->get_n() ) return MB_FAILURE;
686 
687  TupleList tinterp;
688  tinterp.initialize( 5, 0, 0, 1, tl_tmp->get_n() );
689  int t = 0;
690  tinterp.enableWriteAccess();
691  for( int i = 0; i < num_methods; i++ )
692  {
693  for( int j = 0; j < points_per_method[i]; j++ )
694  {
695  tinterp.vi_wr[5 * t] = tl_tmp->vi_rd[3 * t];
696  tinterp.vi_wr[5 * t + 1] = tl_tmp->vi_rd[3 * t + 1];
697  tinterp.vi_wr[5 * t + 2] = tl_tmp->vi_rd[3 * t + 2];
698  tinterp.vi_wr[5 * t + 3] = methods[i];
699  tinterp.vi_wr[5 * t + 4] = i;
700  tinterp.vr_wr[t] = 0.0;
701  tinterp.inc_n();
702  t++;
703  }
704  }
705 
706  // Scatter/gather interpolation points
707  if( myPc )
708  {
709  ( myPc->proc_config().crystal_router() )->gs_transfer( 1, tinterp, 0 );
710 
711  // Perform interpolation on local source mesh; put results into
712  // tinterp.vr_wr[i]
713 
715  for( unsigned int i = 0; i < tinterp.get_n(); i++ )
716  {
717  int mindex = tinterp.vi_rd[5 * i + 2];
718  Method method = (Method)tinterp.vi_rd[5 * i + 3];
719  Tag tag = tags[tinterp.vi_rd[5 * i + 4]];
720 
721  result = MB_FAILURE;
722  if( LINEAR_FE == method || QUADRATIC_FE == method || SPHERICAL == method )
723  {
724  result = interp_field( mappedPts->vul_rd[mindex], CartVect( mappedPts->vr_wr + 3 * mindex ), tag,
725  tinterp.vr_wr[i] );
726  }
727  else if( CONSTANT == method )
728  {
729  result = constant_interp( mappedPts->vul_rd[mindex], tag, tinterp.vr_wr[i] );
730  }
731 
732  if( MB_SUCCESS != result ) return result;
733  }
734 
735  // Scatter/gather interpolation data
736  myPc->proc_config().crystal_router()->gs_transfer( 1, tinterp, 0 );
737  }
738 
739  // Copy the interpolated field as a unit
740  for( unsigned int i = 0; i < tinterp.get_n(); i++ )
741  interp_vals[tinterp.vi_rd[5 * i + 1]] = tinterp.vr_rd[i];
742 
743  // Done
744  return MB_SUCCESS;
745 }
746 
748  std::vector< EntityHandle >& entities,
749  std::vector< CartVect >& nat_coords,
750  double epsilon )
751 {
752  if( !myTree ) return MB_FAILURE;
753 
754  AdaptiveKDTreeIter treeiter;
755  ErrorCode result = myTree->get_tree_iterator( localRoot, treeiter );
756  if( MB_SUCCESS != result )
757  {
758  std::cout << "Problems getting iterator" << std::endl;
759  return result;
760  }
761 
762  EntityHandle closest_leaf;
763  if( epsilon )
764  {
765  std::vector< double > dists;
766  std::vector< EntityHandle > leaves;
767 
768  // Two tolerances
769  result = myTree->distance_search( xyz, epsilon, leaves,
770  /*iter_tol*/ epsilon,
771  /*inside_tol*/ 10 * epsilon, &dists, NULL, &localRoot );
772  if( leaves.empty() )
773  // Not found returns success here, with empty list, just like case with no epsilon
774  return MB_SUCCESS;
775  // Get closest leaf
776  double min_dist = *dists.begin();
777  closest_leaf = *leaves.begin();
778  std::vector< EntityHandle >::iterator vit = leaves.begin() + 1;
779  std::vector< double >::iterator dit = dists.begin() + 1;
780  for( ; vit != leaves.end() && min_dist; ++vit, ++dit )
781  {
782  if( *dit < min_dist )
783  {
784  min_dist = *dit;
785  closest_leaf = *vit;
786  }
787  }
788  }
789  else
790  {
791  result = myTree->point_search( xyz, treeiter, 1.0e-10, 1.0e-6, NULL, &localRoot );
792  if( MB_ENTITY_NOT_FOUND == result ) // Point is outside of myTree's bounding box
793  return MB_SUCCESS;
794  else if( MB_SUCCESS != result )
795  {
796  std::cout << "Problems getting leaf \n";
797  return result;
798  }
799  closest_leaf = treeiter.handle();
800  }
801 
802  // Find natural coordinates of point in element(s) in that leaf
803  CartVect tmp_nat_coords;
804  Range range_leaf;
805  result = mbImpl->get_entities_by_dimension( closest_leaf, max_dim, range_leaf, false );
806  if( MB_SUCCESS != result ) std::cout << "Problem getting leaf in a range" << std::endl;
807 
808  // Loop over the range_leaf
809  for( Range::iterator iter = range_leaf.begin(); iter != range_leaf.end(); ++iter )
810  {
811  // Test to find out in which entity the point is
812  // Get the EntityType and create the appropriate Element::Map subtype
813  // If spectral, do not need coordinates, just the GL points
814  EntityType etype = mbImpl->type_from_handle( *iter );
815  if( NULL != this->_spectralSource && MBHEX == etype )
816  {
817  EntityHandle eh = *iter;
818  const double* xval;
819  const double* yval;
820  const double* zval;
821  ErrorCode rval = mbImpl->tag_get_by_ptr( _xm1Tag, &eh, 1, (const void**)&xval );
822  if( moab::MB_SUCCESS != rval )
823  {
824  std::cout << "Can't get xm1 values \n";
825  return MB_FAILURE;
826  }
827  rval = mbImpl->tag_get_by_ptr( _ym1Tag, &eh, 1, (const void**)&yval );
828  if( moab::MB_SUCCESS != rval )
829  {
830  std::cout << "Can't get ym1 values \n";
831  return MB_FAILURE;
832  }
833  rval = mbImpl->tag_get_by_ptr( _zm1Tag, &eh, 1, (const void**)&zval );
834  if( moab::MB_SUCCESS != rval )
835  {
836  std::cout << "Can't get zm1 values \n";
837  return MB_FAILURE;
838  }
840 
841  spcHex->set_gl_points( (double*)xval, (double*)yval, (double*)zval );
842  try
843  {
844  tmp_nat_coords = spcHex->ievaluate( CartVect( xyz ), epsilon ); // introduce
845  bool inside = spcHex->inside_nat_space( CartVect( tmp_nat_coords ), epsilon );
846  if( !inside )
847  {
848 #ifdef VERBOSE
849  std::cout << "point " << xyz[0] << " " << xyz[1] << " " << xyz[2]
850  << " is not converging inside hex " << mbImpl->id_from_handle( eh ) << "\n";
851 #endif
852  continue; // It is possible that the point is outside, so it will not converge
853  }
854  }
856  {
857  continue;
858  }
859  }
860  else
861  {
862  const EntityHandle* connect;
863  int num_connect;
864 
865  // Get connectivity
866  result = mbImpl->get_connectivity( *iter, connect, num_connect, true );
867  if( MB_SUCCESS != result ) return result;
868 
869  // Get coordinates of the vertices
870  std::vector< CartVect > coords_vert( num_connect );
871  result = mbImpl->get_coords( connect, num_connect, &( coords_vert[0][0] ) );
872  if( MB_SUCCESS != result )
873  {
874  std::cout << "Problems getting coordinates of vertices\n";
875  return result;
876  }
877  CartVect pos( xyz );
878  if( MBHEX == etype )
879  {
880  if( 8 == num_connect )
881  {
882  Element::LinearHex hexmap( coords_vert );
883  if( !hexmap.inside_box( pos, epsilon ) ) continue;
884  try
885  {
886  tmp_nat_coords = hexmap.ievaluate( pos, epsilon );
887  bool inside = hexmap.inside_nat_space( tmp_nat_coords, epsilon );
888  if( !inside ) continue;
889  }
891  {
892  continue;
893  }
894  }
895  else if( 27 == num_connect )
896  {
897  Element::QuadraticHex hexmap( coords_vert );
898  if( !hexmap.inside_box( pos, epsilon ) ) continue;
899  try
900  {
901  tmp_nat_coords = hexmap.ievaluate( pos, epsilon );
902  bool inside = hexmap.inside_nat_space( tmp_nat_coords, epsilon );
903  if( !inside ) continue;
904  }
906  {
907  continue;
908  }
909  }
910  else // TODO this case not treated yet, no interpolation
911  continue;
912  }
913  else if( MBTET == etype )
914  {
915  Element::LinearTet tetmap( coords_vert );
916  // This is just a linear solve; unless degenerate, will not except
917  tmp_nat_coords = tetmap.ievaluate( pos );
918  bool inside = tetmap.inside_nat_space( tmp_nat_coords, epsilon );
919  if( !inside ) continue;
920  }
921  else if( MBQUAD == etype && spherical )
922  {
923  Element::SphericalQuad sphermap( coords_vert );
924  /* skip box test, because it can filter out good elements with high curvature
925  * if (!sphermap.inside_box(pos, epsilon))
926  continue;*/
927  try
928  {
929  tmp_nat_coords = sphermap.ievaluate( pos, epsilon );
930  bool inside = sphermap.inside_nat_space( tmp_nat_coords, epsilon );
931  if( !inside ) continue;
932  }
934  {
935  continue;
936  }
937  }
938  else if( MBTRI == etype && spherical )
939  {
940  Element::SphericalTri sphermap( coords_vert );
941  /* skip box test, because it can filter out good elements with high curvature
942  * if (!sphermap.inside_box(pos, epsilon))
943  continue;*/
944  try
945  {
946  tmp_nat_coords = sphermap.ievaluate( pos, epsilon );
947  bool inside = sphermap.inside_nat_space( tmp_nat_coords, epsilon );
948  if( !inside ) continue;
949  }
951  {
952  continue;
953  }
954  }
955 
956  else if( MBQUAD == etype )
957  {
958  Element::LinearQuad quadmap( coords_vert );
959  if( !quadmap.inside_box( pos, epsilon ) ) continue;
960  try
961  {
962  tmp_nat_coords = quadmap.ievaluate( pos, epsilon );
963  bool inside = quadmap.inside_nat_space( tmp_nat_coords, epsilon );
964  if( !inside ) continue;
965  }
967  {
968  continue;
969  }
970  if( !quadmap.inside_nat_space( tmp_nat_coords, epsilon ) ) continue;
971  }
972  /*
973  else if (etype == MBTRI){
974  Element::LinearTri trimap(coords_vert);
975  if (!trimap.inside_box( pos, epsilon))
976  continue;
977  try {
978  tmp_nat_coords = trimap.ievaluate(pos, epsilon);
979  bool inside = trimap.inside_nat_space(tmp_nat_coords, epsilon);
980  if (!inside) continue;
981  }
982  catch (Element::Map::EvaluationError) {
983  continue;
984  }
985  if (!trimap.inside_nat_space(tmp_nat_coords, epsilon))
986  continue;
987  }
988  */
989  else if( etype == MBEDGE )
990  {
991  Element::LinearEdge edgemap( coords_vert );
992  try
993  {
994  tmp_nat_coords = edgemap.ievaluate( CartVect( xyz ), epsilon );
995  }
997  {
998  continue;
999  }
1000  if( !edgemap.inside_nat_space( tmp_nat_coords, epsilon ) ) continue;
1001  }
1002  else
1003  {
1004  std::cout << "Entity not Hex/Tet/Quad/Tri/Edge. Please verify." << std::endl;
1005  continue;
1006  }
1007  }
1008  // If we get here then we've found the coordinates.
1009  // Save them and the entity and return success.
1010  entities.push_back( *iter );
1011  nat_coords.push_back( tmp_nat_coords );
1012  return MB_SUCCESS;
1013  }
1014 
1015  // Didn't find any elements containing the point
1016  return MB_SUCCESS;
1017 }
1018 
1019 ErrorCode Coupler::interp_field( EntityHandle elem, CartVect nat_coord, Tag tag, double& field )
1020 {
1021  if( _spectralSource )
1022  {
1023  // Get tag values at the GL points for some field (Tag)
1024  const double* vx;
1025  ErrorCode rval = mbImpl->tag_get_by_ptr( tag, &elem, 1, (const void**)&vx );
1026  if( moab::MB_SUCCESS != rval )
1027  {
1028  std::cout << "Can't get field values for the tag \n";
1029  return MB_FAILURE;
1030  }
1032  field = spcHex->evaluate_scalar_field( nat_coord, vx );
1033  }
1034  else
1035  {
1036  double vfields[27]; // Will work for linear hex, quadratic hex or Tets
1037  moab::Element::Map* elemMap = NULL;
1038  int num_verts = 0;
1039  // Get the EntityType
1040  // Get the tag values at the vertices
1041  const EntityHandle* connect;
1042  int num_connect;
1043  ErrorCode result = mbImpl->get_connectivity( elem, connect, num_connect );
1044  if( MB_SUCCESS != result ) return result;
1045  EntityType etype = mbImpl->type_from_handle( elem );
1046  if( MBHEX == etype )
1047  {
1048  if( 8 == num_connect )
1049  {
1050  elemMap = new moab::Element::LinearHex();
1051  num_verts = 8;
1052  }
1053  else
1054  { /* (MBHEX == etype && 27 == num_connect) */
1055  elemMap = new moab::Element::QuadraticHex();
1056  num_verts = 27;
1057  }
1058  }
1059  else if( MBTET == etype )
1060  {
1061  elemMap = new moab::Element::LinearTet();
1062  num_verts = 4;
1063  }
1064  else if( MBQUAD == etype )
1065  {
1066  elemMap = new moab::Element::LinearQuad();
1067  num_verts = 4;
1068  }
1069  else if( MBTRI == etype )
1070  {
1071  elemMap = new moab::Element::LinearTri();
1072  num_verts = 3;
1073  }
1074  else
1075  return MB_FAILURE;
1076 
1077  result = mbImpl->tag_get_data( tag, connect, std::min( num_verts, num_connect ), vfields );
1078  if( MB_SUCCESS != result )
1079  {
1080  delete elemMap;
1081  return result;
1082  }
1083 
1084  // Function for the interpolation
1085  field = 0;
1086 
1087  // Check the number of vertices
1088  assert( num_connect >= num_verts );
1089 
1090  // Calculate the field
1091  try
1092  {
1093  field = elemMap->evaluate_scalar_field( nat_coord, vfields );
1094  }
1096  {
1097  delete elemMap;
1098  return MB_FAILURE;
1099  }
1100 
1101  delete elemMap;
1102  }
1103 
1104  return MB_SUCCESS;
1105 }
1106 
1107 // Simplest "interpolation" for element-based source fields. Set the value of the field
1108 // at the target point to that of the field in the source element it lies in.
1110 {
1111  double tempField;
1112 
1113  // Get the tag values at the vertices
1114  ErrorCode result = mbImpl->tag_get_data( tag, &elem, 1, &tempField );
1115  if( MB_SUCCESS != result ) return result;
1116 
1117  field = tempField;
1118 
1119  return MB_SUCCESS;
1120 }
1121 
1122 // Normalize a field over the entire mesh represented by the root_set.
1124  const char* norm_tag,
1125  Coupler::IntegType integ_type,
1126  int num_integ_pts )
1127 {
1128  ErrorCode err;
1129 
1130  // SLAVE START ****************************************************************
1131  // Search for entities based on tag_handles and tag_values
1132  std::vector< std::vector< EntityHandle > > entity_sets;
1133  std::vector< std::vector< EntityHandle > > entity_groups;
1134 
1135  // put the root_set into entity_sets
1136  std::vector< EntityHandle > ent_set;
1137  ent_set.push_back( root_set );
1138  entity_sets.push_back( ent_set );
1139 
1140  // get all entities from root_set and put into entity_groups
1141  std::vector< EntityHandle > entities;
1142  err = mbImpl->get_entities_by_handle( root_set, entities, true );ERRORR( "Failed to get entities in root_set.", err );
1143 
1144  entity_groups.push_back( entities );
1145 
1146  // Call do_normalization() to continue common normalization processing
1147  err = do_normalization( norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts );ERRORR( "Failure in do_normalization().", err );
1148  // SLAVE END ****************************************************************
1149 
1150  return err;
1151 }
1152 
1153 // Normalize a field over the subset of entities identified by the tags and values passed
1155  const char* norm_tag,
1156  const char** tag_names,
1157  int num_tags,
1158  const char** tag_values,
1159  Coupler::IntegType integ_type,
1160  int num_integ_pts )
1161 {
1162  moab::ErrorCode err;
1163  std::vector< Tag > tag_handles;
1164 
1165  // Lookup tag handles from tag names
1166  for( int t = 0; t < num_tags; t++ )
1167  {
1168  // get tag handle & size
1169  Tag th;
1170  err = mbImpl->tag_get_handle( tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY );ERRORR( "Failed to get tag handle.", err );
1171  tag_handles.push_back( th );
1172  }
1173 
1174  return normalize_subset( root_set, norm_tag, &tag_handles[0], num_tags, tag_values, integ_type, num_integ_pts );
1175 }
1176 
1178  const char* norm_tag,
1179  Tag* tag_handles,
1180  int num_tags,
1181  const char** tag_values,
1182  Coupler::IntegType integ_type,
1183  int num_integ_pts )
1184 {
1185  ErrorCode err;
1186 
1187  // SLAVE START ****************************************************************
1188  // Search for entities based on tag_handles and tag_values
1189  std::vector< std::vector< EntityHandle > > entity_sets;
1190  std::vector< std::vector< EntityHandle > > entity_groups;
1191 
1192  err = get_matching_entities( root_set, tag_handles, tag_values, num_tags, &entity_sets, &entity_groups );ERRORR( "Failed to get matching entities.", err );
1193 
1194  // Call do_normalization() to continue common normalization processing
1195  err = do_normalization( norm_tag, entity_sets, entity_groups, integ_type, num_integ_pts );ERRORR( "Failure in do_normalization().", err );
1196  // SLAVE END ****************************************************************
1197 
1198  return err;
1199 }
1200 
1201 ErrorCode Coupler::do_normalization( const char* norm_tag,
1202  std::vector< std::vector< EntityHandle > >& entity_sets,
1203  std::vector< std::vector< EntityHandle > >& entity_groups,
1204  Coupler::IntegType integ_type,
1205  int num_integ_pts )
1206 {
1207  // SLAVE START ****************************************************************
1208  ErrorCode err;
1209  int ierr = 0;
1210 
1211  // Setup data for parallel computing
1212  int nprocs, rank;
1213  ierr = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
1214  ERRORMPI( "Getting number of procs failed.", ierr );
1215  ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
1216  ERRORMPI( "Getting rank failed.", ierr );
1217 
1218  // Get the integrated field value for each group(vector) of entities.
1219  // If no entities are in a group then a zero will be put in the list
1220  // of return values.
1221  unsigned int num_ent_grps = entity_groups.size();
1222  std::vector< double > integ_vals( num_ent_grps );
1223 
1224  err = get_group_integ_vals( entity_groups, integ_vals, norm_tag, num_integ_pts, integ_type );ERRORR( "Failed to get integrated field values for groups in mesh.", err );
1225  // SLAVE END ****************************************************************
1226 
1227  // SLAVE/MASTER START #########################################################
1228  // Send list of integrated values back to master proc. The ordering of the
1229  // values will match the ordering of the entity groups (i.e. vector of vectors)
1230  // sent from master to slaves earlier. The values for each entity group will
1231  // be summed during the transfer.
1232  std::vector< double > sum_integ_vals( num_ent_grps );
1233 
1234  if( nprocs > 1 )
1235  {
1236  // If parallel then send the values back to the master.
1237  ierr = MPI_Reduce( &integ_vals[0], &sum_integ_vals[0], num_ent_grps, MPI_DOUBLE, MPI_SUM, MASTER_PROC,
1238  myPc->proc_config().proc_comm() );
1239  ERRORMPI( "Transfer and reduction of integrated values failed.", ierr );
1240  }
1241  else
1242  {
1243  // Otherwise just copy the vector
1244  sum_integ_vals = integ_vals;
1245  }
1246  // SLAVE/MASTER END #########################################################
1247 
1248  // MASTER START ***************************************************************
1249  // Calculate the normalization factor for each group by taking the
1250  // inverse of each integrated field value. Put the normalization factor
1251  // for each group back into the list in the same order.
1252  for( unsigned int i = 0; i < num_ent_grps; i++ )
1253  {
1254  double val = sum_integ_vals[i];
1255  if( fabs( val ) > 1e-8 )
1256  sum_integ_vals[i] = 1.0 / val;
1257  else
1258  {
1259  sum_integ_vals[i] = 0.0; /* VSM: not sure what we should do here ? */
1260  /* commenting out error below since if integral(value)=0.0, then normalization
1261  is probably unnecessary to start with ? */
1262  /* ERRORR("Integrating an invalid field -- integral("<<norm_tag<<") = "<<val<<".", err);
1263  */
1264  }
1265  }
1266  // MASTER END ***************************************************************
1267 
1268  // MASTER/SLAVE START #########################################################
1269  if( nprocs > 1 )
1270  {
1271  // If parallel then broadcast the normalization factors to the procs.
1272  ierr = MPI_Bcast( &sum_integ_vals[0], num_ent_grps, MPI_DOUBLE, MASTER_PROC, myPc->proc_config().proc_comm() );
1273  ERRORMPI( "Broadcast of normalization factors failed.", ierr );
1274  }
1275  // MASTER/SLAVE END #########################################################
1276 
1277  // SLAVE START ****************************************************************
1278  // Save the normalization factors to a new tag with name of norm_tag's value
1279  // and the string "_normF" appended. This new tag will be created on the entity
1280  // set that contains all of the entities from a group.
1281 
1282  err = apply_group_norm_factor( entity_sets, sum_integ_vals, norm_tag, integ_type );ERRORR( "Failed to set the normalization factor for groups in mesh.", err );
1283  // SLAVE END ****************************************************************
1284 
1285  return err;
1286 }
1287 
1288 // Functions supporting the subset normalization function
1289 
1290 // Retrieve groups of entities matching tags and values if present
1292  const char** tag_names,
1293  const char** tag_values,
1294  int num_tags,
1295  std::vector< std::vector< EntityHandle > >* entity_sets,
1296  std::vector< std::vector< EntityHandle > >* entity_groups )
1297 {
1298  ErrorCode err;
1299  std::vector< Tag > tag_handles;
1300 
1301  for( int t = 0; t < num_tags; t++ )
1302  {
1303  // Get tag handle & size
1304  Tag th;
1305  err = mbImpl->tag_get_handle( tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY );ERRORR( "Failed to get tag handle.", err );
1306  tag_handles.push_back( th );
1307  }
1308 
1309  return get_matching_entities( root_set, &tag_handles[0], tag_values, num_tags, entity_sets, entity_groups );
1310 }
1311 
1312 // Retrieve groups of entities matching tags and values if present
1314  Tag* tag_handles,
1315  const char** tag_values,
1316  int num_tags,
1317  std::vector< std::vector< EntityHandle > >* entity_sets,
1318  std::vector< std::vector< EntityHandle > >* entity_groups )
1319 {
1320  // SLAVE START ****************************************************************
1321 
1322  // Setup data for parallel computing
1323  ErrorCode err;
1324  int ierr = 0;
1325  int nprocs, rank;
1326  ierr = MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
1327  ERRORMPI( "Getting number of procs failed.", ierr );
1328  ierr = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
1329  ERRORMPI( "Getting rank failed.", ierr );
1330 
1331  Range ent_sets;
1332  err =
1333  mbImpl->get_entities_by_type_and_tag( root_set, moab::MBENTITYSET, tag_handles, (const void* const*)tag_values,
1334  num_tags, ent_sets, Interface::INTERSECT, false );ERRORR( "Core::get_entities_by_type_and_tag failed.", err );
1335 
1336  TupleList* tag_list = NULL;
1337  err = create_tuples( ent_sets, tag_handles, num_tags, &tag_list );ERRORR( "Failed to create tuples from entity sets.", err );
1338 
1339  // Free up range
1340  ent_sets.clear();
1341  // SLAVE END ****************************************************************
1342 
1343  // If we are running in a multi-proc session then send tuple list back to master
1344  // proc for consolidation. Otherwise just copy the pointer to the tuple_list.
1345  TupleList* cons_tuples;
1346  if( nprocs > 1 )
1347  {
1348  // SLAVE/MASTER START #########################################################
1349 
1350  // Pack the tuple_list in a buffer.
1351  uint* tuple_buf;
1352  int tuple_buf_len;
1353  tuple_buf_len = pack_tuples( tag_list, (void**)&tuple_buf );
1354 
1355  // Free tag_list here as its not used again if nprocs > 1
1356  tag_list->reset();
1357 
1358  // Send back the buffer sizes to the master proc
1359  int* recv_cnts = (int*)malloc( nprocs * sizeof( int ) );
1360  int* offsets = (int*)malloc( nprocs * sizeof( int ) );
1361  uint* all_tuples_buf = NULL;
1362 
1363  MPI_Gather( &tuple_buf_len, 1, MPI_INT, recv_cnts, 1, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
1364  ERRORMPI( "Gathering buffer sizes failed.", err );
1365 
1366  // Allocate a buffer large enough for all the data
1367  if( rank == MASTER_PROC )
1368  {
1369  int all_tuples_len = recv_cnts[0];
1370  offsets[0] = 0;
1371  for( int i = 1; i < nprocs; i++ )
1372  {
1373  offsets[i] = offsets[i - 1] + recv_cnts[i - 1];
1374  all_tuples_len += recv_cnts[i];
1375  }
1376 
1377  all_tuples_buf = (uint*)malloc( all_tuples_len * sizeof( uint ) );
1378  }
1379 
1380  // Send all buffers to the master proc for consolidation
1381  MPI_Gatherv( (void*)tuple_buf, tuple_buf_len, MPI_INT, (void*)all_tuples_buf, recv_cnts, offsets, MPI_INT,
1383  ERRORMPI( "Gathering tuple_lists failed.", err );
1384  free( tuple_buf ); // malloc'd in pack_tuples
1385 
1386  if( rank == MASTER_PROC )
1387  {
1388  // Unpack the tuple_list from the buffer.
1389  TupleList** tl_array = (TupleList**)malloc( nprocs * sizeof( TupleList* ) );
1390  for( int i = 0; i < nprocs; i++ )
1391  unpack_tuples( (void*)&all_tuples_buf[offsets[i]], &tl_array[i] );
1392 
1393  // Free all_tuples_buf here as it is only allocated on the MASTER_PROC
1394  free( all_tuples_buf );
1395  // SLAVE/MASTER END #########################################################
1396 
1397  // MASTER START ***************************************************************
1398  // Consolidate all tuple_lists into one tuple_list with no duplicates.
1399  err = consolidate_tuples( tl_array, nprocs, &cons_tuples );ERRORR( "Failed to consolidate tuples.", err );
1400 
1401  for( int i = 0; i < nprocs; i++ )
1402  tl_array[i]->reset();
1403  free( tl_array );
1404  // MASTER END ***************************************************************
1405  }
1406 
1407  // Free offsets and recv_cnts as they are allocated on all procs
1408  free( offsets );
1409  free( recv_cnts );
1410 
1411  // MASTER/SLAVE START #########################################################
1412  // Broadcast condensed tuple list back to all procs.
1413  uint* ctl_buf;
1414  int ctl_buf_sz;
1415  if( rank == MASTER_PROC ) ctl_buf_sz = pack_tuples( cons_tuples, (void**)&ctl_buf );
1416 
1417  // Send buffer size
1418  ierr = MPI_Bcast( &ctl_buf_sz, 1, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
1419  ERRORMPI( "Broadcasting tuple_list size failed.", ierr );
1420 
1421  // Allocate a buffer in the other procs
1422  if( rank != MASTER_PROC ) ctl_buf = (uint*)malloc( ctl_buf_sz * sizeof( uint ) );
1423 
1424  ierr = MPI_Bcast( (void*)ctl_buf, ctl_buf_sz, MPI_INT, MASTER_PROC, myPc->proc_config().proc_comm() );
1425  ERRORMPI( "Broadcasting tuple_list failed.", ierr );
1426 
1427  if( rank != MASTER_PROC ) unpack_tuples( ctl_buf, &cons_tuples );
1428  free( ctl_buf );
1429  // MASTER/SLAVE END #########################################################
1430  }
1431  else
1432  cons_tuples = tag_list;
1433 
1434  // SLAVE START ****************************************************************
1435  // Loop over the tuple list getting the entities with the tags in the tuple_list entry
1436  uint mi, ml, mul, mr;
1437  cons_tuples->getTupleSize( mi, ml, mul, mr );
1438 
1439  for( unsigned int i = 0; i < cons_tuples->get_n(); i++ )
1440  {
1441  // Get Entity Sets that match the tags and values.
1442 
1443  // Convert the data in the tuple_list to an array of pointers to the data
1444  // in the tuple_list as that is what the iMesh API call is expecting.
1445  int** vals = (int**)malloc( mi * sizeof( int* ) );
1446  for( unsigned int j = 0; j < mi; j++ )
1447  vals[j] = (int*)&( cons_tuples->vi_rd[( i * mi ) + j] );
1448 
1449  // Get entities recursively based on type and tag data
1450  err = mbImpl->get_entities_by_type_and_tag( root_set, moab::MBENTITYSET, tag_handles, (const void* const*)vals,
1451  mi, ent_sets, Interface::INTERSECT, false );ERRORR( "Core::get_entities_by_type_and_tag failed.", err );
1452  if( debug ) std::cout << "ent_sets_size=" << ent_sets.size() << std::endl;
1453 
1454  // Free up the array of pointers
1455  free( vals );
1456 
1457  // Loop over the entity sets and then free the memory for ent_sets.
1458  std::vector< EntityHandle > ent_set_hdls;
1459  std::vector< EntityHandle > ent_hdls;
1460  for( unsigned int j = 0; j < ent_sets.size(); j++ )
1461  {
1462  // Save the entity set
1463  ent_set_hdls.push_back( ent_sets[j] );
1464 
1465  // Get all entities for the entity set
1466  Range ents;
1467 
1468  /* VSM: do we need to filter out entity sets ? */
1469  err = mbImpl->get_entities_by_handle( ent_sets[j], ents, false );ERRORR( "Core::get_entities_by_handle failed.", err );
1470  if( debug ) std::cout << "ents_size=" << ents.size() << std::endl;
1471 
1472  // Save all of the entities from the entity set and free the memory for ents.
1473  for( unsigned int k = 0; k < ents.size(); k++ )
1474  {
1475  ent_hdls.push_back( ents[k] );
1476  }
1477  ents.clear();
1478  if( debug ) std::cout << "ent_hdls.size=" << ent_hdls.size() << std::endl;
1479  }
1480 
1481  // Free the entity set list for next tuple iteration.
1482  ent_sets.clear();
1483 
1484  // Push ent_set_hdls onto entity_sets, ent_hdls onto entity_groups
1485  // and clear both ent_set_hdls and ent_hdls.
1486  entity_sets->push_back( ent_set_hdls );
1487  ent_set_hdls.clear();
1488  entity_groups->push_back( ent_hdls );
1489  ent_hdls.clear();
1490  if( debug )
1491  std::cout << "entity_sets->size=" << entity_sets->size()
1492  << ", entity_groups->size=" << entity_groups->size() << std::endl;
1493  }
1494 
1495  cons_tuples->reset();
1496  // SLAVE END ****************************************************************
1497 
1498  return err;
1499 }
1500 
1501 // Return a tuple_list containing tag values for each Entity Set
1502 // The tuple_list will have a column for each tag and a row for each
1503 // Entity Set. It is assumed all of the tags are integer tags.
1505  const char** tag_names,
1506  unsigned int num_tags,
1507  TupleList** tuple_list )
1508 {
1509  ErrorCode err;
1510  std::vector< Tag > tag_handles;
1511 
1512  for( unsigned int t = 0; t < num_tags; t++ )
1513  {
1514  // Get tag handle & size
1515  Tag th;
1516  err = mbImpl->tag_get_handle( tag_names[t], 1, moab::MB_TYPE_DOUBLE, th, moab::MB_TAG_ANY );ERRORR( "Failed to get tag handle.", err );
1517  tag_handles.push_back( th );
1518  }
1519 
1520  return create_tuples( ent_sets, &tag_handles[0], num_tags, tuple_list );
1521 }
1522 
1523 // Return a tuple_list containing tag values for each Entity Set
1524 // The tuple_list will have a column for each tag and a row for each
1525 // Entity Set. It is assumed all of the tags are integer tags.
1526 ErrorCode Coupler::create_tuples( Range& ent_sets, Tag* tag_handles, unsigned int num_tags, TupleList** tuples )
1527 {
1528  // ASSUMPTION: All tags are of type integer. This may need to be expanded in future.
1529  ErrorCode err;
1530 
1531  // Allocate a tuple_list for the number of entity sets passed in
1532  TupleList* tag_tuples = new TupleList( num_tags, 0, 0, 0, (int)ent_sets.size() );
1533  // tag_tuples->initialize(num_tags, 0, 0, 0, num_sets);
1534  uint mi, ml, mul, mr;
1535  tag_tuples->getTupleSize( mi, ml, mul, mr );
1536  tag_tuples->enableWriteAccess();
1537 
1538  if( mi == 0 ) ERRORR( "Failed to initialize tuple_list.", MB_FAILURE );
1539 
1540  // Loop over the filtered entity sets retrieving each matching tag value one by one.
1541  int val;
1542  for( unsigned int i = 0; i < ent_sets.size(); i++ )
1543  {
1544  for( unsigned int j = 0; j < num_tags; j++ )
1545  {
1546  EntityHandle set_handle = ent_sets[i];
1547  err = mbImpl->tag_get_data( tag_handles[j], &set_handle, 1, &val );ERRORR( "Failed to get integer tag data.", err );
1548  tag_tuples->vi_wr[i * mi + j] = val;
1549  }
1550 
1551  // If we get here there was no error so increment n in the tuple_list
1552  tag_tuples->inc_n();
1553  }
1554  tag_tuples->disableWriteAccess();
1555  *tuples = tag_tuples;
1556 
1557  return MB_SUCCESS;
1558 }
1559 
1560 // Consolidate tuple_lists into one list with no duplicates
1561 ErrorCode Coupler::consolidate_tuples( TupleList** all_tuples, unsigned int num_tuples, TupleList** unique_tuples )
1562 {
1563  int total_rcv_tuples = 0;
1564  int offset = 0, copysz = 0;
1565  unsigned num_tags = 0;
1566 
1567  uint ml, mul, mr;
1568  uint* mi = (uint*)malloc( sizeof( uint ) * num_tuples );
1569 
1570  for( unsigned int i = 0; i < num_tuples; i++ )
1571  {
1572  all_tuples[i]->getTupleSize( mi[i], ml, mul, mr );
1573  }
1574 
1575  for( unsigned int i = 0; i < num_tuples; i++ )
1576  {
1577  if( all_tuples[i] != NULL )
1578  {
1579  total_rcv_tuples += all_tuples[i]->get_n();
1580  num_tags = mi[i];
1581  }
1582  }
1583  const unsigned int_size = sizeof( sint );
1584  const unsigned int_width = num_tags * int_size;
1585 
1586  // Get the total size of all of the tuple_lists in all_tuples.
1587  for( unsigned int i = 0; i < num_tuples; i++ )
1588  {
1589  if( all_tuples[i] != NULL ) total_rcv_tuples += all_tuples[i]->get_n();
1590  }
1591 
1592  // Copy the tuple_lists into a single tuple_list.
1593  TupleList* all_tuples_list = new TupleList( num_tags, 0, 0, 0, total_rcv_tuples );
1594  all_tuples_list->enableWriteAccess();
1595  // all_tuples_list->initialize(num_tags, 0, 0, 0, total_rcv_tuples);
1596  for( unsigned int i = 0; i < num_tuples; i++ )
1597  {
1598  if( all_tuples[i] != NULL )
1599  {
1600  copysz = all_tuples[i]->get_n() * int_width;
1601  memcpy( all_tuples_list->vi_wr + offset, all_tuples[i]->vi_rd, copysz );
1602  offset = offset + ( all_tuples[i]->get_n() * mi[i] );
1603  all_tuples_list->set_n( all_tuples_list->get_n() + all_tuples[i]->get_n() );
1604  }
1605  }
1606 
1607  // Sort the new tuple_list. Use a radix type sort, starting with the last (or least significant)
1608  // tag column in the vi array and working towards the first (or most significant) tag column.
1609  TupleList::buffer sort_buffer;
1610  sort_buffer.buffer_init( 2 * total_rcv_tuples * int_width );
1611  for( int i = num_tags - 1; i >= 0; i-- )
1612  {
1613  all_tuples_list->sort( i, &sort_buffer );
1614  }
1615 
1616  // Cycle through the sorted list eliminating duplicates.
1617  // Keep counters to the current end of the tuple_list (w/out dups) and the last tuple examined.
1618  unsigned int end_idx = 0, last_idx = 1;
1619  while( last_idx < all_tuples_list->get_n() )
1620  {
1621  if( memcmp( all_tuples_list->vi_rd + ( end_idx * num_tags ), all_tuples_list->vi_rd + ( last_idx * num_tags ),
1622  int_width ) == 0 )
1623  {
1624  // Values equal - skip
1625  last_idx += 1;
1626  }
1627  else
1628  {
1629  // Values different - copy
1630  // Move up the end index
1631  end_idx += 1;
1632  memcpy( all_tuples_list->vi_wr + ( end_idx * num_tags ), all_tuples_list->vi_rd + ( last_idx * num_tags ),
1633  int_width );
1634  last_idx += 1;
1635  }
1636  }
1637  // Update the count in all_tuples_list
1638  all_tuples_list->set_n( end_idx + 1 );
1639 
1640  // Resize the tuple_list
1641  all_tuples_list->resize( all_tuples_list->get_n() );
1642 
1643  // Set the output parameter
1644  *unique_tuples = all_tuples_list;
1645 
1646  return MB_SUCCESS;
1647 }
1648 
1649 // Calculate integrated field values for groups of entities
1650 ErrorCode Coupler::get_group_integ_vals( std::vector< std::vector< EntityHandle > >& groups,
1651  std::vector< double >& integ_vals,
1652  const char* norm_tag,
1653  int /*num_integ_vals*/,
1654  Coupler::IntegType integ_type )
1655 {
1656  ErrorCode err;
1657 
1658  std::vector< std::vector< EntityHandle > >::iterator iter_i;
1659  std::vector< EntityHandle >::iterator iter_j;
1660  double grp_intrgr_val, intgr_val;
1661 
1662  // Get the tag handle for norm_tag
1663  Tag norm_hdl;
1664  err =
1665  mbImpl->tag_get_handle( norm_tag, 1, moab::MB_TYPE_DOUBLE, norm_hdl, moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT );ERRORR( "Failed to get norm_tag handle.", err );
1666 
1667  // Check size of integ_vals vector
1668  if( integ_vals.size() != groups.size() ) integ_vals.resize( groups.size() );
1669 
1670  // Loop over the groups(vectors) of entities
1671  unsigned int i;
1672  for( i = 0, iter_i = groups.begin(); iter_i != groups.end(); i++, ++iter_i )
1673  {
1674  grp_intrgr_val = 0;
1675 
1676  // Loop over the all the entities in the group, integrating
1677  // the field_fn over the entity in iter_j
1678  for( iter_j = ( *iter_i ).begin(); iter_j != ( *iter_i ).end(); ++iter_j )
1679  {
1680  EntityHandle ehandle = ( *iter_j );
1681 
1682  // Check that the entity in iter_j is of the same dimension as the
1683  // integ_type we are performing
1684  EntityType j_type;
1685  j_type = mbImpl->type_from_handle( ehandle );ERRORR( "Failed to get entity type.", err );
1686  // Skip any entities in the group that are not of the type being considered
1687  if( ( integ_type == VOLUME ) && ( j_type < MBTET || j_type >= MBENTITYSET ) ) continue;
1688 
1689  intgr_val = 0;
1690 
1691  // Retrieve the vertices from the element
1692  const EntityHandle* verts = NULL;
1693  int connectivity_size = 0;
1694 
1695  err = mbImpl->get_connectivity( ehandle, verts, connectivity_size, false );ERRORR( "Failed to get vertices from entity.", err );
1696 
1697  // Get the vertex coordinates and the field values at the vertices.
1698  double* coords = (double*)malloc( sizeof( double ) * ( 3 * connectivity_size ) );
1699  /* TODO: VSM: check if this works for lower dimensions also without problems */
1700  /* if (3 == geom_dim) */
1701  err = mbImpl->get_coords( verts, connectivity_size, coords );ERRORR( "Failed to get vertex coordinates.", err );
1702 
1703  /* allocate the field data array */
1704  double* vfield = (double*)malloc( sizeof( double ) * ( connectivity_size ) );
1705  err = mbImpl->tag_get_data( norm_hdl, verts, connectivity_size, vfield );
1706  if( MB_SUCCESS != err )
1707  {
1708  free( coords );
1709  }
1710  ERRORR( "Failed to get vertex coordinates.", err );
1711 
1712  // Get coordinates of all corner vertices (in normal order) and
1713  // put in array of CartVec.
1714  std::vector< CartVect > vertices( connectivity_size );
1715 
1716  // Put the vertices into a CartVect vector
1717  double* x = coords;
1718  for( int j = 0; j < connectivity_size; j++, x += 3 )
1719  {
1720  vertices[j] = CartVect( x );
1721  }
1722  free( coords );
1723 
1724  moab::Element::Map* elemMap;
1725  if( j_type == MBHEX )
1726  {
1727  if( connectivity_size == 8 )
1728  elemMap = new moab::Element::LinearHex( vertices );
1729  else
1730  elemMap = new moab::Element::QuadraticHex( vertices );
1731  }
1732  else if( j_type == MBTET )
1733  {
1734  elemMap = new moab::Element::LinearTet( vertices );
1735  }
1736  else if( j_type == MBQUAD )
1737  {
1738  elemMap = new moab::Element::LinearQuad( vertices );
1739  }
1740  /*
1741  else if (j_type == MBTRI) {
1742  elemMap = new moab::Element::LinearTri(vertices);
1743  }
1744  */
1745  else if( j_type == MBEDGE )
1746  {
1747  elemMap = new moab::Element::LinearEdge( vertices );
1748  }
1749  else
1750  ERRORR( "Unknown topology type.", MB_UNSUPPORTED_OPERATION );
1751 
1752  // Set the vertices in the Map and perform the integration
1753  try
1754  {
1755  /* VSM: Do we need this call ?? */
1756  // elemMap->set_vertices(vertices);
1757 
1758  // Perform the actual integration over the element
1759  intgr_val = elemMap->integrate_scalar_field( vfield );
1760 
1761  // Combine the result with those of the group
1762  grp_intrgr_val += intgr_val;
1763  }
1765  {
1766  std::cerr << "Failed to set vertices on Element::Map." << std::endl;
1767  }
1769  {
1770  std::cerr << "Failed to get inverse evaluation of coordinate on Element::Map." << std::endl;
1771  }
1772 
1773  delete( elemMap );
1774  free( vfield );
1775  }
1776 
1777  // Set the group integrated value in the vector
1778  integ_vals[i] = grp_intrgr_val;
1779  }
1780 
1781  return err;
1782 }
1783 
1784 // Apply a normalization factor to group of entities
1785 ErrorCode Coupler::apply_group_norm_factor( std::vector< std::vector< EntityHandle > >& entity_sets,
1786  std::vector< double >& norm_factors,
1787  const char* norm_tag,
1788  Coupler::IntegType /*integ_type*/ )
1789 {
1790  ErrorCode err;
1791 
1792  // Construct the new tag for the normalization factor from the norm_tag name
1793  // and "_normf".
1794  int norm_tag_len = strlen( norm_tag );
1795  const char* normf_appd = "_normf";
1796  int normf_appd_len = strlen( normf_appd );
1797 
1798  char* normf_tag = (char*)malloc( norm_tag_len + normf_appd_len + 1 );
1799  char* tmp_ptr = normf_tag;
1800 
1801  memcpy( tmp_ptr, norm_tag, norm_tag_len );
1802  tmp_ptr += norm_tag_len;
1803  memcpy( tmp_ptr, normf_appd, normf_appd_len );
1804  tmp_ptr += normf_appd_len;
1805  *tmp_ptr = '\0';
1806 
1807  Tag normf_hdl;
1808  // Check to see if the tag exists. If not then create it and get the handle.
1809  err = mbImpl->tag_get_handle( normf_tag, 1, moab::MB_TYPE_DOUBLE, normf_hdl,
1810  moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT );ERRORR( "Failed to create normalization factor tag.", err );
1811  if( normf_hdl == NULL )
1812  {
1813  std::string msg( "Failed to create normalization factor tag named '" );
1814  msg += std::string( normf_tag ) + std::string( "'" );ERRORR( msg.c_str(), MB_FAILURE );
1815  }
1816  free( normf_tag );
1817 
1818  std::vector< std::vector< EntityHandle > >::iterator iter_i;
1819  std::vector< EntityHandle >::iterator iter_j;
1820  std::vector< double >::iterator iter_f;
1821  double grp_norm_factor = 0.0;
1822 
1823  // Loop over the entity sets
1824  for( iter_i = entity_sets.begin(), iter_f = norm_factors.begin();
1825  ( iter_i != entity_sets.end() ) && ( iter_f != norm_factors.end() ); ++iter_i, ++iter_f )
1826  {
1827  grp_norm_factor = *iter_f;
1828 
1829  // Loop over the all the entity sets in iter_i and set the
1830  // new normf_tag with the norm factor value on each
1831  for( iter_j = ( *iter_i ).begin(); iter_j != ( *iter_i ).end(); ++iter_j )
1832  {
1833  EntityHandle entset = *iter_j;
1834 
1835  std::cout << "Coupler: applying normalization for entity set=" << entset
1836  << ", normalization_factor=" << grp_norm_factor << std::endl;
1837 
1838  err = mbImpl->tag_set_data( normf_hdl, &entset, 1, &grp_norm_factor );ERRORR( "Failed to set normalization factor on entity set.", err );
1839  }
1840  }
1841 
1842  return MB_SUCCESS;
1843 }
1844 
1845 #define UINT_PER_X( X ) ( ( sizeof( X ) + sizeof( uint ) - 1 ) / sizeof( uint ) )
1846 #define UINT_PER_REAL UINT_PER_X( realType )
1847 #define UINT_PER_LONG UINT_PER_X( slong )
1848 #define UINT_PER_UNSIGNED UINT_PER_X( unsigned )
1849 
1850 // Function for packing tuple_list. Returns number of uints copied into buffer.
1851 int pack_tuples( TupleList* tl, void** ptr )
1852 {
1853  uint mi, ml, mul, mr;
1854  tl->getTupleSize( mi, ml, mul, mr );
1855 
1856  uint n = tl->get_n();
1857 
1858  int sz_buf = 1 + 4 * UINT_PER_UNSIGNED +
1859  tl->get_n() * ( mi + ml * UINT_PER_LONG + mul * UINT_PER_LONG + mr * UINT_PER_REAL );
1860 
1861  uint* buf = (uint*)malloc( sz_buf * sizeof( uint ) );
1862  *ptr = (void*)buf;
1863 
1864  // Copy n
1865  memcpy( buf, &n, sizeof( uint ) ), buf += 1;
1866  // Copy mi
1867  memcpy( buf, &mi, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
1868  // Copy ml
1869  memcpy( buf, &ml, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
1870  // Copy mul
1871  memcpy( buf, &mul, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
1872  // Copy mr
1873  memcpy( buf, &mr, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
1874  // Copy vi_wr
1875  memcpy( buf, tl->vi_rd, tl->get_n() * mi * sizeof( sint ) ), buf += tl->get_n() * mi;
1876  // Copy vl_wr
1877  memcpy( buf, tl->vl_rd, tl->get_n() * ml * sizeof( slong ) ), buf += tl->get_n() * ml * UINT_PER_LONG;
1878  // Copy vul_wr
1879  memcpy( buf, tl->vul_rd, tl->get_n() * mul * sizeof( ulong ) ), buf += tl->get_n() * mul * UINT_PER_LONG;
1880  // Copy vr_wr
1881  memcpy( buf, tl->vr_rd, tl->get_n() * mr * sizeof( realType ) ), buf += tl->get_n() * mr * UINT_PER_REAL;
1882 
1883  return sz_buf;
1884 }
1885 
1886 // Function for packing tuple_list
1887 void unpack_tuples( void* ptr, TupleList** tlp )
1888 {
1889  TupleList* tl = new TupleList();
1890  *tlp = tl;
1891 
1892  uint nt;
1893  unsigned mit, mlt, mult, mrt;
1894  uint* buf = (uint*)ptr;
1895 
1896  // Get n
1897  memcpy( &nt, buf, sizeof( uint ) ), buf += 1;
1898  // Get mi
1899  memcpy( &mit, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
1900  // Get ml
1901  memcpy( &mlt, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
1902  // Get mul
1903  memcpy( &mult, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
1904  // Get mr
1905  memcpy( &mrt, buf, sizeof( unsigned ) ), buf += UINT_PER_UNSIGNED;
1906 
1907  // Initialize tl
1908  tl->initialize( mit, mlt, mult, mrt, nt );
1909  tl->enableWriteAccess();
1910  tl->set_n( nt );
1911 
1912  uint mi, ml, mul, mr;
1913  tl->getTupleSize( mi, ml, mul, mr );
1914 
1915  // Get vi_rd
1916  memcpy( tl->vi_wr, buf, tl->get_n() * mi * sizeof( sint ) ), buf += tl->get_n() * mi;
1917  // Get vl_rd
1918  memcpy( tl->vl_wr, buf, tl->get_n() * ml * sizeof( slong ) ), buf += tl->get_n() * ml * UINT_PER_LONG;
1919  // Get vul_rd
1920  memcpy( tl->vul_wr, buf, tl->get_n() * mul * sizeof( ulong ) ), buf += tl->get_n() * mul * UINT_PER_LONG;
1921  // Get vr_rd
1922  memcpy( tl->vr_wr, buf, tl->get_n() * mr * sizeof( realType ) ), buf += tl->get_n() * mr * UINT_PER_REAL;
1923 
1924  tl->disableWriteAccess();
1925  return;
1926 }
1927 
1928 ErrorCode Coupler::get_gl_points_on_elements( Range& targ_elems, std::vector< double >& vpos, int& numPointsOfInterest )
1929 {
1930  numPointsOfInterest = targ_elems.size() * _ntot;
1931  vpos.resize( 3 * numPointsOfInterest );
1932  int ielem = 0;
1933  for( Range::iterator eit = targ_elems.begin(); eit != targ_elems.end(); ++eit, ielem += _ntot * 3 )
1934  {
1935  EntityHandle eh = *eit;
1936  const double* xval;
1937  const double* yval;
1938  const double* zval;
1939  ErrorCode rval = mbImpl->tag_get_by_ptr( _xm1Tag, &eh, 1, (const void**)&xval );
1940  if( moab::MB_SUCCESS != rval )
1941  {
1942  std::cout << "Can't get xm1 values \n";
1943  return MB_FAILURE;
1944  }
1945  rval = mbImpl->tag_get_by_ptr( _ym1Tag, &eh, 1, (const void**)&yval );
1946  if( moab::MB_SUCCESS != rval )
1947  {
1948  std::cout << "Can't get ym1 values \n";
1949  return MB_FAILURE;
1950  }
1951  rval = mbImpl->tag_get_by_ptr( _zm1Tag, &eh, 1, (const void**)&zval );
1952  if( moab::MB_SUCCESS != rval )
1953  {
1954  std::cout << "Can't get zm1 values \n";
1955  return MB_FAILURE;
1956  }
1957  // Now, in a stride, populate vpos
1958  for( int i = 0; i < _ntot; i++ )
1959  {
1960  vpos[ielem + 3 * i] = xval[i];
1961  vpos[ielem + 3 * i + 1] = yval[i];
1962  vpos[ielem + 3 * i + 2] = zval[i];
1963  }
1964  }
1965 
1966  return MB_SUCCESS;
1967 }
1968 
1969 } // namespace moab