1 /**
2 * \class moab::Coupler
3 * \author Tim Tautges
4 *
5 * \brief This class couples data between meshes.
6 *
7 * The coupler interpolates solution data at a set of points. Data
8 * being interpolated resides on a source mesh, in a tag.
9 * Applications calling this coupler send in entities, usually points
10 * or vertices, and receive back the tag value interpolated at those
11 * points. Entities in the source mesh containing those points
12 * do not have to reside on the same processor.
13 *
14 * To use, an application should:
15 * - instantiate this coupler by calling the constructor collectively
16 * on all processors in the communicator
17 * - call locate_points, which locates the points to be interpolated and
18 * (optionally) caches the results in this object
19 * - call interpolate, which does the interpolation
20 *
21 * Multiple interpolations can be done after locating the points.
22 *
23 */
24 #ifndef COUPLER_HPP
25 #define COUPLER_HPP
26
27 #include "moab/Range.hpp"
28 #include "moab/Interface.hpp"
29 #include "moab/CartVect.hpp"
30 #include "moab/TupleList.hpp"
31
32 #include <sstream>
33
34 namespace moab {
35
36 class ParallelComm;
37
38 class AdaptiveKDTree;
39
40 class TupleList;
41
42 class Coupler
43 {
44 public:
45
46 enum Method {CONSTANT, LINEAR_FE, QUADRATIC_FE, SPECTRAL, SPHERICAL};
47
48 enum IntegType {VOLUME};
49
50 /* Constructor
51 * Constructor, which also optionally initializes the coupler
52 * \param pc ParallelComm object to be used with this coupler, representing the union
53 * of processors containing source and target meshes
54 * \param local_elems Local elements in the source mesh
55 * \param coupler_id Id of this coupler, should be the same over all procs
56 * \param init_tree If true, initializes kdtree inside the constructor
57 */
58 Coupler(Interface *impl,
59 ParallelComm *pc,
60 Range &local_elems,
61 int coupler_id,
62 bool init_tree = true,
63 int max_ent_dim = 3);
64
65 /* Destructor
66 */
67 virtual ~Coupler();
68
69 /**
70 * \brief Initialize the kdtree, locally and across communicator
71 */
72 ErrorCode initialize_tree();
73
74 /* \brief Locate points on the source mesh
75 * This function finds the element/processor/natural coordinates for the
76 * source mesh element containing each point, optionally storing the results
77 * on the target mesh processor. Relative tolerance is compared to bounding
78 * box diagonal length. Tolerance is compared to [-1,1] parametric extent
79 * in the reference element.
80 * \param xyz Point locations (interleaved) being located
81 * \param num_points Number of points in xyz
82 * \param rel_eps Relative tolerance for the non-linear iteration inside a given element
83 * \param abs_eps Absolute tolerance for the non-linear iteration inside a given element
84 * \param tl Tuple list containing the results, with each tuple
85 * consisting of (p, i), p = proc, i = index on that proc
86 * \param store_local If true, stores the tuple list in targetPts
87 */
88 ErrorCode locate_points(double *xyz, unsigned int num_points,
89 double rel_eps = 0.0,
90 double abs_eps = 0.0,
91 TupleList *tl = NULL,
92 bool store_local = true);
93
94 /* \brief Locate entities on the source mesh
95 * This function finds the element/processor/natural coordinates for the
96 * source mesh element containing each entity, optionally storing the results
97 * on the target mesh processor. Location of each target mesh entity passed in
98 * is its centroid (for non-vertices, the avg of its vertex positions).
99 * Relative tolerance is compared to bounding
100 * box diagonal length. Tolerance is compared to [-1,1] parametric extent
101 * in the reference element.
102 * \param ents Entities being located
103 * \param rel_eps Relative tolerance for the non-linear iteration inside a given element
104 * \param abs_eps Absolute tolerance for the non-linear iteration inside a given element
105 * \param tl Tuple list containing the results, with each tuple
106 * consisting of (p, i), p = proc, i = index on that proc
107 * \param store_local If true, stores the tuple list in targetPts
108 *
109 */
110 ErrorCode locate_points(Range &ents,
111 double rel_eps = 0.0,
112 double abs_eps = 0.0,
113 TupleList *tl = NULL,
114 bool store_local = true);
115
116 /* \brief Interpolate data from the source mesh onto points
117 * All entities/points or, if tuple_list is input, only those points
118 * are interpolated from the source mesh. Application should
119 * allocate enough memory in interp_vals to hold interpolation results.
120 *
121 * If normalization is requested, technique used depends on the coupling
122 * method.
123 *
124 * \param method Interpolation/normalization method
125 * \param tag Tag on source mesh holding data to be interpolated
126 * \param interp_vals Memory holding interpolated data
127 * \param tl Tuple list of points to be interpolated, in format used by targetPts
128 * (see documentation for targetPts below); if NULL, all locations
129 * in targetPts are interpolated
130 * \param normalize If true, normalization is done according to method
131 */
132 ErrorCode interpolate(Coupler::Method method,
133 Tag tag,
134 double *interp_vals,
135 TupleList *tl = NULL,
136 bool normalize = true);
137
138 /* \brief Interpolate data from the source mesh onto points
139 * All entities/points or, if tuple_list is input, only those points
140 * are interpolated from the source mesh. Application should
141 * allocate enough memory in interp_vals to hold interpolation results.
142 *
143 * If normalization is requested, technique used depends on the coupling
144 * method.
145 *
146 * \param methods Interpolation/normalization method
147 * \param tag_name Name of tag on source mesh holding data to be interpolated
148 * \param interp_vals Memory holding interpolated data
149 * \param tl Tuple list of points to be interpolated, in format used by targetPts
150 * (see documentation for targetPts below); if NULL, all locations
151 * in targetPts are interpolated
152 * \param normalize If true, normalization is done according to method
153 */
154 ErrorCode interpolate(Coupler::Method method,
155 const std::string &tag_name,
156 double *interp_vals,
157 TupleList *tl = NULL,
158 bool normalize = true);
159
160 /* \brief Interpolate data from multiple tags
161 * All entities/points or, if tuple_list is input, only those points
162 * are interpolated from the source mesh. Application should
163 * allocate enough memory in interp_vals to hold interpolation results.
164 *
165 * In this variant, multiple tags, possibly with multiple interpolation
166 * methods, are specified. Sum of values in points_per_method should be
167 * the number of points in tl or, if NULL, targetPts.
168 *
169 * If normalization is requested, technique used depends on the coupling
170 * method.
171 *
172 * \param methods Vector of Interpolation/normalization methods
173 * \param tag_names Names of tag being interpolated for each method
174 * \param points_per_method Number of points for each method
175 * \param num_methods Length of vectors in previous 3 arguments
176 * \param interp_vals Memory holding interpolated data
177 * \param tl Tuple list of points to be interpolated; if NULL, all locations
178 * stored in this object are interpolated
179 * \param normalize If true, normalization is done according to method
180 */
181 ErrorCode interpolate(Coupler::Method *methods,
182 const std::string *tag_names,
183 int *points_per_method,
184 int num_methods,
185 double *interp_vals,
186 TupleList *tl = NULL,
187 bool normalize = true);
188
189 /* \brief Interpolate data from multiple tags
190 * All entities/points or, if tuple_list is input, only those points
191 * are interpolated from the source mesh. Application should
192 * allocate enough memory in interp_vals to hold interpolation results.
193 *
194 * In this variant, multiple tags, possibly with multiple interpolation
195 * methods, are specified. Sum of values in points_per_method should be
196 * the number of points in tl or, if NULL, targetPts.
197 *
198 * If normalization is requested, technique used depends on the coupling
199 * method.
200 *
201 * \param methods Vector of Interpolation/normalization methods
202 * \param tag_names Names of tag being interpolated for each method
203 * \param points_per_method Number of points for each method
204 * \param num_methods Length of vectors in previous 3 arguments
205 * \param interp_vals Memory holding interpolated data
206 * \param tl Tuple list of points to be interpolated; if NULL, all locations
207 * stored in this object are interpolated
208 * \param normalize If true, normalization is done according to method
209 */
210 ErrorCode interpolate(Coupler::Method *methods,
211 Tag *tag_names,
212 int *points_per_method,
213 int num_methods,
214 double *interp_vals,
215 TupleList *tl = NULL,
216 bool normalize = true);
217
218 /* \brief Normalize a field over an entire mesh
219 * A field existing on the vertices of elements of a mesh is integrated
220 * over all elements in the mesh. The integrated value is normalized
221 * and the normalization factor is saved to a new tag
222 * on the mesh entity set.
223 *
224 * \param root_set Entity Set representing the entire mesh
225 * \param norm_tag Tag containing field data to integrate
226 * \param integ_type Type of integration to perform
227 * \param num_integ_pts The number of Gaussian integration points to use in each dimension
228 */
229 ErrorCode normalize_mesh( EntityHandle root_set,
230 const char *norm_tag,
231 Coupler::IntegType integ_type,
232 int num_integ_pts);
233
234 /* \brief Normalize a field over subsets of entities
235 * A field existing on the vertices of elements of a mesh is integrated
236 * over subsets of elements identified by the tags and values. The integrated
237 * values are normalized and the normalization factor is saved to a new tag
238 * on the entity sets which contain the elements of a subset.
239 *
240 * \param root_set Entity Set from the mesh from which to select subsets
241 * \param norm_tag Tag containing field data to integrate
242 * \param tag_names Array of tag names used for selecting element subsets
243 * \param num_tags Number of tag names
244 * \param tag_values Array of tag values passed as strings; the array will be
245 * the same length as that for tag names however some entries may be
246 * NULL indicating that tag should be matched for existence and not value
247 * \param integ_type Type of integration to perform
248 * \param num_integ_pts The number of Gaussian integration points to use in each dimension
249 */
250 ErrorCode normalize_subset( EntityHandle root_set,
251 const char *norm_tag,
252 const char **tag_names,
253 int num_tags,
254 const char **tag_values,
255 Coupler::IntegType integ_type,
256 int num_integ_pts);
257
258 /* \brief Normalize a field over subsets of entities
259 * A field existing on the vertices of elements of a mesh is integrated
260 * over subsets of elements identified by the tags and values. The integrated
261 * values are normalized and the normalization factor is saved to a new tag
262 * on the entity sets which contain the elements of a subset.
263 *
264 * \param root_set Entity Set from the mesh from which to select subsets
265 * \param norm_tag Tag containing field data to integrate
266 * \param tag_handles Array of tag handles used for selecting element subsets
267 * \param num_tags Number of tag handles
268 * \param tag_values Array of tag values passed as strings; the array will be
269 * the same length as that for tag handles however some entries may be
270 * NULL indicating that tag should be matched for existence and not value
271 * \param integ_type Type of integration to perform
272 * \param num_integ_pts The number of Gaussian integration points to use in each dimension
273 */
274 ErrorCode normalize_subset( EntityHandle root_set,
275 const char *norm_tag,
276 Tag *tag_handles,
277 int num_tags,
278 const char **tag_values,
279 Coupler::IntegType integ_type,
280 int num_integ_pts);
281
282 /* \brief Retrieve groups of entities matching tags and values(if present)
283 * Retrieve a vector of vectors of entity handles matching the
284 * tags and values. The entity set passed is used as the search domain.
285 *
286 * \param norm_tag Tag containing field data to integrate
287 * \param entity_sets Pointer to vector of vectors of entity set handles
288 * \param entity_groups Pointer to vector of vectors of entity handles from each entity set
289 * \param integ_type Type of integration to perform
290 * \param num_integ_pts The number of Gaussian integration points to use in each dimension
291 */
292 ErrorCode do_normalization( const char *norm_tag,
293 std::vector< std::vector<EntityHandle> > &entity_sets,
294 std::vector< std::vector<EntityHandle> > &entity_groups,
295 Coupler::IntegType integ_type,
296 int num_integ_pts);
297
298 /* \brief Retrieve groups of entities matching tags and values(if present)
299 * Retrieve a vector of vectors of entity handles matching the
300 * tags and values. The entity set passed is used as the search domain.
301 *
302 * \param root_set Set from which to search for matching entities
303 * \param tag_names Array of tag names used to select entities
304 * \param tag_values Array of tag values used to select entities
305 * \param num_tags Number of tag names
306 * \param entity_sets Pointer to vector of vectors of entity set handles found in the search
307 * \param entity_groups Pointer to vector of vectors of entity handles from each entity set
308 */
309 ErrorCode get_matching_entities( EntityHandle root_set,
310 const char **tag_names,
311 const char **tag_values,
312 int num_tags,
313 std::vector< std::vector<EntityHandle> > *entity_sets,
314 std::vector< std::vector<EntityHandle> > *entity_groups);
315
316 /* \brief Retrieve groups of entities matching tags and values(if present)
317 * Retrieve a vector of vectors of entity handles matching the
318 * tags and values. The entity set passed is used as the search domain.
319 *
320 * \param root_set Set from which to search for matching entities
321 * \param tag_handles Array of tag handles used to select entities
322 * \param tag_values Array of tag values used to select entities
323 * \param num_tags Number of tag handles
324 * \param entity_sets Pointer to vector of vectors of entity set handles found in the search
325 * \param entity_groups Pointer to vector of vectors of entity handles from each entity set
326 */
327 ErrorCode get_matching_entities( EntityHandle root_set,
328 Tag *tag_handles,
329 const char **tag_values,
330 int num_tags,
331 std::vector< std::vector<EntityHandle> > *entity_sets,
332 std::vector< std::vector<EntityHandle> > *entity_groups);
333
334 /* \brief Return an array of tuples of tag values for each Entity Set
335 * A list of n-tuples will be constructed with 1 n-tuple for each Entity Set.
336 * The n-tuple will have an component for each tag given. It is assumed all
337 * of the tags are integer tags.
338 *
339 * \param ent_sets Array of Entity Set handles to use for retrieving tag data
340 * \param num_sets Number of Entity Sets
341 * \param tag_names Array of tag names
342 * \param num_tags Number of tag names
343 * \param tuples The returned tuple_list structure
344 */
345 ErrorCode create_tuples( Range &ent_sets,
346 const char **tag_names,
347 unsigned int num_tags,
348 TupleList **tuples);
349
350 /* \brief Return an array of tuples of tag values for each Entity Set
351 * A list of n-tuples will be constructed with 1 n-tuple for each Entity Set.
352 * The n-tuple will have an component for each tag given. It is assumed all
353 * of the tags are integer tags.
354 *
355 * \param ent_sets Array of Entity Set handles to use for retrieving tag data
356 * \param num_sets Number of Entity Sets
357 * \param tag_handles Array of tag handles
358 * \param num_tags Number of tag handles
359 * \param tuples The returned tuple_list structure
360 */
361 ErrorCode create_tuples( Range &ent_sets,
362 Tag *tag_handles,
363 unsigned int num_tags,
364 TupleList **tuples);
365
366 /* \brief Consolidate an array of n-tuples lists into one n-tuple list with no duplicates
367 * An array of list of n-tuples are consolidated into a single list of n-tuples
368 * with all duplicates removed. Only integer columns in the tuple_list are assumed to
369 * be used.
370 *
371 * \param all_tuples Array of tuple_lists to consolidate to one
372 * \param num_tuples Number of tuple_lists
373 * \param unique_tuples The consolidated tuple_list with no duplicates
374 */
375 ErrorCode consolidate_tuples( TupleList **all_tuples,
376 unsigned int num_tuples,
377 TupleList **unique_tuples);
378
379 /* \brief Calculate integrated field values for groups of entities
380 * An integrated field value, as defined by the field function,
381 * is calculated for each group of entities passed in.
382 *
383 * \param groups The vector contains vectors of entity handles, each representing a group
384 * \param integ_vals The integrated field values for each group
385 * \param norm_tag The tag name of the vertex-based field to be integrated
386 * \param num_integ_pts The number of Gaussian integration points to use in each dimension
387 * \param integ_type Type of integration to perform
388 */
389 ErrorCode get_group_integ_vals( std::vector<std::vector<EntityHandle> > &groups,
390 std::vector<double> &integ_vals,
391 const char *norm_tag,
392 int num_integ_pts,
393 Coupler::IntegType integ_type);
394
395 /* \brief Apply a normalization factor to group of entities
396 * Multiply a normalization factor with the value of norm_tag for each vertex
397 * of each entity in a group. Save the value back to norm_tag on each vertex.
398 *
399 * \param entity_sets The vector contains vectors of entity set handles, each containing the members of a group
400 * \param norm_factors The normalization factors for each group
401 * \param norm_tag The tag to be normalized on each group
402 * \param integ_type Type of integration to perform
403 */
404 ErrorCode apply_group_norm_factor( std::vector<std::vector<EntityHandle> > &entity_sets,
405 std::vector<double> &norm_factors,
406 const char *norm_tag,
407 Coupler::IntegType integ_type);
408
409 /*
410 * this method will look at source (and target sets?) sets, and look for the SEM_DIMS tag
411 * if it exists, it will trigger a spectral element caching, with the order specified
412 */
413 ErrorCode initialize_spectral_elements(EntityHandle rootSource, EntityHandle rootTarget,
414 bool &specSou, bool &specTar);
415
416 /*
417 * this method will put in an array, interleaved, the points of interest for coupling
418 * with a target mesh (so where do we need to compute the field of interest)
419 */
420 ErrorCode get_gl_points_on_elements(Range &targ_elems, std::vector<double> &vpos, int &numPointsOfInterest);
421
422 /* Get functions */
423
mb_impl() const424 inline Interface *mb_impl() const { return mbImpl; }
my_tree() const425 inline AdaptiveKDTree *my_tree() const { return myTree; }
local_root() const426 inline EntityHandle local_root() const { return localRoot; }
all_boxes() const427 inline const std::vector<double> &all_boxes() const { return allBoxes; }
my_pc() const428 inline ParallelComm *my_pc() const { return myPc; }
target_ents() const429 inline const Range &target_ents() const { return targetEnts; }
my_id() const430 inline int my_id() const { return myId; }
my_range() const431 inline const Range &my_range() const { return myRange; }
mapped_pts() const432 inline TupleList *mapped_pts() const { return mappedPts; }
num_its() const433 inline int num_its() const { return numIts; }
434 // used for spherical tests
set_spherical(bool arg1=true)435 inline void set_spherical (bool arg1=true) {spherical=arg1;}
436
437 private:
438
439 // Given a coordinate position, find all entities containing
440 // the point and the natural coords in those ents
441 ErrorCode nat_param(double xyz[3],
442 std::vector<EntityHandle> &entities,
443 std::vector<CartVect> &nat_coords,
444 double epsilon = 0.0);
445
446 ErrorCode interp_field(EntityHandle elem,
447 CartVect nat_coord,
448 Tag tag,
449 double &field);
450
451 ErrorCode constant_interp(EntityHandle elem,
452 Tag tag,
453 double &field);
454
455 ErrorCode test_local_box(double *xyz,
456 int from_proc, int remote_index, int index,
457 bool &point_located,
458 double rel_eps = 0.0,
459 double abs_eps = 0.0,
460 TupleList *tl = NULL);
461
462 /* \brief MOAB instance
463 */
464 Interface *mbImpl;
465
466 /* \brief Kdtree for local mesh
467 */
468 AdaptiveKDTree *myTree;
469
470 /* \brief Local root of the kdtree
471 */
472 EntityHandle localRoot;
473
474 /* \brief Min/max bounding boxes for all proc tree roots
475 */
476 std::vector<double> allBoxes;
477
478 /* \brief ParallelComm object for this coupler
479 */
480 ParallelComm *myPc;
481
482 /* \brief Id of this coupler
483 */
484 int myId;
485
486 /* \brief Range of source elements
487 */
488 Range myRange;
489
490 /* \brief Range of target entities
491 */
492 Range targetEnts;
493
494 /* \brief List of locally mapped tuples
495 * Tuples contain the following:
496 * n = # mapped points
497 * vul[i] = local handle of mapped entity
498 * vr[3*i..3*i+2] = natural coordinates in mapped entity
499 */
500 TupleList *mappedPts;
501
502 /* \brief Tuple list of target points and interpolated data
503 * Tuples contain the following:
504 * n = # target points
505 * vi[3*i] = remote proc mapping target point
506 * vi[3*i+1] = local index of target point
507 * vi[3*i+2] = remote index of target point
508 */
509 TupleList *targetPts;
510
511 /* \brief Number of iterations of tree building before failing
512 *
513 */
514 int numIts;
515
516 // Entity dimension
517 int max_dim;
518
519 // A cached spectral element for source and target, separate
520 // Assume that their number of GL points (order + 1) does not change
521 // If it does change, we need to reinitialize it
522 void * _spectralSource;
523 void * _spectralTarget;
524 moab::Tag _xm1Tag, _ym1Tag, _zm1Tag;
525 int _ntot;
526
527 // spherical coupling
528 bool spherical;
529 };
530
interpolate(Coupler::Method method,Tag tag,double * interp_vals,TupleList * tl,bool normalize)531 inline ErrorCode Coupler::interpolate(Coupler::Method method,
532 Tag tag,
533 double *interp_vals,
534 TupleList *tl,
535 bool normalize)
536 {
537 int num_pts = (tl ? tl->get_n() : targetPts->get_n());
538 return interpolate(&method, &tag, &num_pts, 1,
539 interp_vals, tl, normalize);
540 }
541
542 } // namespace moab
543
544 #endif
545