Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
Coupler.hpp
Go to the documentation of this file.
00001 /**
00002  * \class moab::Coupler
00003  * \author Tim Tautges
00004  *
00005  * \brief This class couples data between meshes.
00006  *
00007  * The coupler interpolates solution data at a set of points.  Data
00008  * being interpolated resides on a source mesh, in a tag.
00009  * Applications calling this coupler send in entities, usually points
00010  * or vertices, and receive back the tag value interpolated at those
00011  * points.  Entities in the source mesh containing those points
00012  * do not have to reside on the same processor.
00013  *
00014  * To use, an application should:
00015  * - instantiate this coupler by calling the constructor collectively
00016  *   on all processors in the communicator
00017  * - call locate_points, which locates the points to be interpolated and
00018  *   (optionally) caches the results in this object
00019  * - call interpolate, which does the interpolation
00020  *
00021  * Multiple interpolations can be done after locating the points.
00022  *
00023  */
00024 #ifndef COUPLER_HPP
00025 #define COUPLER_HPP
00026 
00027 #include "moab/Range.hpp"
00028 #include "moab/Interface.hpp"
00029 #include "moab/CartVect.hpp"
00030 #include "moab/TupleList.hpp"
00031 
00032 #include <sstream>
00033 
00034 namespace moab
00035 {
00036 
00037 class ParallelComm;
00038 
00039 class AdaptiveKDTree;
00040 
00041 class TupleList;
00042 
00043 class Coupler
00044 {
00045   public:
00046     enum Method
00047     {
00048         CONSTANT,
00049         LINEAR_FE,
00050         QUADRATIC_FE,
00051         SPECTRAL,
00052         SPHERICAL
00053     };
00054 
00055     enum IntegType
00056     {
00057         VOLUME
00058     };
00059 
00060     /* Constructor
00061      * Constructor, which also optionally initializes the coupler
00062      * \param pc ParallelComm object to be used with this coupler, representing the union
00063      *    of processors containing source and target meshes
00064      * \param local_elems Local elements in the source mesh
00065      * \param coupler_id Id of this coupler, should be the same over all procs
00066      * \param init_tree If true, initializes kdtree inside the constructor
00067      */
00068     Coupler( Interface* impl,
00069              ParallelComm* pc,
00070              Range& local_elems,
00071              int coupler_id,
00072              bool init_tree  = true,
00073              int max_ent_dim = 3 );
00074 
00075     /* Destructor
00076      */
00077     virtual ~Coupler();
00078 
00079     /**
00080      *  \brief Initialize the kdtree, locally and across communicator
00081      */
00082     ErrorCode initialize_tree();
00083 
00084     /* \brief Locate points on the source mesh
00085      * This function finds the element/processor/natural coordinates for the
00086      * source mesh element containing each point, optionally storing the results
00087      * on the target mesh processor.  Relative tolerance is compared to bounding
00088      * box diagonal length.  Tolerance is compared to [-1,1] parametric extent
00089      * in the reference element.
00090      * \param xyz Point locations (interleaved) being located
00091      * \param num_points Number of points in xyz
00092      * \param rel_eps Relative tolerance for the non-linear iteration inside a given element
00093      * \param abs_eps Absolute tolerance for the non-linear iteration inside a given element
00094      * \param tl Tuple list containing the results, with each tuple
00095      *           consisting of (p, i), p = proc, i = index on that proc
00096      * \param store_local If true, stores the tuple list in targetPts
00097      */
00098     ErrorCode locate_points( double* xyz,
00099                              unsigned int num_points,
00100                              double rel_eps   = 0.0,
00101                              double abs_eps   = 0.0,
00102                              TupleList* tl    = NULL,
00103                              bool store_local = true );
00104 
00105     /* \brief Locate entities on the source mesh
00106      * This function finds the element/processor/natural coordinates for the
00107      * source mesh element containing each entity, optionally storing the results
00108      * on the target mesh processor.  Location of each target mesh entity passed in
00109      * is its centroid (for non-vertices, the avg of its vertex positions).
00110      * Relative tolerance is compared to bounding
00111      * box diagonal length.  Tolerance is compared to [-1,1] parametric extent
00112      * in the reference element.
00113      * \param ents Entities being located
00114      * \param rel_eps Relative tolerance for the non-linear iteration inside a given element
00115      * \param abs_eps Absolute tolerance for the non-linear iteration inside a given element
00116      * \param tl Tuple list containing the results, with each tuple
00117      *           consisting of (p, i), p = proc, i = index on that proc
00118      * \param store_local If true, stores the tuple list in targetPts
00119      *
00120      */
00121     ErrorCode locate_points( Range& ents,
00122                              double rel_eps   = 0.0,
00123                              double abs_eps   = 0.0,
00124                              TupleList* tl    = NULL,
00125                              bool store_local = true );
00126 
00127     /* \brief Interpolate data from the source mesh onto points
00128      * All entities/points or, if tuple_list is input, only those points
00129      * are interpolated from the source mesh.  Application should
00130      * allocate enough memory in interp_vals to hold interpolation results.
00131      *
00132      * If normalization is requested, technique used depends on the coupling
00133      * method.
00134      *
00135      * \param method Interpolation/normalization method
00136      * \param tag Tag on source mesh holding data to be interpolated
00137      * \param interp_vals Memory holding interpolated data
00138      * \param tl Tuple list of points to be interpolated, in format used by targetPts
00139      *    (see documentation for targetPts below); if NULL, all locations
00140      *    in targetPts are interpolated
00141      * \param normalize If true, normalization is done according to method
00142      */
00143     ErrorCode interpolate( Coupler::Method method,
00144                            Tag tag,
00145                            double* interp_vals,
00146                            TupleList* tl  = NULL,
00147                            bool normalize = true );
00148 
00149     /* \brief Interpolate data from the source mesh onto points
00150      * All entities/points or, if tuple_list is input, only those points
00151      * are interpolated from the source mesh.  Application should
00152      * allocate enough memory in interp_vals to hold interpolation results.
00153      *
00154      * If normalization is requested, technique used depends on the coupling
00155      * method.
00156      *
00157      * \param methods Interpolation/normalization method
00158      * \param tag_name Name of tag on source mesh holding data to be interpolated
00159      * \param interp_vals Memory holding interpolated data
00160      * \param tl Tuple list of points to be interpolated, in format used by targetPts
00161      *    (see documentation for targetPts below); if NULL, all locations
00162      *    in targetPts are interpolated
00163      * \param normalize If true, normalization is done according to method
00164      */
00165     ErrorCode interpolate( Coupler::Method method,
00166                            const std::string& tag_name,
00167                            double* interp_vals,
00168                            TupleList* tl  = NULL,
00169                            bool normalize = true );
00170 
00171     /* \brief Interpolate data from multiple tags
00172      * All entities/points or, if tuple_list is input, only those points
00173      * are interpolated from the source mesh.  Application should
00174      * allocate enough memory in interp_vals to hold interpolation results.
00175      *
00176      * In this variant, multiple tags, possibly with multiple interpolation
00177      * methods, are specified.  Sum of values in points_per_method should be
00178      * the number of points in tl or, if NULL, targetPts.
00179      *
00180      * If normalization is requested, technique used depends on the coupling
00181      * method.
00182      *
00183      * \param methods Vector of Interpolation/normalization methods
00184      * \param tag_names Names of tag being interpolated for each method
00185      * \param points_per_method Number of points for each method
00186      * \param num_methods Length of vectors in previous 3 arguments
00187      * \param interp_vals Memory holding interpolated data
00188      * \param tl Tuple list of points to be interpolated; if NULL, all locations
00189      *       stored in this object are interpolated
00190      * \param normalize If true, normalization is done according to method
00191      */
00192     ErrorCode interpolate( Coupler::Method* methods,
00193                            const std::string* tag_names,
00194                            int* points_per_method,
00195                            int num_methods,
00196                            double* interp_vals,
00197                            TupleList* tl  = NULL,
00198                            bool normalize = true );
00199 
00200     /* \brief Interpolate data from multiple tags
00201      * All entities/points or, if tuple_list is input, only those points
00202      * are interpolated from the source mesh.  Application should
00203      * allocate enough memory in interp_vals to hold interpolation results.
00204      *
00205      * In this variant, multiple tags, possibly with multiple interpolation
00206      * methods, are specified.  Sum of values in points_per_method should be
00207      * the number of points in tl or, if NULL, targetPts.
00208      *
00209      * If normalization is requested, technique used depends on the coupling
00210      * method.
00211      *
00212      * \param methods Vector of Interpolation/normalization methods
00213      * \param tag_names Names of tag being interpolated for each method
00214      * \param points_per_method Number of points for each method
00215      * \param num_methods Length of vectors in previous 3 arguments
00216      * \param interp_vals Memory holding interpolated data
00217      * \param tl Tuple list of points to be interpolated; if NULL, all locations
00218      *       stored in this object are interpolated
00219      * \param normalize If true, normalization is done according to method
00220      */
00221     ErrorCode interpolate( Coupler::Method* methods,
00222                            Tag* tag_names,
00223                            int* points_per_method,
00224                            int num_methods,
00225                            double* interp_vals,
00226                            TupleList* tl  = NULL,
00227                            bool normalize = true );
00228 
00229     /* \brief Normalize a field over an entire mesh
00230      * A field existing on the vertices of elements of a mesh is integrated
00231      * over all elements in the mesh.  The integrated value is normalized
00232      * and the normalization factor is saved to a new tag
00233      * on the mesh entity set.
00234      *
00235      * \param root_set Entity Set representing the entire mesh
00236      * \param norm_tag Tag containing field data to integrate
00237      * \param integ_type Type of integration to perform
00238      * \param num_integ_pts The number of Gaussian integration points to use in each dimension
00239      */
00240     ErrorCode normalize_mesh( EntityHandle root_set,
00241                               const char* norm_tag,
00242                               Coupler::IntegType integ_type,
00243                               int num_integ_pts );
00244 
00245     /* \brief Normalize a field over subsets of entities
00246      * A field existing on the vertices of elements of a mesh is integrated
00247      * over subsets of elements identified by the tags and values.  The integrated
00248      * values are normalized and the normalization factor is saved to a new tag
00249      * on the entity sets which contain the elements of a subset.
00250      *
00251      * \param root_set Entity Set from the mesh from which to select subsets
00252      * \param norm_tag Tag containing field data to integrate
00253      * \param tag_names Array of tag names used for selecting element subsets
00254      * \param num_tags Number of tag names
00255      * \param tag_values Array of tag values passed as strings; the array will be
00256      *       the same length as that for tag names however some entries may be
00257      *       NULL indicating that tag should be matched for existence and not value
00258      * \param integ_type Type of integration to perform
00259      * \param num_integ_pts The number of Gaussian integration points to use in each dimension
00260      */
00261     ErrorCode normalize_subset( EntityHandle root_set,
00262                                 const char* norm_tag,
00263                                 const char** tag_names,
00264                                 int num_tags,
00265                                 const char** tag_values,
00266                                 Coupler::IntegType integ_type,
00267                                 int num_integ_pts );
00268 
00269     /* \brief Normalize a field over subsets of entities
00270      * A field existing on the vertices of elements of a mesh is integrated
00271      * over subsets of elements identified by the tags and values.  The integrated
00272      * values are normalized and the normalization factor is saved to a new tag
00273      * on the entity sets which contain the elements of a subset.
00274      *
00275      * \param root_set Entity Set from the mesh from which to select subsets
00276      * \param norm_tag Tag containing field data to integrate
00277      * \param tag_handles Array of tag handles used for selecting element subsets
00278      * \param num_tags Number of tag handles
00279      * \param tag_values Array of tag values passed as strings; the array will be
00280      *       the same length as that for tag handles however some entries may be
00281      *       NULL indicating that tag should be matched for existence and not value
00282      * \param integ_type Type of integration to perform
00283      * \param num_integ_pts The number of Gaussian integration points to use in each dimension
00284      */
00285     ErrorCode normalize_subset( EntityHandle root_set,
00286                                 const char* norm_tag,
00287                                 Tag* tag_handles,
00288                                 int num_tags,
00289                                 const char** tag_values,
00290                                 Coupler::IntegType integ_type,
00291                                 int num_integ_pts );
00292 
00293     /* \brief Retrieve groups of entities matching tags and values(if present)
00294      * Retrieve a vector of vectors of entity handles matching the
00295      * tags and values.  The entity set passed is used as the search domain.
00296      *
00297      * \param norm_tag Tag containing field data to integrate
00298      * \param entity_sets Pointer to vector of vectors of entity set handles
00299      * \param entity_groups Pointer to vector of vectors of entity handles from each entity set
00300      * \param integ_type Type of integration to perform
00301      * \param num_integ_pts The number of Gaussian integration points to use in each dimension
00302      */
00303     ErrorCode do_normalization( const char* norm_tag,
00304                                 std::vector< std::vector< EntityHandle > >& entity_sets,
00305                                 std::vector< std::vector< EntityHandle > >& entity_groups,
00306                                 Coupler::IntegType integ_type,
00307                                 int num_integ_pts );
00308 
00309     /* \brief Retrieve groups of entities matching tags and values(if present)
00310      * Retrieve a vector of vectors of entity handles matching the
00311      * tags and values.  The entity set passed is used as the search domain.
00312      *
00313      * \param root_set Set from which to search for matching entities
00314      * \param tag_names Array of tag names used to select entities
00315      * \param tag_values Array of tag values used to select entities
00316      * \param num_tags Number of tag names
00317      * \param entity_sets Pointer to vector of vectors of entity set handles found in the search
00318      * \param entity_groups Pointer to vector of vectors of entity handles from each entity set
00319      */
00320     ErrorCode get_matching_entities( EntityHandle root_set,
00321                                      const char** tag_names,
00322                                      const char** tag_values,
00323                                      int num_tags,
00324                                      std::vector< std::vector< EntityHandle > >* entity_sets,
00325                                      std::vector< std::vector< EntityHandle > >* entity_groups );
00326 
00327     /* \brief Retrieve groups of entities matching tags and values(if present)
00328      * Retrieve a vector of vectors of entity handles matching the
00329      * tags and values.  The entity set passed is used as the search domain.
00330      *
00331      * \param root_set Set from which to search for matching entities
00332      * \param tag_handles Array of tag handles used to select entities
00333      * \param tag_values Array of tag values used to select entities
00334      * \param num_tags Number of tag handles
00335      * \param entity_sets Pointer to vector of vectors of entity set handles found in the search
00336      * \param entity_groups Pointer to vector of vectors of entity handles from each entity set
00337      */
00338     ErrorCode get_matching_entities( EntityHandle root_set,
00339                                      Tag* tag_handles,
00340                                      const char** tag_values,
00341                                      int num_tags,
00342                                      std::vector< std::vector< EntityHandle > >* entity_sets,
00343                                      std::vector< std::vector< EntityHandle > >* entity_groups );
00344 
00345     /* \brief Return an array of tuples of tag values for each Entity Set
00346      * A list of n-tuples will be constructed with 1 n-tuple for each Entity Set.
00347      * The n-tuple will have an component for each tag given.  It is assumed all
00348      * of the tags are integer tags.
00349      *
00350      * \param ent_sets Array of Entity Set handles to use for retrieving tag data
00351      * \param num_sets Number of Entity Sets
00352      * \param tag_names Array of tag names
00353      * \param num_tags Number of tag names
00354      * \param tuples The returned tuple_list structure
00355      */
00356     ErrorCode create_tuples( Range& ent_sets, const char** tag_names, unsigned int num_tags, TupleList** tuples );
00357 
00358     /* \brief Return an array of tuples of tag values for each Entity Set
00359      * A list of n-tuples will be constructed with 1 n-tuple for each Entity Set.
00360      * The n-tuple will have an component for each tag given.  It is assumed all
00361      * of the tags are integer tags.
00362      *
00363      * \param ent_sets Array of Entity Set handles to use for retrieving tag data
00364      * \param num_sets Number of Entity Sets
00365      * \param tag_handles Array of tag handles
00366      * \param num_tags Number of tag handles
00367      * \param tuples The returned tuple_list structure
00368      */
00369     ErrorCode create_tuples( Range& ent_sets, Tag* tag_handles, unsigned int num_tags, TupleList** tuples );
00370 
00371     /* \brief Consolidate an array of n-tuples lists into one n-tuple list with no duplicates
00372      * An array of list of n-tuples are consolidated into a single list of n-tuples
00373      * with all duplicates removed.  Only integer columns in the tuple_list are assumed to
00374      * be used.
00375      *
00376      * \param all_tuples Array of tuple_lists to consolidate to one
00377      * \param num_tuples Number of tuple_lists
00378      * \param unique_tuples The consolidated tuple_list with no duplicates
00379      */
00380     ErrorCode consolidate_tuples( TupleList** all_tuples, unsigned int num_tuples, TupleList** unique_tuples );
00381 
00382     /* \brief Calculate integrated field values for groups of entities
00383      * An integrated field value, as defined by the field function,
00384      * is calculated for each group of entities passed in.
00385      *
00386      * \param groups The vector contains vectors of entity handles, each representing a group
00387      * \param integ_vals The integrated field values for each group
00388      * \param norm_tag The tag name of the vertex-based field to be integrated
00389      * \param num_integ_pts The number of Gaussian integration points to use in each dimension
00390      * \param integ_type Type of integration to perform
00391      */
00392     ErrorCode get_group_integ_vals( std::vector< std::vector< EntityHandle > >& groups,
00393                                     std::vector< double >& integ_vals,
00394                                     const char* norm_tag,
00395                                     int num_integ_pts,
00396                                     Coupler::IntegType integ_type );
00397 
00398     /* \brief Apply a normalization factor to group of entities
00399      * Multiply a normalization factor with the value of norm_tag for each vertex
00400      * of each entity in a group.  Save the value back to norm_tag on each vertex.
00401      *
00402      * \param entity_sets The vector contains vectors of entity set handles, each containing the
00403      * members of a group \param norm_factors The normalization factors for each group \param
00404      * norm_tag The tag to be normalized on each group \param integ_type Type of integration to
00405      * perform
00406      */
00407     ErrorCode apply_group_norm_factor( std::vector< std::vector< EntityHandle > >& entity_sets,
00408                                        std::vector< double >& norm_factors,
00409                                        const char* norm_tag,
00410                                        Coupler::IntegType integ_type );
00411 
00412     /*
00413      * this method will look at source (and target sets?) sets, and look for the SEM_DIMS tag
00414      * if it exists, it will trigger a spectral element caching, with the order specified
00415      */
00416     ErrorCode initialize_spectral_elements( EntityHandle rootSource,
00417                                             EntityHandle rootTarget,
00418                                             bool& specSou,
00419                                             bool& specTar );
00420 
00421     /*
00422      * this method will put in an array, interleaved, the points of interest for coupling
00423      * with a target mesh (so where do we need to compute the field of interest)
00424      */
00425     ErrorCode get_gl_points_on_elements( Range& targ_elems, std::vector< double >& vpos, int& numPointsOfInterest );
00426 
00427     /* Get functions */
00428 
00429     inline Interface* mb_impl() const
00430     {
00431         return mbImpl;
00432     }
00433     inline AdaptiveKDTree* my_tree() const
00434     {
00435         return myTree;
00436     }
00437     inline EntityHandle local_root() const
00438     {
00439         return localRoot;
00440     }
00441     inline const std::vector< double >& all_boxes() const
00442     {
00443         return allBoxes;
00444     }
00445     inline ParallelComm* my_pc() const
00446     {
00447         return myPc;
00448     }
00449     inline const Range& target_ents() const
00450     {
00451         return targetEnts;
00452     }
00453     inline int my_id() const
00454     {
00455         return myId;
00456     }
00457     inline const Range& my_range() const
00458     {
00459         return myRange;
00460     }
00461     inline TupleList* mapped_pts() const
00462     {
00463         return mappedPts;
00464     }
00465     inline int num_its() const
00466     {
00467         return numIts;
00468     }
00469     // used for spherical tests
00470     inline void set_spherical( bool arg1 = true )
00471     {
00472         spherical = arg1;
00473     }
00474 
00475   private:
00476     // Given a coordinate position, find all entities containing
00477     // the point and the natural coords in those ents
00478     ErrorCode nat_param( double xyz[3],
00479                          std::vector< EntityHandle >& entities,
00480                          std::vector< CartVect >& nat_coords,
00481                          double epsilon = 0.0 );
00482 
00483     ErrorCode interp_field( EntityHandle elem, CartVect nat_coord, Tag tag, double& field );
00484 
00485     ErrorCode constant_interp( EntityHandle elem, Tag tag, double& field );
00486 
00487     ErrorCode test_local_box( double* xyz,
00488                               int from_proc,
00489                               int remote_index,
00490                               int index,
00491                               bool& point_located,
00492                               double rel_eps = 0.0,
00493                               double abs_eps = 0.0,
00494                               TupleList* tl  = NULL );
00495 
00496     /* \brief MOAB instance
00497      */
00498     Interface* mbImpl;
00499 
00500     /* \brief Kdtree for local mesh
00501      */
00502     AdaptiveKDTree* myTree;
00503 
00504     /* \brief Local root of the kdtree
00505      */
00506     EntityHandle localRoot;
00507 
00508     /* \brief Min/max bounding boxes for all proc tree roots
00509      */
00510     std::vector< double > allBoxes;
00511 
00512     /* \brief ParallelComm object for this coupler
00513      */
00514     ParallelComm* myPc;
00515 
00516     /* \brief Id of this coupler
00517      */
00518     int myId;
00519 
00520     /* \brief Range of source elements
00521      */
00522     Range myRange;
00523 
00524     /* \brief Range of target entities
00525      */
00526     Range targetEnts;
00527 
00528     /* \brief List of locally mapped tuples
00529      * Tuples contain the following:
00530      * n = # mapped points
00531      * vul[i] = local handle of mapped entity
00532      * vr[3*i..3*i+2] = natural coordinates in mapped entity
00533      */
00534     TupleList* mappedPts;
00535 
00536     /* \brief Tuple list of target points and interpolated data
00537      * Tuples contain the following:
00538      * n = # target points
00539      * vi[3*i]   = remote proc mapping target point
00540      * vi[3*i+1] = local index of target point
00541      * vi[3*i+2] = remote index of target point
00542      */
00543     TupleList* targetPts;
00544 
00545     /* \brief Number of iterations of tree building before failing
00546      *
00547      */
00548     int numIts;
00549 
00550     // Entity dimension
00551     int max_dim;
00552 
00553     // A cached spectral element for source and target, separate
00554     // Assume that their number of GL points (order + 1) does not change
00555     // If it does change, we need to reinitialize it
00556     void* _spectralSource;
00557     void* _spectralTarget;
00558     moab::Tag _xm1Tag, _ym1Tag, _zm1Tag;
00559     int _ntot;
00560 
00561     // spherical coupling
00562     bool spherical;
00563 };
00564 
00565 inline ErrorCode Coupler::interpolate( Coupler::Method method,
00566                                        Tag tag,
00567                                        double* interp_vals,
00568                                        TupleList* tl,
00569                                        bool normalize )
00570 {
00571     int num_pts = ( tl ? tl->get_n() : targetPts->get_n() );
00572     return interpolate( &method, &tag, &num_pts, 1, interp_vals, tl, normalize );
00573 }
00574 
00575 }  // namespace moab
00576 
00577 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines