Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
NestedRefine.hpp
Go to the documentation of this file.
1 /*! \file NestedRefine.hpp
2  * This class implements the generation of a hierarchy of meshes via uniform refinement from an
3  * input mesh. The internal upper bound on the number of levels is set to 20.
4  */
5 
6 #ifndef NESTED_REFINE_HPP
7 #define NESTED_REFINE_HPP
8 
9 #include "moab/MOABConfig.h"
10 #include "moab/Range.hpp"
11 #include "moab/CN.hpp"
12 #include <map>
13 #include <set>
14 
15 namespace moab
16 {
17 
18 #define MAX_DEGREE 3
19 #define MAX_VERTS 64
20 #define MAX_CHILDRENS 27
21 #define MAX_HE 12
22 #define MAX_HF 6
23 #define MAX_CONN 8
24 #define MAX_VHF 20
25 #define MAX_LEVELS 20
26 #define MAX_PROCS 64
27 
28 // Forward Declarations
29 class Core;
30 class HalfFacetRep;
31 class ParallelComm;
32 class CpuTimer;
33 
35 {
36 
37  public:
38  NestedRefine( Core* impl, ParallelComm* comm = 0, EntityHandle rset = 0 );
39 
40  ~NestedRefine();
41 
43 
44  //! \brief Generate a mesh hierarchy.
45  /** Given a mesh, generate a sequence of meshes via uniform refinement. The inputs are: a) an
46  * array(level_degrees) storing the degrees which will be used to refine the previous level mesh
47  * to generate a new level and b) the total number of levels(should be same length as that of
48  * the array in a). Each mesh level in the hierarchy are stored in different meshsets whose
49  * handles are returned after the hierarchy generation. These handles can be used to work with a
50  * specific mesh level. \param level_degrees Integer array storing the degrees used in each
51  * level. \param num_level The total number of levels in the hierarchy. \param hm_set
52  * EntityHandle STL vector that returns the handles of the sets created for each mesh level.
53  */
54 
55  ErrorCode generate_mesh_hierarchy( int num_level,
56  int* level_degrees,
57  std::vector< EntityHandle >& level_sets,
58  bool optimize = false );
59 
60  //! Given an entity and its level, return its connectivity.
61  /** Given an entity at a certain level, it finds the connectivity via direct access to a stored
62  * internal pointer to the memory to connectivity sequence for the given level. \param ent
63  * EntityHandle of the entity \param level Integer level of the entity for which connectivity is
64  * requested \param conn std::vector returning the connectivity of the entity
65  */
66 
67  ErrorCode get_connectivity( EntityHandle ent, int level, std::vector< EntityHandle >& conn );
68 
69  //! Given a vector of vertices and their level, return its coordinates.
70  /** Given a vector of vertices at a certain level, it finds the coordinates via direct access to
71  * a stored internal pointer to the memory to coordinate sequence for the given level. \param
72  * verts std::vector of the entity handles of the vertices \param num_verts The number of
73  * vertices \param level Integer level of the entity for which connectivity is requested \param
74  * coords double pointer returning the coordinates of the vertices
75  */
76 
77  ErrorCode get_coordinates( EntityHandle* verts, int num_verts, int level, double* coords );
78 
79  //! Get the adjacencies associated with an entity.
80  /** Given an entity of dimension <em>d</em>, gather all the adjacent <em>D</em> dimensional
81  * entities where <em>D >, = , < d </em>.
82  *
83  * \param source_entity EntityHandle to which adjacent entities have to be found.
84  * \param target_dimension Int Dimension of the desired adjacent entities.
85  * \param target_entities Vector in which the adjacent EntityHandle are returned.
86  */
87 
88  ErrorCode get_adjacencies( const EntityHandle source_entity,
89  const unsigned int target_dimension,
90  std::vector< EntityHandle >& target_entities );
91 
92  // Interlevel parent-child or vice-versa queries
93  /** Given an entity from a certain level, it returns a pointer to its parent at the requested
94  * parent level. NOTE: This query does not support vertices.
95  *
96  * \param child EntityHandle of the entity whose parent is requested
97  * \param child_level Mesh level where the child exists
98  * \param parent_level Mesh level from which parent is requested
99  * \param parent Pointer to the parent in the requested parent_level
100  */
101 
102  ErrorCode child_to_parent( EntityHandle child, int child_level, int parent_level, EntityHandle* parent );
103 
104  /** Given an entity from a certain level, it returns a std::vector of all its children from the
105  * requested child level. NOTE: This query does not support vertices.
106  *
107  * \param parent EntityHandle of the entity whose children in subsequent level is requested
108  * \param parent_level Mesh level where the parent exists
109  * \param child_level Mesh level from which its children are requested
110  * \param children Vector containing all childrens from the requested child_level
111  */
112 
114  int parent_level,
115  int child_level,
116  std::vector< EntityHandle >& children );
117 
118  /** Given a vertex from a certain level, it returns a std::vector of all entities from any
119  * previous levels that contains it.
120  *
121  * \param vertex EntityHandle of the vertex
122  * \param vert_level Mesh level of the vertex
123  * \param parent_level Mesh level from which entities containing vertex is requested
124  * \param incident_entities Vector containing entities from the parent level incident on the
125  * vertex
126  */
127 
129  int vert_level,
130  int parent_level,
131  std::vector< EntityHandle >& incident_entities );
132 
133  /** Given a vertex from a certain level, it returns a std::vector of all children entities of
134  * incident entities to vertex from any subsequent levels
135  *
136  * \param vertex EntityHandle of the vertex
137  * \param vert_level Mesh level of the vertex
138  * \param child_level Mesh level from which child entities are requested
139  * \param incident_entities Vector containing entities from the child level
140  */
141 
143  int vert_level,
144  int child_level,
145  std::vector< EntityHandle >& incident_entities );
146 
148 
149  /** Given an entity at a certain level, it returns a boolean value true if it lies on the domain
150  * boundary. \param entity
151  */
152 
153  bool is_entity_on_boundary( const EntityHandle& entity );
154 
155  ErrorCode exchange_ghosts( std::vector< EntityHandle >& lsets, int num_glayers );
156 
157  ErrorCode update_special_tags( int level, EntityHandle& lset );
158 
159  struct codeperf
160  {
161  double tm_total;
162  double tm_refine;
163  double tm_resolve;
164  };
165 
167 
168  protected:
174 
176 
177  EntityType elementype;
180  std::map< int, int > deg_index;
181  bool hasghost;
182 
183  /*! \struct refPatterns
184  * Refinement patterns w.r.t the reference element. It consists of a locally indexed vertex list
185  * along with their natural coordinates, the connectivity of the subdivided entities with local
186  * indices, their local AHF maps along with other helper fields to aid in general book keeping
187  * such as avoiding vertex duplication during refinement. The entity and degree specific values
188  * are stored in the Templates.hpp. \sa Templates.hpp
189  */
190 
191  //! refPatterns
192  struct refPatterns
193  {
194  //! Number of new vertices on edge
195  short int nv_edge;
196  //! Number of new vertices on face, does not include those on edge
197  short int nv_face;
198  //! Number of new vertices in cell
199  short int nv_cell;
200  //! Total number of new vertices per entity
201  short int total_new_verts;
202  //! Total number of new child entities
203  short int total_new_ents;
204  //! Lower and upper indices of the new vertices
206  //! Natural coordinates of the new vertices w.r.t reference
208  //! Connectivity of the new entities
210  //! Vertex to half-facet map of the new vertices
211  int v2hf[MAX_VERTS][2];
212  //! Opposite half-facet map of the new entities
214  //! Helper: storing the local ids of vertices on each local edge
216  //! Helper: storing local ids of verts on each local face, doesnt include those on edges of
217  //! the face
219  //! Helper: stores child half-facets incident on parent half-facet. First column contain the
220  //! number of such children
222  //! Helper: stores child ents incident on new verts on edge.
223  // Each triad in the column consists of :
224  // 1) a local child incident on the vertex on the edge
225  // 2) the local face id from the child
226  // 3) the local vertex id wrt the child connectivity
228  };
229  //! refPatterns
230 
232 
233  //! Helper
234  struct intFEdge
235  {
236  //! Number of edges interior to a face
237  short int nie;
238  //! Local connectivity of the interior edges
239  int ieconn[12][2];
240  };
241  static const intFEdge intFacEdg[2][2];
242 
243  int get_index_from_degree( int degree );
244 
245  // HM Storage Helper
247  {
250  std::vector< double* > coordinates;
253  };
254 
256 
257  // Basic Functions
258 
259  // Estimate and create storage for the levels
260  ErrorCode estimate_hm_storage( EntityHandle set, int level_degree, int cur_level, int hmest[4] );
261  ErrorCode create_hm_storage_single_level( EntityHandle* set, int cur_level, int estL[4] );
262 
263  // Generate HM : Construct the hierarchical mesh: 1D, 2D, 3D
264  ErrorCode generate_hm( int* level_degrees, int num_level, EntityHandle* hm_set, bool optimize );
265  ErrorCode construct_hm_entities( int cur_level, int deg );
266  ErrorCode construct_hm_1D( int cur_level, int deg );
267  ErrorCode construct_hm_1D( int cur_level, int deg, EntityType type, std::vector< EntityHandle >& trackverts );
268  ErrorCode construct_hm_2D( int cur_level, int deg );
269  ErrorCode construct_hm_2D( int cur_level,
270  int deg,
271  EntityType type,
272  std::vector< EntityHandle >& trackvertsC_edg,
273  std::vector< EntityHandle >& trackvertsF );
274  ErrorCode construct_hm_3D( int cur_level, int deg );
275 
276  ErrorCode subdivide_cells( EntityType type, int cur_level, int deg );
277  ErrorCode subdivide_tets( int cur_level, int deg );
278 
279  // General helper functions
280  ErrorCode copy_vertices_from_prev_level( int cur_level );
281  ErrorCode count_subentities( EntityHandle set, int cur_level, int* nedges, int* nfaces );
282  ErrorCode get_octahedron_corner_coords( int cur_level, int deg, EntityHandle* vbuffer, double* ocoords );
283  int find_shortest_diagonal_octahedron( int cur_level, int deg, EntityHandle* vbuffer );
284  int get_local_vid( EntityHandle vid, EntityHandle ent, int level );
285 
286  // Book-keeping functions
288  int cur_level,
289  int deg,
290  std::vector< EntityHandle >& trackvertsC_edg,
291  std::vector< EntityHandle >& trackvertsC_face,
292  EntityHandle* vbuffer );
293  ErrorCode reorder_indices( int cur_level,
294  int deg,
295  EntityHandle cell,
296  int lfid,
297  EntityHandle sib_cell,
298  int sib_lfid,
299  int index,
300  int* id_sib );
301  ErrorCode reorder_indices( int deg,
302  EntityHandle* face1_conn,
303  EntityHandle* face2_conn,
304  int nvF,
305  std::vector< int >& lemap,
306  std::vector< int >& vidx,
307  int* leorient = NULL );
308  ErrorCode reorder_indices( int deg, int nvF, int comb, int* childfid_map );
310  EntityHandle* face2_conn,
311  int nvF,
312  int* conn_map,
313  int& comb,
314  int* orient = NULL );
315  ErrorCode get_lid_inci_child( EntityType type,
316  int deg,
317  int lfid,
318  int leid,
319  std::vector< int >& child_ids,
320  std::vector< int >& child_lvids );
321 
322  // Permutation matrices
323  struct pmat
324  {
325  short int num_comb; // Number of combinations
326  int comb[MAX_HE][MAX_HE]; // Combinations
327  int lemap[MAX_HE][MAX_HE]; // Local edge map
328  int orient[MAX_HE]; // Orientation
329  int porder2[MAX_HE][MAX_HE]; // Permuted order degree 2
330  int porder3[MAX_HE][MAX_HE]; // Permuted order degree 3
331  };
332 
333  static const pmat permutation[2];
334 
335  // Print functions
336  ErrorCode print_maps_1D( int level );
337  ErrorCode print_maps_2D( int level, EntityType type );
338  ErrorCode print_maps_3D( int level, EntityType type );
339 
340  // Coordinates
341  ErrorCode compute_coordinates( int cur_level,
342  int deg,
343  EntityType type,
344  EntityHandle* vbuffer,
345  int vtotal,
346  double* corner_coords,
347  std::vector< int >& vflag,
348  int nverts_prev );
349 
350  // Update the ahf maps
351 
352  ErrorCode update_local_ahf( int deg, EntityType type, EntityHandle* vbuffer, EntityHandle* ent_buffer, int etotal );
353 
354  ErrorCode update_local_ahf( int deg,
355  EntityType type,
356  int pat_id,
357  EntityHandle* vbuffer,
358  EntityHandle* ent_buffer,
359  int etotal );
360 
361  // ErrorCode update_global_ahf(EntityType type, int cur_level, int deg);
362 
363  // ErrorCode update_global_ahf(int cur_level, int deg, std::vector<int> &pattern_ids);
364 
365  ErrorCode update_global_ahf( EntityType type, int cur_level, int deg, std::vector< int >* pattern_ids = NULL );
366 
367  ErrorCode update_global_ahf_1D( int cur_level, int deg );
368 
369  ErrorCode update_global_ahf_1D_sub( int cur_level, int deg );
370 
371  ErrorCode update_ahf_1D( int cur_level );
372 
373  ErrorCode update_global_ahf_2D( int cur_level, int deg );
374 
375  ErrorCode update_global_ahf_2D_sub( int cur_level, int deg );
376 
377  ErrorCode update_global_ahf_3D( int cur_level, int deg, std::vector< int >* pattern_ids = NULL );
378 
379  // ErrorCode update_global_ahf_3D(int cur_level, int deg);
380 
381  // ErrorCode update_global_ahf_3D(int cur_level, int deg, std::vector<int> &pattern_ids);
382 
383  /** Boundary extraction functions
384  * Given a vertex at a certain level, it returns a boolean value true if it lies on the domain
385  * boundary. Note: This is a specialization of the NestedRefine::is_entity_on_boundary function
386  * and applies only to vertex queries. \param entity
387  */
388  bool is_vertex_on_boundary( const EntityHandle& entity );
389  bool is_edge_on_boundary( const EntityHandle& entity );
390  bool is_face_on_boundary( const EntityHandle& entity );
391  bool is_cell_on_boundary( const EntityHandle& entity );
392 
393  // ErrorCode find_skin_faces(EntityHandle set, int level, int nskinF);
394 
395  /** Parallel communication routines
396  * We implement two strategies to resolve the shared entities of the newly created entities.
397  * The first strategy is to use the existing parallel merge capability which essentially uses
398  * a coordinate-based matching of vertices and subsequently the entity handles through
399  * their connectivities. The second strategy is an optimized and a new algorithm. It uses
400  * the existing shared information from the coarse entities and propagates the parallel
401  * information appropriately.
402  */
403 
404  // Send/Recv Buffer storage
405  /* struct pbuffer{
406  int rank;
407  std::vector<int> msgsize;
408  std::vector<EntityHandle> localBuff;
409  std::vector<EntityHandle> remoteBuff;
410  };
411 
412  pbuffer commBuffers[MAX_PROCS];
413 
414  //Parallel tag values
415  struct parinfo{
416  std::multimap<EntityHandle, int> remoteps;
417  std::multimap<EntityHandle, EntityHandle> remotehs;
418  };
419 
420  parinfo parallelInfo[MAX_LEVELS];*/
421 
422 #ifdef MOAB_HAVE_MPI
423 
424  ErrorCode resolve_shared_ents_parmerge( int level, EntityHandle levelset );
425  ErrorCode resolve_shared_ents_opt( EntityHandle* hm_set, int num_levels );
426 
427  ErrorCode collect_shared_entities_by_dimension( Range sharedEnts, Range& allEnts );
428  ErrorCode collect_FList( int to_proc,
429  Range faces,
430  std::vector< EntityHandle >& FList,
431  std::vector< EntityHandle >& RList );
432  ErrorCode collect_EList( int to_proc,
433  Range edges,
434  std::vector< EntityHandle >& EList,
435  std::vector< EntityHandle >& RList );
436  ErrorCode collect_VList( int to_proc,
437  Range verts,
438  std::vector< EntityHandle >& VList,
439  std::vector< EntityHandle >& RList );
440 
441  ErrorCode decipher_remote_handles( std::vector< int >& sharedprocs,
442  std::vector< std::vector< int > >& auxinfo,
443  std::vector< std::vector< EntityHandle > >& localbuffers,
444  std::vector< std::vector< EntityHandle > >& remotebuffers,
445  std::multimap< EntityHandle, int >& remProcs,
446  std::multimap< EntityHandle, EntityHandle >& remHandles );
447 
448  ErrorCode decipher_remote_handles_face( int shared_proc,
449  int numfaces,
450  std::vector< EntityHandle >& localFaceList,
451  std::vector< EntityHandle >& remFaceList,
452  std::multimap< EntityHandle, int >& remProcs,
453  std::multimap< EntityHandle, EntityHandle >& remHandles );
454 
455  ErrorCode decipher_remote_handles_edge( int shared_proc,
456  int numedges,
457  std::vector< EntityHandle >& localEdgeList,
458  std::vector< EntityHandle >& remEdgeList,
459  std::multimap< EntityHandle, int >& remProcs,
460  std::multimap< EntityHandle, EntityHandle >& remHandles );
461 
462  ErrorCode decipher_remote_handles_vertex( int shared_proc,
463  int numverts,
464  std::vector< EntityHandle >& localVertexList,
465  std::vector< EntityHandle >& remVertexList,
466  std::multimap< EntityHandle, int >& remProcs,
467  std::multimap< EntityHandle, EntityHandle >& remHandles );
468 
469  ErrorCode update_parallel_tags( std::multimap< EntityHandle, int >& remProcs,
470  std::multimap< EntityHandle, EntityHandle >& remHandles );
471 
472  ErrorCode get_data_from_buff( int dim,
473  int type,
474  int level,
475  int entityidx,
476  int nentities,
477  std::vector< EntityHandle >& buffer,
478  std::vector< EntityHandle >& data );
479 
480  bool check_for_parallelinfo( EntityHandle entity, int proc, std::multimap< EntityHandle, int >& remProcs );
481 
482  ErrorCode check_for_parallelinfo( EntityHandle entity,
483  int proc,
484  std::multimap< EntityHandle, EntityHandle >& remHandles,
485  std::multimap< EntityHandle, int >& remProcs,
486  EntityHandle& rhandle );
487 
488 #endif
489 };
490 
491 } // namespace moab
492 #endif