Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
DataCoupler.cpp
Go to the documentation of this file.
1 #include "DataCoupler.hpp"
2 #include "moab/Tree.hpp"
3 #include "moab/TupleList.hpp"
5 #include "moab/ElemEvaluator.hpp"
6 #include "moab/Error.hpp"
7 
8 #ifdef MOAB_HAVE_MPI
9 #include "moab/ParallelComm.hpp"
10 #endif
11 
12 #include "iostream"
13 #include <algorithm>
14 #include <sstream>
15 
16 #include <cstdio>
17 #include <cassert>
18 
19 namespace moab
20 {
21 
23  Range& source_ents,
24  int coupler_id,
25  ParallelComm* pc,
26  bool init_locator,
27  int dim )
28  : mbImpl( impl ), myPcomm( pc ), myId( coupler_id ), myDim( dim )
29 {
30  assert( NULL != mbImpl && ( myPcomm || !source_ents.empty() ) );
31 
32  // Now initialize the tree
33  if( init_locator )
34  {
35  myLocator = new SpatialLocator( mbImpl, source_ents );
37 
38  // Initialize element evaluator with the default for the entity types in source_ents;
39  // can be replaced later by application if desired
40  if( !source_ents.empty() )
41  {
42  Range::pair_iterator pit = source_ents.pair_begin();
43  EntityType last_type = MBMAXTYPE;
44  for( ; pit != source_ents.pair_end(); ++pit )
45  {
46  EntityType this_type = mbImpl->type_from_handle( pit->first );
47  if( last_type == this_type ) continue;
48  ErrorCode rval = myLocator->elem_eval()->set_eval_set( pit->first );
49  if( MB_SUCCESS != rval ) throw( rval );
50  last_type = this_type;
51  }
52  }
53  }
54 
55  if( -1 == dim && !source_ents.empty() ) myDim = mbImpl->dimension_from_handle( *source_ents.rbegin() );
56 }
57 
58 /* Destructor
59  */
61 {
62  delete myLocator->elem_eval();
63  delete myLocator;
64 }
65 
67  const double rel_iter_tol,
68  const double abs_iter_tol,
69  const double inside_tol )
70 {
71  targetEnts = targ_ents;
72 
73 #ifdef MOAB_HAVE_MPI
74  if( myPcomm && myPcomm->size() > 1 )
75  return myLocator->par_locate_points( myPcomm, targ_ents, rel_iter_tol, abs_iter_tol, inside_tol );
76 #endif
77 
78  return myLocator->locate_points( targ_ents, rel_iter_tol, abs_iter_tol, inside_tol );
79 }
80 
82  int num_points,
83  const double rel_iter_tol,
84  const double abs_iter_tol,
85  const double inside_tol )
86 {
87 #ifdef MOAB_HAVE_MPI
88  if( myPcomm && myPcomm->size() > 1 )
89  return myLocator->par_locate_points( myPcomm, xyz, num_points, rel_iter_tol, abs_iter_tol, inside_tol );
90 #endif
91 
92  return myLocator->locate_points( xyz, num_points, rel_iter_tol, abs_iter_tol, inside_tol );
93 }
94 
95 ErrorCode DataCoupler::interpolate( /*DataCoupler::Method*/ int method,
96  const std::string& interp_tag,
97  double* interp_vals,
98  std::vector< int >* point_indices,
99  bool normalize )
100 {
101  // Tag name input, translate to tag handle and pass down the chain
102 
103  Tag tag;
104  ErrorCode result = mbImpl->tag_get_handle( interp_tag.c_str(), tag );MB_CHK_SET_ERR( result, "Failed to get handle for interpolation tag \"" << interp_tag << "\"" );
105 
106  return interpolate( method, tag, interp_vals, point_indices, normalize );
107 }
108 
109 ErrorCode DataCoupler::interpolate( /*DataCoupler::Method*/ int* /* methods */,
110  Tag* tags,
111  int* points_per_method,
112  int num_methods,
113  double* interp_vals,
114  std::vector< int >* point_indices,
115  bool /*normalize*/ )
116 {
117  // Lowest-level interpolate function, does actual interpolation using calls to ElemEvaluator
118 
119  ErrorCode result = MB_SUCCESS;
120 
121  unsigned int pts_total = 0;
122  for( int i = 0; i < num_methods; i++ )
123  pts_total += ( points_per_method ? points_per_method[i] : targetEnts.size() );
124 
125  unsigned int num_indices = ( point_indices ? point_indices->size() : targetEnts.size() );
126  // # points and indices should be identical
127  if( pts_total != num_indices ) return MB_FAILURE;
128 
129  // Since each tuple contains one interpolated tag, if we're interpolating multiple tags, every
130  // tuple needs to be able to store up to max tag size
131  int max_tsize = -1;
132  for( int i = 0; i < num_methods; i++ )
133  {
134  int tmp_tsize;
135  result = mbImpl->tag_get_length( tags[i], tmp_tsize );
136  if( MB_SUCCESS != result ) return MB_FAILURE;
137  max_tsize = std::max( max_tsize, tmp_tsize );
138  }
139 
140  if( myPcomm && myPcomm->size() > 1 )
141  {
142  // Build up the tuple list to distribute from my target points; assume that
143  // all procs use same method/tag input
144  TupleList TLob; // TLob structure: (pto_i, ridx_i, lidx_i, meth_tagidx_i)
145  TLob.initialize( 4, 0, 0, 0, num_indices );
146  int tn = 0; // The tuple number we're currently on
147  TLob.enableWriteAccess();
148  for( int m = 0; m < num_methods; m++ )
149  {
150  int num_points = ( points_per_method ? points_per_method[m] : targetEnts.size() );
151  for( int j = 0; j < num_points; j++ )
152  {
153  int idx = ( point_indices ? ( *point_indices )[j]
154  : j ); // The index in my targetEnts for this interpolation point
155 
156  // Remote proc/idx from myLocator->parLocTable
157  TLob.vi_wr[4 * tn] = myLocator->par_loc_table().vi_rd[2 * idx]; // proc
158  TLob.vi_wr[4 * tn + 1] = myLocator->par_loc_table().vi_rd[2 * idx + 1]; // Remote idx
159 
160  // Local entity index, tag/method index from my data
161  TLob.vi_wr[4 * tn + 2] = idx;
162  TLob.vi_wr[4 * tn + 3] = m;
163  TLob.inc_n();
164  tn++;
165  }
166  }
167 
168  // Scatter/gather interpolation points
169  myPcomm->proc_config().crystal_router()->gs_transfer( 1, TLob, 0 );
170 
171  // Perform interpolation on local source mesh; put results into TLinterp
172  TupleList TLinterp; // TLinterp structure: (pto_i, ridx_i, vals[max_tsize]_d)
173  TLinterp.initialize( 2, 0, 0, max_tsize, TLob.get_n() );
174  TLinterp.set_n( TLob.get_n() ); // Set the size right away
175  TLinterp.enableWriteAccess();
176 
177  for( unsigned int i = 0; i < TLob.get_n(); i++ )
178  {
179  int lidx = TLob.vi_rd[4 * i + 1]; // Index into myLocator's local table
180  // Method and tag indexed with same index
181  ///*Method*/ int method = methods[TLob.vi_rd[4*i + 3]];
182  Tag tag = tags[TLob.vi_rd[4 * i + 3]];
183  // Copy proc/remote index from TLob
184  TLinterp.vi_wr[2 * i] = TLob.vi_rd[4 * i];
185  TLinterp.vi_wr[2 * i + 1] = TLob.vi_rd[4 * i + 2];
186 
187  // Set up the evaluator for the tag and entity, then interpolate, putting result in
188  // TLinterp
191  result =
192  myLocator->elem_eval()->eval( myLocator->loc_table().vr_rd + 3 * lidx, TLinterp.vr_rd + i * max_tsize );
193  if( MB_SUCCESS != result ) return result;
194  }
195 
196  // Scatter/gather interpolation data
197  myPcomm->proc_config().crystal_router()->gs_transfer( 1, TLinterp, 0 );
198 
199  // Copy the interpolated field as a unit
200  std::copy( TLinterp.vr_rd, TLinterp.vr_rd + TLinterp.get_n() * max_tsize, interp_vals );
201  }
202  else
203  {
204  // If called serially, interpolate directly from local locations/entities,
205  // into either interp_vals or by setting data directly on tags on those entities
206  std::vector< double > tmp_vals;
207  std::vector< EntityHandle > tmp_ents;
208  double* tmp_dbl = interp_vals;
209  for( int i = 0; i < num_methods; i++ )
210  {
211  int num_points = ( points_per_method ? points_per_method[i] : targetEnts.size() );
212 
213  // Interpolated data is tsize long, which is either max size (if data passed back to
214  // caller in interp_vals) or tag size (if data will be set on entities, in which case it
215  // shouldn't have padding) set sizes here, inside loop over methods, to reduce memory
216  // usage
217  int tsize = max_tsize, tsize_bytes = 0;
218  if( !interp_vals )
219  {
220  tmp_vals.resize( num_points * max_tsize );
221  tmp_dbl = &tmp_vals[0];
222  tmp_ents.resize( num_points );
223  result = mbImpl->tag_get_length( tags[i], tsize );
224  if( MB_SUCCESS != result ) return result;
225  result = mbImpl->tag_get_bytes( tags[i], tsize_bytes );
226  if( MB_SUCCESS != result ) return result;
227  }
228 
229  for( int j = 0; j < num_points; j++ )
230  {
231  int lidx;
232  if( point_indices )
233  lidx = ( *point_indices )[j];
234  else
235  lidx = j;
236 
237  myLocator->elem_eval()->set_tag_handle( tags[i] );
239  if( !interp_vals ) tmp_ents[j] = targetEnts[lidx]; // Could be performance-sensitive, thus the if test
240  result = myLocator->elem_eval()->eval( myLocator->loc_table().vr_rd + 3 * lidx, tmp_dbl );
241  tmp_dbl += tsize;
242  if( MB_SUCCESS != result ) return result;
243  } // for j over points
244 
245  if( !interp_vals )
246  {
247  // Set tags on tmp_ents; data is already w/o padding, due to tsize setting above
248  result = mbImpl->tag_set_data( tags[i], &tmp_ents[0], tmp_ents.size(), &tmp_vals[0] );
249  if( MB_SUCCESS != result ) return result;
250  }
251 
252  } // for i over methods
253  } // if myPcomm
254 
255  // Done
256  return MB_SUCCESS;
257 }
258 
259 } // namespace moab