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