Coupler.hpp

Go to the documentation of this file.

/**
 * \class moab::Coupler
 * \author Tim Tautges
 *
 * \brief This class couples data between meshes.
 *
 * The coupler interpolates solution data at a set of points. Data
 * being interpolated resides on a source mesh, in a tag.
 * Applications calling this coupler send in entities, usually points
 * or vertices, and receive back the tag value interpolated at those
 * points. Entities in the source mesh containing those points
 * do not have to reside on the same processor.
 *
 * To use, an application should:
 * - instantiate this coupler by calling the constructor collectively
 * on all processors in the communicator
 * - call locate_points, which locates the points to be interpolated and
 * (optionally) caches the results in this object
 * - call interpolate, which does the interpolation
 *
 * Multiple interpolations can be done after locating the points.
 *
 */
#ifndef COUPLER_HPP
#define COUPLER_HPP
 
#include "moab/Range.hpp"
#include "moab/Interface.hpp"
#include "moab/CartVect.hpp"
#include "moab/TupleList.hpp"
 
#include <sstream>
 
namespace moab
{
 
class ParallelComm;
 
class AdaptiveKDTree;
 
class TupleList;
 
class Coupler
{
 public:
 enum Method
 {
 CONSTANT,
 LINEAR_FE,
 QUADRATIC_FE,
 SPECTRAL,
 SPHERICAL
 };
 
 enum IntegType
 {
 VOLUME
 };
 
 /* Constructor
 * Constructor, which also optionally initializes the coupler
 * \param pc ParallelComm object to be used with this coupler, representing the union
 * of processors containing source and target meshes
 * \param local_elems Local elements in the source mesh
 * \param coupler_id Id of this coupler, should be the same over all procs
 * \param init_tree If true, initializes kdtree inside the constructor
 */
 Coupler( Interface* impl,
 ParallelComm* pc,
 Range& local_elems,
 int coupler_id,
 bool init_tree = true,
 int max_ent_dim = 3 );
 
 /* Destructor
 */
 virtual ~Coupler();
 
 /**
 * \brief Initialize the kdtree, locally and across communicator
 */
 ErrorCode initialize_tree();
 
 /* \brief Locate points on the source mesh
 * This function finds the element/processor/natural coordinates for the
 * source mesh element containing each point, optionally storing the results
 * on the target mesh processor. Relative tolerance is compared to bounding
 * box diagonal length. Tolerance is compared to [-1,1] parametric extent
 * in the reference element.
 * \param xyz Point locations (interleaved) being located
 * \param num_points Number of points in xyz
 * \param rel_eps Relative tolerance for the non-linear iteration inside a given element
 * \param abs_eps Absolute tolerance for the non-linear iteration inside a given element
 * \param tl Tuple list containing the results, with each tuple
 * consisting of (p, i), p = proc, i = index on that proc
 * \param store_local If true, stores the tuple list in targetPts
 */
 ErrorCode locate_points( double* xyz,
 unsigned int num_points,
 double rel_eps = 0.0,
 double abs_eps = 0.0,
 TupleList* tl = NULL,
 bool store_local = true );
 
 /* \brief Locate entities on the source mesh
 * This function finds the element/processor/natural coordinates for the
 * source mesh element containing each entity, optionally storing the results
 * on the target mesh processor. Location of each target mesh entity passed in
 * is its centroid (for non-vertices, the avg of its vertex positions).
 * Relative tolerance is compared to bounding
 * box diagonal length. Tolerance is compared to [-1,1] parametric extent
 * in the reference element.
 * \param ents Entities being located
 * \param rel_eps Relative tolerance for the non-linear iteration inside a given element
 * \param abs_eps Absolute tolerance for the non-linear iteration inside a given element
 * \param tl Tuple list containing the results, with each tuple
 * consisting of (p, i), p = proc, i = index on that proc
 * \param store_local If true, stores the tuple list in targetPts
 *
 */
 ErrorCode locate_points( Range& ents,
 double rel_eps = 0.0,
 double abs_eps = 0.0,
 TupleList* tl = NULL,
 bool store_local = true );
 
 /* \brief Interpolate data from the source mesh onto points
 * All entities/points or, if tuple_list is input, only those points
 * are interpolated from the source mesh. Application should
 * allocate enough memory in interp_vals to hold interpolation results.
 *
 * If normalization is requested, technique used depends on the coupling
 * method.
 *
 * \param method Interpolation/normalization method
 * \param tag Tag on source mesh holding data to be interpolated
 * \param interp_vals Memory holding interpolated data
 * \param tl Tuple list of points to be interpolated, in format used by targetPts
 * (see documentation for targetPts below); if NULL, all locations
 * in targetPts are interpolated
 * \param normalize If true, normalization is done according to method
 */
 ErrorCode interpolate( Coupler::Method method,
 Tag tag,
 double* interp_vals,
 TupleList* tl = NULL,
 bool normalize = true );
 
 /* \brief Interpolate data from the source mesh onto points
 * All entities/points or, if tuple_list is input, only those points
 * are interpolated from the source mesh. Application should
 * allocate enough memory in interp_vals to hold interpolation results.
 *
 * If normalization is requested, technique used depends on the coupling
 * method.
 *
 * \param methods Interpolation/normalization method
 * \param tag_name Name of tag on source mesh holding data to be interpolated
 * \param interp_vals Memory holding interpolated data
 * \param tl Tuple list of points to be interpolated, in format used by targetPts
 * (see documentation for targetPts below); if NULL, all locations
 * in targetPts are interpolated
 * \param normalize If true, normalization is done according to method
 */
 ErrorCode interpolate( Coupler::Method method,
 const std::string& tag_name,
 double* interp_vals,
 TupleList* tl = NULL,
 bool normalize = true );
 
 /* \brief Interpolate data from multiple tags
 * All entities/points or, if tuple_list is input, only those points
 * are interpolated from the source mesh. Application should
 * allocate enough memory in interp_vals to hold interpolation results.
 *
 * In this variant, multiple tags, possibly with multiple interpolation
 * methods, are specified. Sum of values in points_per_method should be
 * the number of points in tl or, if NULL, targetPts.
 *
 * If normalization is requested, technique used depends on the coupling
 * method.
 *
 * \param methods Vector of Interpolation/normalization methods
 * \param tag_names Names of tag being interpolated for each method
 * \param points_per_method Number of points for each method
 * \param num_methods Length of vectors in previous 3 arguments
 * \param interp_vals Memory holding interpolated data
 * \param tl Tuple list of points to be interpolated; if NULL, all locations
 * stored in this object are interpolated
 * \param normalize If true, normalization is done according to method
 */
 ErrorCode interpolate( Coupler::Method* methods,
 const std::string* tag_names,
 int* points_per_method,
 int num_methods,
 double* interp_vals,
 TupleList* tl = NULL,
 bool normalize = true );
 
 /* \brief Interpolate data from multiple tags
 * All entities/points or, if tuple_list is input, only those points
 * are interpolated from the source mesh. Application should
 * allocate enough memory in interp_vals to hold interpolation results.
 *
 * In this variant, multiple tags, possibly with multiple interpolation
 * methods, are specified. Sum of values in points_per_method should be
 * the number of points in tl or, if NULL, targetPts.
 *
 * If normalization is requested, technique used depends on the coupling
 * method.
 *
 * \param methods Vector of Interpolation/normalization methods
 * \param tag_names Names of tag being interpolated for each method
 * \param points_per_method Number of points for each method
 * \param num_methods Length of vectors in previous 3 arguments
 * \param interp_vals Memory holding interpolated data
 * \param tl Tuple list of points to be interpolated; if NULL, all locations
 * stored in this object are interpolated
 * \param normalize If true, normalization is done according to method
 */
 ErrorCode interpolate( Coupler::Method* methods,
 Tag* tag_names,
 int* points_per_method,
 int num_methods,
 double* interp_vals,
 TupleList* tl = NULL,
 bool normalize = true );
 
 /* \brief Normalize a field over an entire mesh
 * A field existing on the vertices of elements of a mesh is integrated
 * over all elements in the mesh. The integrated value is normalized
 * and the normalization factor is saved to a new tag
 * on the mesh entity set.
 *
 * \param root_set Entity Set representing the entire mesh
 * \param norm_tag Tag containing field data to integrate
 * \param integ_type Type of integration to perform
 * \param num_integ_pts The number of Gaussian integration points to use in each dimension
 */
 ErrorCode normalize_mesh( EntityHandle root_set,
 const char* norm_tag,
 Coupler::IntegType integ_type,
 int num_integ_pts );
 
 /* \brief Normalize a field over subsets of entities
 * A field existing on the vertices of elements of a mesh is integrated
 * over subsets of elements identified by the tags and values. The integrated
 * values are normalized and the normalization factor is saved to a new tag
 * on the entity sets which contain the elements of a subset.
 *
 * \param root_set Entity Set from the mesh from which to select subsets
 * \param norm_tag Tag containing field data to integrate
 * \param tag_names Array of tag names used for selecting element subsets
 * \param num_tags Number of tag names
 * \param tag_values Array of tag values passed as strings; the array will be
 * the same length as that for tag names however some entries may be
 * NULL indicating that tag should be matched for existence and not value
 * \param integ_type Type of integration to perform
 * \param num_integ_pts The number of Gaussian integration points to use in each dimension
 */
 ErrorCode normalize_subset( EntityHandle root_set,
 const char* norm_tag,
 const char** tag_names,
 int num_tags,
 const char** tag_values,
 Coupler::IntegType integ_type,
 int num_integ_pts );
 
 /* \brief Normalize a field over subsets of entities
 * A field existing on the vertices of elements of a mesh is integrated
 * over subsets of elements identified by the tags and values. The integrated
 * values are normalized and the normalization factor is saved to a new tag
 * on the entity sets which contain the elements of a subset.
 *
 * \param root_set Entity Set from the mesh from which to select subsets
 * \param norm_tag Tag containing field data to integrate
 * \param tag_handles Array of tag handles used for selecting element subsets
 * \param num_tags Number of tag handles
 * \param tag_values Array of tag values passed as strings; the array will be
 * the same length as that for tag handles however some entries may be
 * NULL indicating that tag should be matched for existence and not value
 * \param integ_type Type of integration to perform
 * \param num_integ_pts The number of Gaussian integration points to use in each dimension
 */
 ErrorCode normalize_subset( EntityHandle root_set,
 const char* norm_tag,
 Tag* tag_handles,
 int num_tags,
 const char** tag_values,
 Coupler::IntegType integ_type,
 int num_integ_pts );
 
 /* \brief Retrieve groups of entities matching tags and values(if present)
 * Retrieve a vector of vectors of entity handles matching the
 * tags and values. The entity set passed is used as the search domain.
 *
 * \param norm_tag Tag containing field data to integrate
 * \param entity_sets Pointer to vector of vectors of entity set handles
 * \param entity_groups Pointer to vector of vectors of entity handles from each entity set
 * \param integ_type Type of integration to perform
 * \param num_integ_pts The number of Gaussian integration points to use in each dimension
 */
 ErrorCode do_normalization( const char* norm_tag,
 std::vector< std::vector< EntityHandle > >& entity_sets,
 std::vector< std::vector< EntityHandle > >& entity_groups,
 Coupler::IntegType integ_type,
 int num_integ_pts );
 
 /* \brief Retrieve groups of entities matching tags and values(if present)
 * Retrieve a vector of vectors of entity handles matching the
 * tags and values. The entity set passed is used as the search domain.
 *
 * \param root_set Set from which to search for matching entities
 * \param tag_names Array of tag names used to select entities
 * \param tag_values Array of tag values used to select entities
 * \param num_tags Number of tag names
 * \param entity_sets Pointer to vector of vectors of entity set handles found in the search
 * \param entity_groups Pointer to vector of vectors of entity handles from each entity set
 */
 ErrorCode get_matching_entities( EntityHandle root_set,
 const char** tag_names,
 const char** tag_values,
 int num_tags,
 std::vector< std::vector< EntityHandle > >* entity_sets,
 std::vector< std::vector< EntityHandle > >* entity_groups );
 
 /* \brief Retrieve groups of entities matching tags and values(if present)
 * Retrieve a vector of vectors of entity handles matching the
 * tags and values. The entity set passed is used as the search domain.
 *
 * \param root_set Set from which to search for matching entities
 * \param tag_handles Array of tag handles used to select entities
 * \param tag_values Array of tag values used to select entities
 * \param num_tags Number of tag handles
 * \param entity_sets Pointer to vector of vectors of entity set handles found in the search
 * \param entity_groups Pointer to vector of vectors of entity handles from each entity set
 */
 ErrorCode get_matching_entities( EntityHandle root_set,
 Tag* tag_handles,
 const char** tag_values,
 int num_tags,
 std::vector< std::vector< EntityHandle > >* entity_sets,
 std::vector< std::vector< EntityHandle > >* entity_groups );
 
 /* \brief Return an array of tuples of tag values for each Entity Set
 * A list of n-tuples will be constructed with 1 n-tuple for each Entity Set.
 * The n-tuple will have an component for each tag given. It is assumed all
 * of the tags are integer tags.
 *
 * \param ent_sets Array of Entity Set handles to use for retrieving tag data
 * \param num_sets Number of Entity Sets
 * \param tag_names Array of tag names
 * \param num_tags Number of tag names
 * \param tuples The returned tuple_list structure
 */
 ErrorCode create_tuples( Range& ent_sets, const char** tag_names, unsigned int num_tags, TupleList** tuples );
 
 /* \brief Return an array of tuples of tag values for each Entity Set
 * A list of n-tuples will be constructed with 1 n-tuple for each Entity Set.
 * The n-tuple will have an component for each tag given. It is assumed all
 * of the tags are integer tags.
 *
 * \param ent_sets Array of Entity Set handles to use for retrieving tag data
 * \param num_sets Number of Entity Sets
 * \param tag_handles Array of tag handles
 * \param num_tags Number of tag handles
 * \param tuples The returned tuple_list structure
 */
 ErrorCode create_tuples( Range& ent_sets, Tag* tag_handles, unsigned int num_tags, TupleList** tuples );
 
 /* \brief Consolidate an array of n-tuples lists into one n-tuple list with no duplicates
 * An array of list of n-tuples are consolidated into a single list of n-tuples
 * with all duplicates removed. Only integer columns in the tuple_list are assumed to
 * be used.
 *
 * \param all_tuples Array of tuple_lists to consolidate to one
 * \param num_tuples Number of tuple_lists
 * \param unique_tuples The consolidated tuple_list with no duplicates
 */
 ErrorCode consolidate_tuples( TupleList** all_tuples, unsigned int num_tuples, TupleList** unique_tuples );
 
 /* \brief Calculate integrated field values for groups of entities
 * An integrated field value, as defined by the field function,
 * is calculated for each group of entities passed in.
 *
 * \param groups The vector contains vectors of entity handles, each representing a group
 * \param integ_vals The integrated field values for each group
 * \param norm_tag The tag name of the vertex-based field to be integrated
 * \param num_integ_pts The number of Gaussian integration points to use in each dimension
 * \param integ_type Type of integration to perform
 */
 ErrorCode get_group_integ_vals( std::vector< std::vector< EntityHandle > >& groups,
 std::vector< double >& integ_vals,
 const char* norm_tag,
 int num_integ_pts,
 Coupler::IntegType integ_type );
 
 /* \brief Apply a normalization factor to group of entities
 * Multiply a normalization factor with the value of norm_tag for each vertex
 * of each entity in a group. Save the value back to norm_tag on each vertex.
 *
 * \param entity_sets The vector contains vectors of entity set handles, each containing the
 * members of a group \param norm_factors The normalization factors for each group \param
 * norm_tag The tag to be normalized on each group \param integ_type Type of integration to
 * perform
 */
 ErrorCode apply_group_norm_factor( std::vector< std::vector< EntityHandle > >& entity_sets,
 std::vector< double >& norm_factors,
 const char* norm_tag,
 Coupler::IntegType integ_type );
 
 /*
 * this method will look at source (and target sets?) sets, and look for the SEM_DIMS tag
 * if it exists, it will trigger a spectral element caching, with the order specified
 */
 ErrorCode initialize_spectral_elements( EntityHandle rootSource,
 EntityHandle rootTarget,
 bool& specSou,
 bool& specTar );
 
 /*
 * this method will put in an array, interleaved, the points of interest for coupling
 * with a target mesh (so where do we need to compute the field of interest)
 */
 ErrorCode get_gl_points_on_elements( Range& targ_elems, std::vector< double >& vpos, int& numPointsOfInterest );
 
 /* Get functions */
 
 inline Interface* mb_impl() const
 {
 return mbImpl;
 }
 inline AdaptiveKDTree* my_tree() const
 {
 return myTree;
 }
 inline EntityHandle local_root() const
 {
 return localRoot;
 }
 inline const std::vector< double >& all_boxes() const
 {
 return allBoxes;
 }
 inline ParallelComm* my_pc() const
 {
 return myPc;
 }
 inline const Range& target_ents() const
 {
 return targetEnts;
 }
 inline int my_id() const
 {
 return myId;
 }
 inline const Range& my_range() const
 {
 return myRange;
 }
 inline TupleList* mapped_pts() const
 {
 return mappedPts;
 }
 inline int num_its() const
 {
 return numIts;
 }
 // used for spherical tests
 inline void set_spherical( bool arg1 = true )
 {
 spherical = arg1;
 }
 
 private:
 // Given a coordinate position, find all entities containing
 // the point and the natural coords in those ents
 ErrorCode nat_param( double xyz[3],
 std::vector< EntityHandle >& entities,
 std::vector< CartVect >& nat_coords,
 double epsilon = 0.0 );
 
 ErrorCode interp_field( EntityHandle elem, CartVect nat_coord, Tag tag, double& field );
 
 ErrorCode constant_interp( EntityHandle elem, Tag tag, double& field );
 
 ErrorCode test_local_box( double* xyz,
 int from_proc,
 int remote_index,
 int index,
 bool& point_located,
 double rel_eps = 0.0,
 double abs_eps = 0.0,
 TupleList* tl = NULL );
 
 /* \brief MOAB instance
 */
 Interface* mbImpl;
 
 /* \brief Kdtree for local mesh
 */
 AdaptiveKDTree* myTree;
 
 /* \brief Local root of the kdtree
 */
 EntityHandle localRoot;
 
 /* \brief Min/max bounding boxes for all proc tree roots
 */
 std::vector< double > allBoxes;
 
 /* \brief ParallelComm object for this coupler
 */
 ParallelComm* myPc;
 
 /* \brief Id of this coupler
 */
 int myId;
 
 /* \brief Range of source elements
 */
 Range myRange;
 
 /* \brief Range of target entities
 */
 Range targetEnts;
 
 /* \brief List of locally mapped tuples
 * Tuples contain the following:
 * n = # mapped points
 * vul[i] = local handle of mapped entity
 * vr[3*i..3*i+2] = natural coordinates in mapped entity
 */
 TupleList* mappedPts;
 
 /* \brief Tuple list of target points and interpolated data
 * Tuples contain the following:
 * n = # target points
 * vi[3*i] = remote proc mapping target point
 * vi[3*i+1] = local index of target point
 * vi[3*i+2] = remote index of target point
 */
 TupleList* targetPts;
 
 /* \brief Number of iterations of tree building before failing
 *
 */
 int numIts;
 
 // Entity dimension
 int max_dim;
 
 // A cached spectral element for source and target, separate
 // Assume that their number of GL points (order + 1) does not change
 // If it does change, we need to reinitialize it
 void* _spectralSource;
 void* _spectralTarget;
 moab::Tag _xm1Tag, _ym1Tag, _zm1Tag;
 int _ntot;
 
 // spherical coupling
 bool spherical;
};
 
inline ErrorCode Coupler::interpolate( Coupler::Method method,
 Tag tag,
 double* interp_vals,
 TupleList* tl,
 bool normalize )
{
 int num_pts = ( tl ? tl->get_n() : targetPts->get_n() );
 return interpolate( &method, &tag, &num_pts, 1, interp_vals, tl, normalize );
}
 
} // namespace moab
 
#endif