MOAB: Mesh Oriented datABase  (version 5.2.1)
NestedRefine.hpp
Go to the documentation of this file.
00001 /*! \file NestedRefine.hpp
00002  * This class implements the generation of a hierarchy of meshes via uniform refinement from an
00003  * input mesh. The internal upper bound on the number of levels is set to 20.
00004  */
00005 
00006 #ifndef NESTED_REFINE_HPP
00007 #define NESTED_REFINE_HPP
00008 
00009 #include "moab/MOABConfig.h"
00010 #include "moab/Range.hpp"
00011 #include "moab/CN.hpp"
00012 #include <map>
00013 #include <set>
00014 
00015 namespace moab
00016 {
00017 
00018 #define MAX_DEGREE    3
00019 #define MAX_VERTS     64
00020 #define MAX_CHILDRENS 27
00021 #define MAX_HE        12
00022 #define MAX_HF        6
00023 #define MAX_CONN      8
00024 #define MAX_VHF       20
00025 #define MAX_LEVELS    20
00026 #define MAX_PROCS     64
00027 
00028 // Forward Declarations
00029 class Core;
00030 class HalfFacetRep;
00031 class ParallelComm;
00032 class CpuTimer;
00033 
00034 class NestedRefine
00035 {
00036 
00037   public:
00038     NestedRefine( Core* impl, ParallelComm* comm = 0, EntityHandle rset = 0 );
00039 
00040     ~NestedRefine();
00041 
00042     ErrorCode initialize();
00043 
00044     //! \brief Generate a mesh hierarchy.
00045     /** Given a mesh, generate a sequence of meshes via uniform refinement. The inputs are: a) an
00046      * array(level_degrees) storing the degrees which will be used to refine the previous level mesh
00047      * to generate a new level and b) the total number of levels(should be same length as that of
00048      * the array in a). Each mesh level in the hierarchy are stored in different meshsets whose
00049      * handles are returned after the hierarchy generation. These handles can be used to work with a
00050      * specific mesh level. \param level_degrees Integer array storing the degrees used in each
00051      * level. \param num_level The total number of levels in the hierarchy. \param hm_set
00052      * EntityHandle STL vector that returns the handles of the sets created for each mesh level.
00053      */
00054 
00055     ErrorCode generate_mesh_hierarchy( int num_level, int* level_degrees, std::vector< EntityHandle >& level_sets,
00056                                        bool optimize = false );
00057 
00058     //! Given an entity and its level, return its connectivity.
00059     /** Given an entity at a certain level, it finds the connectivity via direct access to a stored
00060      * internal pointer to the memory to connectivity sequence for the given level. \param ent
00061      * EntityHandle of the entity \param level Integer level of the entity for which connectivity is
00062      * requested \param conn std::vector returning the connectivity of the entity
00063      */
00064 
00065     ErrorCode get_connectivity( EntityHandle ent, int level, std::vector< EntityHandle >& conn );
00066 
00067     //! Given a vector of vertices and their level, return its coordinates.
00068     /** Given a vector of vertices at a certain level, it finds the coordinates via direct access to
00069      * a stored internal pointer to the memory to coordinate sequence for the given level. \param
00070      * verts std::vector of the entity handles of the vertices \param num_verts The number of
00071      * vertices \param level Integer level of the entity for which connectivity is requested \param
00072      * coords double pointer returning the coordinates of the vertices
00073      */
00074 
00075     ErrorCode get_coordinates( EntityHandle* verts, int num_verts, int level, double* coords );
00076 
00077     //! Get the adjacencies associated with an entity.
00078     /** Given an entity of dimension <em>d</em>, gather all the adjacent <em>D</em> dimensional
00079      * entities where <em>D >, = , < d </em>.
00080      *
00081      * \param source_entity EntityHandle to which adjacent entities have to be found.
00082      * \param target_dimension Int Dimension of the desired adjacent entities.
00083      * \param target_entities Vector in which the adjacent EntityHandle are returned.
00084      */
00085 
00086     ErrorCode get_adjacencies( const EntityHandle source_entity, const unsigned int target_dimension,
00087                                std::vector< EntityHandle >& target_entities );
00088 
00089     // Interlevel parent-child or vice-versa queries
00090     /** Given an entity from a certain level, it returns a pointer to its parent at the requested
00091      * parent level. NOTE: This query does not support vertices.
00092      *
00093      * \param child EntityHandle of the entity whose parent is requested
00094      * \param child_level Mesh level where the child exists
00095      * \param parent_level Mesh level from which parent is requested
00096      * \param parent Pointer to the parent in the requested parent_level
00097      */
00098 
00099     ErrorCode child_to_parent( EntityHandle child, int child_level, int parent_level, EntityHandle* parent );
00100 
00101     /** Given an entity from a certain level, it returns a std::vector of all its children from the
00102      * requested child level. NOTE: This query does not support vertices.
00103      *
00104      * \param parent EntityHandle of the entity whose children in subsequent level is requested
00105      * \param parent_level Mesh level where the parent exists
00106      * \param child_level Mesh level from which its children are requested
00107      * \param children Vector containing all childrens from the requested child_level
00108      */
00109 
00110     ErrorCode parent_to_child( EntityHandle parent, int parent_level, int child_level,
00111                                std::vector< EntityHandle >& children );
00112 
00113     /** Given a vertex from a certain level, it returns a std::vector of all entities from any
00114      * previous levels that contains it.
00115      *
00116      * \param vertex EntityHandle of the vertex
00117      * \param vert_level Mesh level of the vertex
00118      * \param parent_level Mesh level from which entities containing vertex is requested
00119      * \param incident_entities Vector containing entities from the parent level incident on the
00120      * vertex
00121      */
00122 
00123     ErrorCode vertex_to_entities_up( EntityHandle vertex, int vert_level, int parent_level,
00124                                      std::vector< EntityHandle >& incident_entities );
00125 
00126     /** Given a vertex from a certain level, it returns a std::vector of all children entities of
00127      * incident entities to vertex from any subsequent levels
00128      *
00129      * \param vertex EntityHandle of the vertex
00130      * \param vert_level Mesh level of the vertex
00131      * \param child_level Mesh level from which child entities are requested
00132      * \param incident_entities Vector containing entities from the child level
00133      */
00134 
00135     ErrorCode vertex_to_entities_down( EntityHandle vertex, int vert_level, int child_level,
00136                                        std::vector< EntityHandle >& incident_entities );
00137 
00138     ErrorCode get_vertex_duplicates( EntityHandle vertex, int level, EntityHandle& dupvertex );
00139 
00140     /** Given an entity at a certain level, it returns a boolean value true if it lies on the domain
00141      * boundary. \param entity
00142      */
00143 
00144     bool is_entity_on_boundary( const EntityHandle& entity );
00145 
00146     ErrorCode exchange_ghosts( std::vector< EntityHandle >& lsets, int num_glayers );
00147 
00148     ErrorCode update_special_tags( int level, EntityHandle& lset );
00149 
00150     struct codeperf
00151     {
00152         double tm_total;
00153         double tm_refine;
00154         double tm_resolve;
00155     };
00156 
00157     codeperf timeall;
00158 
00159   protected:
00160     Core* mbImpl;
00161     ParallelComm* pcomm;
00162     HalfFacetRep* ahf;
00163     CpuTimer* tm;
00164     EntityHandle _rset;
00165 
00166     Range _inverts, _inedges, _infaces, _incells;
00167 
00168     EntityType elementype;
00169     int meshdim, nlevels;
00170     int level_dsequence[MAX_LEVELS];
00171     std::map< int, int > deg_index;
00172     bool hasghost;
00173 
00174     /*! \struct refPatterns
00175      * Refinement patterns w.r.t the reference element. It consists of a locally indexed vertex list
00176      * along with their natural coordinates, the connectivity of the subdivided entities with local
00177      * indices, their local AHF maps along with other helper fields to aid in general book keeping
00178      * such as avoiding vertex duplication during refinement. The entity and degree specific values
00179      * are stored in the Templates.hpp. \sa Templates.hpp
00180      */
00181 
00182     //! refPatterns
00183     struct refPatterns
00184     {
00185         //! Number of new vertices on edge
00186         short int nv_edge;
00187         //! Number of new vertices on face, does not include those on edge
00188         short int nv_face;
00189         //! Number of new vertices in cell
00190         short int nv_cell;
00191         //! Total number of new vertices per entity
00192         short int total_new_verts;
00193         //! Total number of new child entities
00194         short int total_new_ents;
00195         //! Lower and upper indices of the new vertices
00196         int vert_index_bnds[2];
00197         //! Natural coordinates of the new vertices w.r.t reference
00198         double vert_nat_coord[MAX_VERTS][3];
00199         //! Connectivity of the new entities
00200         int ents_conn[MAX_CHILDRENS][MAX_CONN];
00201         //! Vertex to half-facet map of the new vertices
00202         int v2hf[MAX_VERTS][2];
00203         //! Opposite half-facet map of the new entities
00204         int ents_opphfs[MAX_CHILDRENS][2 * MAX_CONN];
00205         //! Helper: storing the local ids of vertices on each local edge
00206         int vert_on_edges[MAX_HE][MAX_VHF];
00207         //!  Helper: storing local ids of verts on each local face, doesnt include those on edges of
00208         //!  the face
00209         int vert_on_faces[MAX_HF][MAX_VHF];
00210         //! Helper: stores child half-facets incident on parent half-facet. First column contain the
00211         //! number of such children
00212         int ents_on_pent[MAX_HF][MAX_CHILDRENS];
00213         //! Helper: stores child ents incident on new verts on edge.
00214         // Each triad in the column consists of :
00215         // 1) a local child incident on the vertex on the edge
00216         // 2) the local face id from the child
00217         // 3) the local vertex id wrt the child connectivity
00218         int ents_on_vedge[MAX_HE][MAX_VHF * 3];
00219     };
00220     //! refPatterns
00221 
00222     static const refPatterns refTemplates[9][MAX_DEGREE];
00223 
00224     //! Helper
00225     struct intFEdge
00226     {
00227         //! Number of edges interior to a face
00228         short int nie;
00229         //! Local connectivity of the interior edges
00230         int ieconn[12][2];
00231     };
00232     static const intFEdge intFacEdg[2][2];
00233 
00234     int get_index_from_degree( int degree );
00235 
00236     // HM Storage Helper
00237     struct level_memory
00238     {
00239         int num_verts, num_edges, num_faces, num_cells;
00240         EntityHandle start_vertex, start_edge, start_face, start_cell;
00241         std::vector< double* > coordinates;
00242         EntityHandle *edge_conn, *face_conn, *cell_conn;
00243         Range verts, edges, faces, cells;
00244     };
00245 
00246     level_memory level_mesh[MAX_LEVELS];
00247 
00248     // Basic Functions
00249 
00250     // Estimate and create storage for the levels
00251     ErrorCode estimate_hm_storage( EntityHandle set, int level_degree, int cur_level, int hmest[4] );
00252     ErrorCode create_hm_storage_single_level( EntityHandle* set, int cur_level, int estL[4] );
00253 
00254     // Generate HM : Construct the hierarchical mesh: 1D, 2D, 3D
00255     ErrorCode generate_hm( int* level_degrees, int num_level, EntityHandle* hm_set, bool optimize );
00256     ErrorCode construct_hm_entities( int cur_level, int deg );
00257     ErrorCode construct_hm_1D( int cur_level, int deg );
00258     ErrorCode construct_hm_1D( int cur_level, int deg, EntityType type, std::vector< EntityHandle >& trackverts );
00259     ErrorCode construct_hm_2D( int cur_level, int deg );
00260     ErrorCode construct_hm_2D( int cur_level, int deg, EntityType type, std::vector< EntityHandle >& trackvertsC_edg,
00261                                std::vector< EntityHandle >& trackvertsF );
00262     ErrorCode construct_hm_3D( int cur_level, int deg );
00263 
00264     ErrorCode subdivide_cells( EntityType type, int cur_level, int deg );
00265     ErrorCode subdivide_tets( int cur_level, int deg );
00266 
00267     // General helper functions
00268     ErrorCode copy_vertices_from_prev_level( int cur_level );
00269     ErrorCode count_subentities( EntityHandle set, int cur_level, int* nedges, int* nfaces );
00270     ErrorCode get_octahedron_corner_coords( int cur_level, int deg, EntityHandle* vbuffer, double* ocoords );
00271     int find_shortest_diagonal_octahedron( int cur_level, int deg, EntityHandle* vbuffer );
00272     int get_local_vid( EntityHandle vid, EntityHandle ent, int level );
00273 
00274     // Book-keeping functions
00275     ErrorCode update_tracking_verts( EntityHandle cid, int cur_level, int deg,
00276                                      std::vector< EntityHandle >& trackvertsC_edg,
00277                                      std::vector< EntityHandle >& trackvertsC_face, EntityHandle* vbuffer );
00278     ErrorCode reorder_indices( int cur_level, int deg, EntityHandle cell, int lfid, EntityHandle sib_cell, int sib_lfid,
00279                                int index, int* id_sib );
00280     ErrorCode reorder_indices( int deg, EntityHandle* face1_conn, EntityHandle* face2_conn, int nvF,
00281                                std::vector< int >& lemap, std::vector< int >& vidx, int* leorient = NULL );
00282     ErrorCode reorder_indices( int deg, int nvF, int comb, int* childfid_map );
00283     ErrorCode reorder_indices( EntityHandle* face1_conn, EntityHandle* face2_conn, int nvF, int* conn_map, int& comb,
00284                                int* orient = NULL );
00285     ErrorCode get_lid_inci_child( EntityType type, int deg, int lfid, int leid, std::vector< int >& child_ids,
00286                                   std::vector< int >& child_lvids );
00287 
00288     // Permutation matrices
00289     struct pmat
00290     {
00291         short int num_comb;           // Number of combinations
00292         int comb[MAX_HE][MAX_HE];     // Combinations
00293         int lemap[MAX_HE][MAX_HE];    // Local edge map
00294         int orient[MAX_HE];           // Orientation
00295         int porder2[MAX_HE][MAX_HE];  // Permuted order degree 2
00296         int porder3[MAX_HE][MAX_HE];  // Permuted order degree 3
00297     };
00298 
00299     static const pmat permutation[2];
00300 
00301     // Print functions
00302     ErrorCode print_maps_1D( int level );
00303     ErrorCode print_maps_2D( int level, EntityType type );
00304     ErrorCode print_maps_3D( int level, EntityType type );
00305 
00306     // Coordinates
00307     ErrorCode compute_coordinates( int cur_level, int deg, EntityType type, EntityHandle* vbuffer, int vtotal,
00308                                    double* corner_coords, std::vector< int >& vflag, int nverts_prev );
00309 
00310     // Update the ahf maps
00311 
00312     ErrorCode update_local_ahf( int deg, EntityType type, EntityHandle* vbuffer, EntityHandle* ent_buffer, int etotal );
00313 
00314     ErrorCode update_local_ahf( int deg, EntityType type, int pat_id, EntityHandle* vbuffer, EntityHandle* ent_buffer,
00315                                 int etotal );
00316 
00317     //  ErrorCode update_global_ahf(EntityType type, int cur_level, int deg);
00318 
00319     //  ErrorCode update_global_ahf(int cur_level, int deg, std::vector<int> &pattern_ids);
00320 
00321     ErrorCode update_global_ahf( EntityType type, int cur_level, int deg, std::vector< int >* pattern_ids = NULL );
00322 
00323     ErrorCode update_global_ahf_1D( int cur_level, int deg );
00324 
00325     ErrorCode update_global_ahf_1D_sub( int cur_level, int deg );
00326 
00327     ErrorCode update_ahf_1D( int cur_level );
00328 
00329     ErrorCode update_global_ahf_2D( int cur_level, int deg );
00330 
00331     ErrorCode update_global_ahf_2D_sub( int cur_level, int deg );
00332 
00333     ErrorCode update_global_ahf_3D( int cur_level, int deg, std::vector< int >* pattern_ids = NULL );
00334 
00335     //    ErrorCode update_global_ahf_3D(int cur_level, int deg);
00336 
00337     //    ErrorCode update_global_ahf_3D(int cur_level, int deg, std::vector<int> &pattern_ids);
00338 
00339     /** Boundary extraction functions
00340      * Given a vertex at a certain level, it returns a boolean value true if it lies on the domain
00341      * boundary. Note: This is a specialization of the NestedRefine::is_entity_on_boundary function
00342      * and applies only to vertex queries. \param entity
00343      */
00344     bool is_vertex_on_boundary( const EntityHandle& entity );
00345     bool is_edge_on_boundary( const EntityHandle& entity );
00346     bool is_face_on_boundary( const EntityHandle& entity );
00347     bool is_cell_on_boundary( const EntityHandle& entity );
00348 
00349     // ErrorCode find_skin_faces(EntityHandle set, int level, int nskinF);
00350 
00351     /** Parallel communication routines
00352      * We implement two strategies to resolve the shared entities of the newly created entities.
00353      * The first strategy is to use the existing parallel merge capability which essentially uses
00354      * a coordinate-based matching of vertices and subsequently the entity handles through
00355      * their connectivities. The second strategy is an optimized and a new algorithm. It uses
00356      * the existing shared information from the coarse entities and propagates the parallel
00357      *  information appropriately.
00358      */
00359 
00360     // Send/Recv Buffer storage
00361     /*   struct pbuffer{
00362          int rank;
00363          std::vector<int> msgsize;
00364          std::vector<EntityHandle> localBuff;
00365          std::vector<EntityHandle> remoteBuff;
00366        };
00367 
00368        pbuffer commBuffers[MAX_PROCS];
00369 
00370        //Parallel tag values
00371        struct parinfo{
00372          std::multimap<EntityHandle, int> remoteps;
00373          std::multimap<EntityHandle, EntityHandle> remotehs;
00374        };
00375 
00376        parinfo parallelInfo[MAX_LEVELS];*/
00377 
00378 #ifdef MOAB_HAVE_MPI
00379 
00380     ErrorCode resolve_shared_ents_parmerge( int level, EntityHandle levelset );
00381     ErrorCode resolve_shared_ents_opt( EntityHandle* hm_set, int num_levels );
00382 
00383     ErrorCode collect_shared_entities_by_dimension( Range sharedEnts, Range& allEnts );
00384     ErrorCode collect_FList( int to_proc, Range faces, std::vector< EntityHandle >& FList,
00385                              std::vector< EntityHandle >& RList );
00386     ErrorCode collect_EList( int to_proc, Range edges, std::vector< EntityHandle >& EList,
00387                              std::vector< EntityHandle >& RList );
00388     ErrorCode collect_VList( int to_proc, Range verts, std::vector< EntityHandle >& VList,
00389                              std::vector< EntityHandle >& RList );
00390 
00391     ErrorCode decipher_remote_handles( std::vector< int >& sharedprocs, std::vector< std::vector< int > >& auxinfo,
00392                                        std::vector< std::vector< EntityHandle > >& localbuffers,
00393                                        std::vector< std::vector< EntityHandle > >& remotebuffers,
00394                                        std::multimap< EntityHandle, int >& remProcs,
00395                                        std::multimap< EntityHandle, EntityHandle >& remHandles );
00396 
00397     ErrorCode decipher_remote_handles_face( int shared_proc, int numfaces, std::vector< EntityHandle >& localFaceList,
00398                                             std::vector< EntityHandle >& remFaceList,
00399                                             std::multimap< EntityHandle, int >& remProcs,
00400                                             std::multimap< EntityHandle, EntityHandle >& remHandles );
00401 
00402     ErrorCode decipher_remote_handles_edge( int shared_proc, int numedges, std::vector< EntityHandle >& localEdgeList,
00403                                             std::vector< EntityHandle >& remEdgeList,
00404                                             std::multimap< EntityHandle, int >& remProcs,
00405                                             std::multimap< EntityHandle, EntityHandle >& remHandles );
00406 
00407     ErrorCode decipher_remote_handles_vertex( int shared_proc, int numverts,
00408                                               std::vector< EntityHandle >& localVertexList,
00409                                               std::vector< EntityHandle >& remVertexList,
00410                                               std::multimap< EntityHandle, int >& remProcs,
00411                                               std::multimap< EntityHandle, EntityHandle >& remHandles );
00412 
00413     ErrorCode update_parallel_tags( std::multimap< EntityHandle, int >& remProcs,
00414                                     std::multimap< EntityHandle, EntityHandle >& remHandles );
00415 
00416     ErrorCode get_data_from_buff( int dim, int type, int level, int entityidx, int nentities,
00417                                   std::vector< EntityHandle >& buffer, std::vector< EntityHandle >& data );
00418 
00419     bool check_for_parallelinfo( EntityHandle entity, int proc, std::multimap< EntityHandle, int >& remProcs );
00420 
00421     ErrorCode check_for_parallelinfo( EntityHandle entity, int proc,
00422                                       std::multimap< EntityHandle, EntityHandle >& remHandles,
00423                                       std::multimap< EntityHandle, int >& remProcs, EntityHandle& rhandle );
00424 
00425 #endif
00426 };
00427 
00428 }  // namespace moab
00429 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines