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