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