 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include "Coupler.hpp"
2 #include "moab/ParallelComm.hpp"
3 #include "moab/AdaptiveKDTree.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
47 Coupler::Coupler( Interface* impl,
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;
68 _spectralSource = _spectralTarget = NULL;
69 }
70
71 /* Destructor
72 */
73 Coupler::~Coupler()
74 {
75 // This will clear the cache
76 delete(moab::Element::SpectralHex*)_spectralSource;
77 delete(moab::Element::SpectralHex*)_spectralTarget;
78 delete myTree;
79 delete targetPts;
80 delete mappedPts;
81 }
82
83 ErrorCode Coupler::initialize_tree()
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
196 ErrorCode Coupler::initialize_spectral_elements( EntityHandle rootSource,
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];
256 rval = mbImpl->tag_get_handle( "SEM_X", _ntot, moab::MB_TYPE_DOUBLE, _xm1Tag );
257 if( moab::MB_SUCCESS != rval )
258 {
259 std::cout << "Can't get xm1tag \n";
260 return MB_FAILURE;
261 }
262 rval = mbImpl->tag_get_handle( "SEM_Y", _ntot, moab::MB_TYPE_DOUBLE, _ym1Tag );
263 if( moab::MB_SUCCESS != rval )
264 {
265 std::cout << "Can't get ym1tag \n";
266 return MB_FAILURE;
267 }
268 rval = mbImpl->tag_get_handle( "SEM_Z", _ntot, moab::MB_TYPE_DOUBLE, _zm1Tag );
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
319 ErrorCode Coupler::locate_points( double* xyz,
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() );
339 mappedPts->enableWriteAccess();
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
562 ErrorCode Coupler::test_local_box( double* xyz,
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;
588 myTree->get_bounding_box( box, &localRoot );
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
616 mappedPts->enableWriteAccess();
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
643 ErrorCode Coupler::interpolate( Coupler::Method method,
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
663 ErrorCode Coupler::interpolate( Coupler::Method* methods,
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
714 mappedPts->enableWriteAccess();
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
747 ErrorCode Coupler::nat_param( double xyz[3],
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 }
839 Element::SpectralHex* spcHex = (Element::SpectralHex*)_spectralSource;
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 }
855 catch( Element::Map::EvaluationError& )
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 }
890 catch( Element::Map::EvaluationError& )
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 }
905 catch( Element::Map::EvaluationError& )
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 }
933 catch( Element::Map::EvaluationError& )
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 }
950 catch( Element::Map::EvaluationError& )
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 }
966 catch( Element::Map::EvaluationError& )
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 }
996 catch( Element::Map::EvaluationError )
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 }
1031 Element::SpectralHex* spcHex = (Element::SpectralHex*)_spectralSource;
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 }
1095 catch( moab::Element::Map::EvaluationError& )
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.
1109 ErrorCode Coupler::constant_interp( EntityHandle elem, Tag tag, double& field )
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.
1123 ErrorCode Coupler::normalize_mesh( EntityHandle 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
1154 ErrorCode Coupler::normalize_subset( EntityHandle root_set,
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
1177 ErrorCode Coupler::normalize_subset( EntityHandle root_set,
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
1291 ErrorCode Coupler::get_matching_entities( EntityHandle root_set,
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
1313 ErrorCode Coupler::get_matching_entities( EntityHandle root_set,
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,
1382 MASTER_PROC, myPc->proc_config().proc_comm() );
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.
1504 ErrorCode Coupler::create_tuples( Range& ent_sets,
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 }
1764 catch( moab::Element::Map::ArgError& )
1765 {
1766 std::cerr << "Failed to set vertices on Element::Map." << std::endl;
1767 }
1768 catch( moab::Element::Map::EvaluationError& )
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
1.9.1.