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