NestedRefine.hpp

Go to the documentation of this file.

/*! \file NestedRefine.hpp
 * This class implements the generation of a hierarchy of meshes via uniform refinement from an
 * input mesh. The internal upper bound on the number of levels is set to 20.
 */

#ifndef NESTED_REFINE_HPP
#define NESTED_REFINE_HPP

#include "moab/MOABConfig.h"
#include "moab/Range.hpp"
#include "moab/CN.hpp"
#include <map>
#include <set>

namespace moab
{

#define MAX_DEGREE    3
#define MAX_VERTS     64
#define MAX_CHILDRENS 27
#define MAX_HE        12
#define MAX_HF        6
#define MAX_CONN      8
#define MAX_VHF       20
#define MAX_LEVELS    20
#define MAX_PROCS     64

// Forward Declarations
class Core;
class HalfFacetRep;
class ParallelComm;
class CpuTimer;

class NestedRefine
{

  public:
    NestedRefine( Core* impl, ParallelComm* comm = 0, EntityHandle rset = 0 );

    ~NestedRefine();

    ErrorCode initialize();

    //! \brief Generate a mesh hierarchy.
    /** Given a mesh, generate a sequence of meshes via uniform refinement. The inputs are: a) an
     * array(level_degrees) storing the degrees which will be used to refine the previous level mesh
     * to generate a new level and b) the total number of levels(should be same length as that of
     * the array in a). Each mesh level in the hierarchy are stored in different meshsets whose
     * handles are returned after the hierarchy generation. These handles can be used to work with a
     * specific mesh level. \param level_degrees Integer array storing the degrees used in each
     * level. \param num_level The total number of levels in the hierarchy. \param hm_set
     * EntityHandle STL vector that returns the handles of the sets created for each mesh level.
     */

    ErrorCode generate_mesh_hierarchy( int num_level,
                                       int* level_degrees,
                                       std::vector< EntityHandle >& level_sets,
                                       bool optimize = false );

    //! Given an entity and its level, return its connectivity.
    /** Given an entity at a certain level, it finds the connectivity via direct access to a stored
     * internal pointer to the memory to connectivity sequence for the given level. \param ent
     * EntityHandle of the entity \param level Integer level of the entity for which connectivity is
     * requested \param conn std::vector returning the connectivity of the entity
     */

    ErrorCode get_connectivity( EntityHandle ent, int level, std::vector< EntityHandle >& conn );

    //! Given a vector of vertices and their level, return its coordinates.
    /** Given a vector of vertices at a certain level, it finds the coordinates via direct access to
     * a stored internal pointer to the memory to coordinate sequence for the given level. \param
     * verts std::vector of the entity handles of the vertices \param num_verts The number of
     * vertices \param level Integer level of the entity for which connectivity is requested \param
     * coords double pointer returning the coordinates of the vertices
     */

    ErrorCode get_coordinates( EntityHandle* verts, int num_verts, int level, double* coords );

    //! Get the adjacencies associated with an entity.
    /** Given an entity of dimension <em>d</em>, gather all the adjacent <em>D</em> dimensional
     * entities where <em>D >, = , < d </em>.
     *
     * \param source_entity EntityHandle to which adjacent entities have to be found.
     * \param target_dimension Int Dimension of the desired adjacent entities.
     * \param target_entities Vector in which the adjacent EntityHandle are returned.
     */

    ErrorCode get_adjacencies( const EntityHandle source_entity,
                               const unsigned int target_dimension,
                               std::vector< EntityHandle >& target_entities );

    // Interlevel parent-child or vice-versa queries
    /** Given an entity from a certain level, it returns a pointer to its parent at the requested
     * parent level. NOTE: This query does not support vertices.
     *
     * \param child EntityHandle of the entity whose parent is requested
     * \param child_level Mesh level where the child exists
     * \param parent_level Mesh level from which parent is requested
     * \param parent Pointer to the parent in the requested parent_level
     */

    ErrorCode child_to_parent( EntityHandle child, int child_level, int parent_level, EntityHandle* parent );

    /** Given an entity from a certain level, it returns a std::vector of all its children from the
     * requested child level. NOTE: This query does not support vertices.
     *
     * \param parent EntityHandle of the entity whose children in subsequent level is requested
     * \param parent_level Mesh level where the parent exists
     * \param child_level Mesh level from which its children are requested
     * \param children Vector containing all childrens from the requested child_level
     */

    ErrorCode parent_to_child( EntityHandle parent,
                               int parent_level,
                               int child_level,
                               std::vector< EntityHandle >& children );

    /** Given a vertex from a certain level, it returns a std::vector of all entities from any
     * previous levels that contains it.
     *
     * \param vertex EntityHandle of the vertex
     * \param vert_level Mesh level of the vertex
     * \param parent_level Mesh level from which entities containing vertex is requested
     * \param incident_entities Vector containing entities from the parent level incident on the
     * vertex
     */

    ErrorCode vertex_to_entities_up( EntityHandle vertex,
                                     int vert_level,
                                     int parent_level,
                                     std::vector< EntityHandle >& incident_entities );

    /** Given a vertex from a certain level, it returns a std::vector of all children entities of
     * incident entities to vertex from any subsequent levels
     *
     * \param vertex EntityHandle of the vertex
     * \param vert_level Mesh level of the vertex
     * \param child_level Mesh level from which child entities are requested
     * \param incident_entities Vector containing entities from the child level
     */

    ErrorCode vertex_to_entities_down( EntityHandle vertex,
                                       int vert_level,
                                       int child_level,
                                       std::vector< EntityHandle >& incident_entities );

    ErrorCode get_vertex_duplicates( EntityHandle vertex, int level, EntityHandle& dupvertex );

    /** Given an entity at a certain level, it returns a boolean value true if it lies on the domain
     * boundary. \param entity
     */

    bool is_entity_on_boundary( const EntityHandle& entity );

    ErrorCode exchange_ghosts( std::vector< EntityHandle >& lsets, int num_glayers );

    ErrorCode update_special_tags( int level, EntityHandle& lset );

    struct codeperf
    {
        double tm_total;
        double tm_refine;
        double tm_resolve;
    };

    codeperf timeall;

  protected:
    Core* mbImpl;
    ParallelComm* pcomm;
    HalfFacetRep* ahf;
    CpuTimer* tm;
    EntityHandle _rset;

    Range _inverts, _inedges, _infaces, _incells;

    EntityType elementype;
    int meshdim, nlevels;
    int level_dsequence[MAX_LEVELS];
    std::map< int, int > deg_index;
    bool hasghost;

    /*! \struct refPatterns
     * Refinement patterns w.r.t the reference element. It consists of a locally indexed vertex list
     * along with their natural coordinates, the connectivity of the subdivided entities with local
     * indices, their local AHF maps along with other helper fields to aid in general book keeping
     * such as avoiding vertex duplication during refinement. The entity and degree specific values
     * are stored in the Templates.hpp. \sa Templates.hpp
     */

    //! refPatterns
    struct refPatterns
    {
        //! Number of new vertices on edge
        short int nv_edge;
        //! Number of new vertices on face, does not include those on edge
        short int nv_face;
        //! Number of new vertices in cell
        short int nv_cell;
        //! Total number of new vertices per entity
        short int total_new_verts;
        //! Total number of new child entities
        short int total_new_ents;
        //! Lower and upper indices of the new vertices
        int vert_index_bnds[2];
        //! Natural coordinates of the new vertices w.r.t reference
        double vert_nat_coord[MAX_VERTS][3];
        //! Connectivity of the new entities
        int ents_conn[MAX_CHILDRENS][MAX_CONN];
        //! Vertex to half-facet map of the new vertices
        int v2hf[MAX_VERTS][2];
        //! Opposite half-facet map of the new entities
        int ents_opphfs[MAX_CHILDRENS][2 * MAX_CONN];
        //! Helper: storing the local ids of vertices on each local edge
        int vert_on_edges[MAX_HE][MAX_VHF];
        //!  Helper: storing local ids of verts on each local face, doesnt include those on edges of
        //!  the face
        int vert_on_faces[MAX_HF][MAX_VHF];
        //! Helper: stores child half-facets incident on parent half-facet. First column contain the
        //! number of such children
        int ents_on_pent[MAX_HF][MAX_CHILDRENS];
        //! Helper: stores child ents incident on new verts on edge.
        // Each triad in the column consists of :
        // 1) a local child incident on the vertex on the edge
        // 2) the local face id from the child
        // 3) the local vertex id wrt the child connectivity
        int ents_on_vedge[MAX_HE][MAX_VHF * 3];
    };
    //! refPatterns

    static const refPatterns refTemplates[9][MAX_DEGREE];

    //! Helper
    struct intFEdge
    {
        //! Number of edges interior to a face
        short int nie;
        //! Local connectivity of the interior edges
        int ieconn[12][2];
    };
    static const intFEdge intFacEdg[2][2];

    int get_index_from_degree( int degree );

    // HM Storage Helper
    struct level_memory
    {
        int num_verts, num_edges, num_faces, num_cells;
        EntityHandle start_vertex, start_edge, start_face, start_cell;
        std::vector< double* > coordinates;
        EntityHandle *edge_conn, *face_conn, *cell_conn;
        Range verts, edges, faces, cells;
    };

    level_memory level_mesh[MAX_LEVELS];

    // Basic Functions

    // Estimate and create storage for the levels
    ErrorCode estimate_hm_storage( EntityHandle set, int level_degree, int cur_level, int hmest[4] );
    ErrorCode create_hm_storage_single_level( EntityHandle* set, int cur_level, int estL[4] );

    // Generate HM : Construct the hierarchical mesh: 1D, 2D, 3D
    ErrorCode generate_hm( int* level_degrees, int num_level, EntityHandle* hm_set, bool optimize );
    ErrorCode construct_hm_entities( int cur_level, int deg );
    ErrorCode construct_hm_1D( int cur_level, int deg );
    ErrorCode construct_hm_1D( int cur_level, int deg, EntityType type, std::vector< EntityHandle >& trackverts );
    ErrorCode construct_hm_2D( int cur_level, int deg );
    ErrorCode construct_hm_2D( int cur_level,
                               int deg,
                               EntityType type,
                               std::vector< EntityHandle >& trackvertsC_edg,
                               std::vector< EntityHandle >& trackvertsF );
    ErrorCode construct_hm_3D( int cur_level, int deg );

    ErrorCode subdivide_cells( EntityType type, int cur_level, int deg );
    ErrorCode subdivide_tets( int cur_level, int deg );

    // General helper functions
    ErrorCode copy_vertices_from_prev_level( int cur_level );
    ErrorCode count_subentities( EntityHandle set, int cur_level, int* nedges, int* nfaces );
    ErrorCode get_octahedron_corner_coords( int cur_level, int deg, EntityHandle* vbuffer, double* ocoords );
    int find_shortest_diagonal_octahedron( int cur_level, int deg, EntityHandle* vbuffer );
    int get_local_vid( EntityHandle vid, EntityHandle ent, int level );

    // Book-keeping functions
    ErrorCode update_tracking_verts( EntityHandle cid,
                                     int cur_level,
                                     int deg,
                                     std::vector< EntityHandle >& trackvertsC_edg,
                                     std::vector< EntityHandle >& trackvertsC_face,
                                     EntityHandle* vbuffer );
    ErrorCode reorder_indices( int cur_level,
                               int deg,
                               EntityHandle cell,
                               int lfid,
                               EntityHandle sib_cell,
                               int sib_lfid,
                               int index,
                               int* id_sib );
    ErrorCode reorder_indices( int deg,
                               EntityHandle* face1_conn,
                               EntityHandle* face2_conn,
                               int nvF,
                               std::vector< int >& lemap,
                               std::vector< int >& vidx,
                               int* leorient = NULL );
    ErrorCode reorder_indices( int deg, int nvF, int comb, int* childfid_map );
    ErrorCode reorder_indices( EntityHandle* face1_conn,
                               EntityHandle* face2_conn,
                               int nvF,
                               int* conn_map,
                               int& comb,
                               int* orient = NULL );
    ErrorCode get_lid_inci_child( EntityType type,
                                  int deg,
                                  int lfid,
                                  int leid,
                                  std::vector< int >& child_ids,
                                  std::vector< int >& child_lvids );

    // Permutation matrices
    struct pmat
    {
        short int num_comb;           // Number of combinations
        int comb[MAX_HE][MAX_HE];     // Combinations
        int lemap[MAX_HE][MAX_HE];    // Local edge map
        int orient[MAX_HE];           // Orientation
        int porder2[MAX_HE][MAX_HE];  // Permuted order degree 2
        int porder3[MAX_HE][MAX_HE];  // Permuted order degree 3
    };

    static const pmat permutation[2];

    // Print functions
    ErrorCode print_maps_1D( int level );
    ErrorCode print_maps_2D( int level, EntityType type );
    ErrorCode print_maps_3D( int level, EntityType type );

    // Coordinates
    ErrorCode compute_coordinates( int cur_level,
                                   int deg,
                                   EntityType type,
                                   EntityHandle* vbuffer,
                                   int vtotal,
                                   double* corner_coords,
                                   std::vector< int >& vflag,
                                   int nverts_prev );

    // Update the ahf maps

    ErrorCode update_local_ahf( int deg, EntityType type, EntityHandle* vbuffer, EntityHandle* ent_buffer, int etotal );

    ErrorCode update_local_ahf( int deg,
                                EntityType type,
                                int pat_id,
                                EntityHandle* vbuffer,
                                EntityHandle* ent_buffer,
                                int etotal );

    //  ErrorCode update_global_ahf(EntityType type, int cur_level, int deg);

    //  ErrorCode update_global_ahf(int cur_level, int deg, std::vector<int> &pattern_ids);

    ErrorCode update_global_ahf( EntityType type, int cur_level, int deg, std::vector< int >* pattern_ids = NULL );

    ErrorCode update_global_ahf_1D( int cur_level, int deg );

    ErrorCode update_global_ahf_1D_sub( int cur_level, int deg );

    ErrorCode update_ahf_1D( int cur_level );

    ErrorCode update_global_ahf_2D( int cur_level, int deg );

    ErrorCode update_global_ahf_2D_sub( int cur_level, int deg );

    ErrorCode update_global_ahf_3D( int cur_level, int deg, std::vector< int >* pattern_ids = NULL );

    //    ErrorCode update_global_ahf_3D(int cur_level, int deg);

    //    ErrorCode update_global_ahf_3D(int cur_level, int deg, std::vector<int> &pattern_ids);

    /** Boundary extraction functions
     * Given a vertex at a certain level, it returns a boolean value true if it lies on the domain
     * boundary. Note: This is a specialization of the NestedRefine::is_entity_on_boundary function
     * and applies only to vertex queries. \param entity
     */
    bool is_vertex_on_boundary( const EntityHandle& entity );
    bool is_edge_on_boundary( const EntityHandle& entity );
    bool is_face_on_boundary( const EntityHandle& entity );
    bool is_cell_on_boundary( const EntityHandle& entity );

    // ErrorCode find_skin_faces(EntityHandle set, int level, int nskinF);

    /** Parallel communication routines
     * We implement two strategies to resolve the shared entities of the newly created entities.
     * The first strategy is to use the existing parallel merge capability which essentially uses
     * a coordinate-based matching of vertices and subsequently the entity handles through
     * their connectivities. The second strategy is an optimized and a new algorithm. It uses
     * the existing shared information from the coarse entities and propagates the parallel
     *  information appropriately.
     */

    // Send/Recv Buffer storage
    /*   struct pbuffer{
         int rank;
         std::vector<int> msgsize;
         std::vector<EntityHandle> localBuff;
         std::vector<EntityHandle> remoteBuff;
       };

       pbuffer commBuffers[MAX_PROCS];

       //Parallel tag values
       struct parinfo{
         std::multimap<EntityHandle, int> remoteps;
         std::multimap<EntityHandle, EntityHandle> remotehs;
       };

       parinfo parallelInfo[MAX_LEVELS];*/

#ifdef MOAB_HAVE_MPI

    ErrorCode resolve_shared_ents_parmerge( int level, EntityHandle levelset );
    ErrorCode resolve_shared_ents_opt( EntityHandle* hm_set, int num_levels );

    ErrorCode collect_shared_entities_by_dimension( Range sharedEnts, Range& allEnts );
    ErrorCode collect_FList( int to_proc,
                             Range faces,
                             std::vector< EntityHandle >& FList,
                             std::vector< EntityHandle >& RList );
    ErrorCode collect_EList( int to_proc,
                             Range edges,
                             std::vector< EntityHandle >& EList,
                             std::vector< EntityHandle >& RList );
    ErrorCode collect_VList( int to_proc,
                             Range verts,
                             std::vector< EntityHandle >& VList,
                             std::vector< EntityHandle >& RList );

    ErrorCode decipher_remote_handles( std::vector< int >& sharedprocs,
                                       std::vector< std::vector< int > >& auxinfo,
                                       std::vector< std::vector< EntityHandle > >& localbuffers,
                                       std::vector< std::vector< EntityHandle > >& remotebuffers,
                                       std::multimap< EntityHandle, int >& remProcs,
                                       std::multimap< EntityHandle, EntityHandle >& remHandles );

    ErrorCode decipher_remote_handles_face( int shared_proc,
                                            int numfaces,
                                            std::vector< EntityHandle >& localFaceList,
                                            std::vector< EntityHandle >& remFaceList,
                                            std::multimap< EntityHandle, int >& remProcs,
                                            std::multimap< EntityHandle, EntityHandle >& remHandles );

    ErrorCode decipher_remote_handles_edge( int shared_proc,
                                            int numedges,
                                            std::vector< EntityHandle >& localEdgeList,
                                            std::vector< EntityHandle >& remEdgeList,
                                            std::multimap< EntityHandle, int >& remProcs,
                                            std::multimap< EntityHandle, EntityHandle >& remHandles );

    ErrorCode decipher_remote_handles_vertex( int shared_proc,
                                              int numverts,
                                              std::vector< EntityHandle >& localVertexList,
                                              std::vector< EntityHandle >& remVertexList,
                                              std::multimap< EntityHandle, int >& remProcs,
                                              std::multimap< EntityHandle, EntityHandle >& remHandles );

    ErrorCode update_parallel_tags( std::multimap< EntityHandle, int >& remProcs,
                                    std::multimap< EntityHandle, EntityHandle >& remHandles );

    ErrorCode get_data_from_buff( int dim,
                                  int type,
                                  int level,
                                  int entityidx,
                                  int nentities,
                                  std::vector< EntityHandle >& buffer,
                                  std::vector< EntityHandle >& data );

    bool check_for_parallelinfo( EntityHandle entity, int proc, std::multimap< EntityHandle, int >& remProcs );

    ErrorCode check_for_parallelinfo( EntityHandle entity,
                                      int proc,
                                      std::multimap< EntityHandle, EntityHandle >& remHandles,
                                      std::multimap< EntityHandle, int >& remProcs,
                                      EntityHandle& rhandle );

#endif
};

}  // namespace moab
#endif