Branch data Line data Source code
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 : :
34 : : class NestedRefine
35 : : {
36 : :
37 : : public:
38 : : NestedRefine( Core* impl, ParallelComm* comm = 0, EntityHandle rset = 0 );
39 : :
40 : : ~NestedRefine();
41 : :
42 : : ErrorCode initialize();
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, int* level_degrees, std::vector< EntityHandle >& level_sets,
56 : : bool optimize = false );
57 : :
58 : : //! Given an entity and its level, return its connectivity.
59 : : /** Given an entity at a certain level, it finds the connectivity via direct access to a stored
60 : : * internal pointer to the memory to connectivity sequence for the given level. \param ent
61 : : * EntityHandle of the entity \param level Integer level of the entity for which connectivity is
62 : : * requested \param conn std::vector returning the connectivity of the entity
63 : : */
64 : :
65 : : ErrorCode get_connectivity( EntityHandle ent, int level, std::vector< EntityHandle >& conn );
66 : :
67 : : //! Given a vector of vertices and their level, return its coordinates.
68 : : /** Given a vector of vertices at a certain level, it finds the coordinates via direct access to
69 : : * a stored internal pointer to the memory to coordinate sequence for the given level. \param
70 : : * verts std::vector of the entity handles of the vertices \param num_verts The number of
71 : : * vertices \param level Integer level of the entity for which connectivity is requested \param
72 : : * coords double pointer returning the coordinates of the vertices
73 : : */
74 : :
75 : : ErrorCode get_coordinates( EntityHandle* verts, int num_verts, int level, double* coords );
76 : :
77 : : //! Get the adjacencies associated with an entity.
78 : : /** Given an entity of dimension <em>d</em>, gather all the adjacent <em>D</em> dimensional
79 : : * entities where <em>D >, = , < d </em>.
80 : : *
81 : : * \param source_entity EntityHandle to which adjacent entities have to be found.
82 : : * \param target_dimension Int Dimension of the desired adjacent entities.
83 : : * \param target_entities Vector in which the adjacent EntityHandle are returned.
84 : : */
85 : :
86 : : ErrorCode get_adjacencies( const EntityHandle source_entity, const unsigned int target_dimension,
87 : : std::vector< EntityHandle >& target_entities );
88 : :
89 : : // Interlevel parent-child or vice-versa queries
90 : : /** Given an entity from a certain level, it returns a pointer to its parent at the requested
91 : : * parent level. NOTE: This query does not support vertices.
92 : : *
93 : : * \param child EntityHandle of the entity whose parent is requested
94 : : * \param child_level Mesh level where the child exists
95 : : * \param parent_level Mesh level from which parent is requested
96 : : * \param parent Pointer to the parent in the requested parent_level
97 : : */
98 : :
99 : : ErrorCode child_to_parent( EntityHandle child, int child_level, int parent_level, EntityHandle* parent );
100 : :
101 : : /** Given an entity from a certain level, it returns a std::vector of all its children from the
102 : : * requested child level. NOTE: This query does not support vertices.
103 : : *
104 : : * \param parent EntityHandle of the entity whose children in subsequent level is requested
105 : : * \param parent_level Mesh level where the parent exists
106 : : * \param child_level Mesh level from which its children are requested
107 : : * \param children Vector containing all childrens from the requested child_level
108 : : */
109 : :
110 : : ErrorCode parent_to_child( EntityHandle parent, int parent_level, int child_level,
111 : : std::vector< EntityHandle >& children );
112 : :
113 : : /** Given a vertex from a certain level, it returns a std::vector of all entities from any
114 : : * previous levels that contains it.
115 : : *
116 : : * \param vertex EntityHandle of the vertex
117 : : * \param vert_level Mesh level of the vertex
118 : : * \param parent_level Mesh level from which entities containing vertex is requested
119 : : * \param incident_entities Vector containing entities from the parent level incident on the
120 : : * vertex
121 : : */
122 : :
123 : : ErrorCode vertex_to_entities_up( EntityHandle vertex, int vert_level, int parent_level,
124 : : std::vector< EntityHandle >& incident_entities );
125 : :
126 : : /** Given a vertex from a certain level, it returns a std::vector of all children entities of
127 : : * incident entities to vertex from any subsequent levels
128 : : *
129 : : * \param vertex EntityHandle of the vertex
130 : : * \param vert_level Mesh level of the vertex
131 : : * \param child_level Mesh level from which child entities are requested
132 : : * \param incident_entities Vector containing entities from the child level
133 : : */
134 : :
135 : : ErrorCode vertex_to_entities_down( EntityHandle vertex, int vert_level, int child_level,
136 : : std::vector< EntityHandle >& incident_entities );
137 : :
138 : : ErrorCode get_vertex_duplicates( EntityHandle vertex, int level, EntityHandle& dupvertex );
139 : :
140 : : /** Given an entity at a certain level, it returns a boolean value true if it lies on the domain
141 : : * boundary. \param entity
142 : : */
143 : :
144 : : bool is_entity_on_boundary( const EntityHandle& entity );
145 : :
146 : : ErrorCode exchange_ghosts( std::vector< EntityHandle >& lsets, int num_glayers );
147 : :
148 : : ErrorCode update_special_tags( int level, EntityHandle& lset );
149 : :
150 : : struct codeperf
151 : : {
152 : : double tm_total;
153 : : double tm_refine;
154 : : double tm_resolve;
155 : : };
156 : :
157 : : codeperf timeall;
158 : :
159 : : protected:
160 : : Core* mbImpl;
161 : : ParallelComm* pcomm;
162 : : HalfFacetRep* ahf;
163 : : CpuTimer* tm;
164 : : EntityHandle _rset;
165 : :
166 : : Range _inverts, _inedges, _infaces, _incells;
167 : :
168 : : EntityType elementype;
169 : : int meshdim, nlevels;
170 : : int level_dsequence[MAX_LEVELS];
171 : : std::map< int, int > deg_index;
172 : : bool hasghost;
173 : :
174 : : /*! \struct refPatterns
175 : : * Refinement patterns w.r.t the reference element. It consists of a locally indexed vertex list
176 : : * along with their natural coordinates, the connectivity of the subdivided entities with local
177 : : * indices, their local AHF maps along with other helper fields to aid in general book keeping
178 : : * such as avoiding vertex duplication during refinement. The entity and degree specific values
179 : : * are stored in the Templates.hpp. \sa Templates.hpp
180 : : */
181 : :
182 : : //! refPatterns
183 : : struct refPatterns
184 : : {
185 : : //! Number of new vertices on edge
186 : : short int nv_edge;
187 : : //! Number of new vertices on face, does not include those on edge
188 : : short int nv_face;
189 : : //! Number of new vertices in cell
190 : : short int nv_cell;
191 : : //! Total number of new vertices per entity
192 : : short int total_new_verts;
193 : : //! Total number of new child entities
194 : : short int total_new_ents;
195 : : //! Lower and upper indices of the new vertices
196 : : int vert_index_bnds[2];
197 : : //! Natural coordinates of the new vertices w.r.t reference
198 : : double vert_nat_coord[MAX_VERTS][3];
199 : : //! Connectivity of the new entities
200 : : int ents_conn[MAX_CHILDRENS][MAX_CONN];
201 : : //! Vertex to half-facet map of the new vertices
202 : : int v2hf[MAX_VERTS][2];
203 : : //! Opposite half-facet map of the new entities
204 : : int ents_opphfs[MAX_CHILDRENS][2 * MAX_CONN];
205 : : //! Helper: storing the local ids of vertices on each local edge
206 : : int vert_on_edges[MAX_HE][MAX_VHF];
207 : : //! Helper: storing local ids of verts on each local face, doesnt include those on edges of
208 : : //! the face
209 : : int vert_on_faces[MAX_HF][MAX_VHF];
210 : : //! Helper: stores child half-facets incident on parent half-facet. First column contain the
211 : : //! number of such children
212 : : int ents_on_pent[MAX_HF][MAX_CHILDRENS];
213 : : //! Helper: stores child ents incident on new verts on edge.
214 : : // Each triad in the column consists of :
215 : : // 1) a local child incident on the vertex on the edge
216 : : // 2) the local face id from the child
217 : : // 3) the local vertex id wrt the child connectivity
218 : : int ents_on_vedge[MAX_HE][MAX_VHF * 3];
219 : : };
220 : : //! refPatterns
221 : :
222 : : static const refPatterns refTemplates[9][MAX_DEGREE];
223 : :
224 : : //! Helper
225 : : struct intFEdge
226 : : {
227 : : //! Number of edges interior to a face
228 : : short int nie;
229 : : //! Local connectivity of the interior edges
230 : : int ieconn[12][2];
231 : : };
232 : : static const intFEdge intFacEdg[2][2];
233 : :
234 : : int get_index_from_degree( int degree );
235 : :
236 : : // HM Storage Helper
237 [ + - ][ + - ]: 1020 : struct level_memory
[ + - ][ + - ]
238 : : {
239 : : int num_verts, num_edges, num_faces, num_cells;
240 : : EntityHandle start_vertex, start_edge, start_face, start_cell;
241 : : std::vector< double* > coordinates;
242 : : EntityHandle *edge_conn, *face_conn, *cell_conn;
243 : : Range verts, edges, faces, cells;
244 : : };
245 : :
246 : : level_memory level_mesh[MAX_LEVELS];
247 : :
248 : : // Basic Functions
249 : :
250 : : // Estimate and create storage for the levels
251 : : ErrorCode estimate_hm_storage( EntityHandle set, int level_degree, int cur_level, int hmest[4] );
252 : : ErrorCode create_hm_storage_single_level( EntityHandle* set, int cur_level, int estL[4] );
253 : :
254 : : // Generate HM : Construct the hierarchical mesh: 1D, 2D, 3D
255 : : ErrorCode generate_hm( int* level_degrees, int num_level, EntityHandle* hm_set, bool optimize );
256 : : ErrorCode construct_hm_entities( int cur_level, int deg );
257 : : ErrorCode construct_hm_1D( int cur_level, int deg );
258 : : ErrorCode construct_hm_1D( int cur_level, int deg, EntityType type, std::vector< EntityHandle >& trackverts );
259 : : ErrorCode construct_hm_2D( int cur_level, int deg );
260 : : ErrorCode construct_hm_2D( int cur_level, int deg, EntityType type, std::vector< EntityHandle >& trackvertsC_edg,
261 : : std::vector< EntityHandle >& trackvertsF );
262 : : ErrorCode construct_hm_3D( int cur_level, int deg );
263 : :
264 : : ErrorCode subdivide_cells( EntityType type, int cur_level, int deg );
265 : : ErrorCode subdivide_tets( int cur_level, int deg );
266 : :
267 : : // General helper functions
268 : : ErrorCode copy_vertices_from_prev_level( int cur_level );
269 : : ErrorCode count_subentities( EntityHandle set, int cur_level, int* nedges, int* nfaces );
270 : : ErrorCode get_octahedron_corner_coords( int cur_level, int deg, EntityHandle* vbuffer, double* ocoords );
271 : : int find_shortest_diagonal_octahedron( int cur_level, int deg, EntityHandle* vbuffer );
272 : : int get_local_vid( EntityHandle vid, EntityHandle ent, int level );
273 : :
274 : : // Book-keeping functions
275 : : ErrorCode update_tracking_verts( EntityHandle cid, int cur_level, int deg,
276 : : std::vector< EntityHandle >& trackvertsC_edg,
277 : : std::vector< EntityHandle >& trackvertsC_face, EntityHandle* vbuffer );
278 : : ErrorCode reorder_indices( int cur_level, int deg, EntityHandle cell, int lfid, EntityHandle sib_cell, int sib_lfid,
279 : : int index, int* id_sib );
280 : : ErrorCode reorder_indices( int deg, EntityHandle* face1_conn, EntityHandle* face2_conn, int nvF,
281 : : std::vector< int >& lemap, std::vector< int >& vidx, int* leorient = NULL );
282 : : ErrorCode reorder_indices( int deg, int nvF, int comb, int* childfid_map );
283 : : ErrorCode reorder_indices( EntityHandle* face1_conn, EntityHandle* face2_conn, int nvF, int* conn_map, int& comb,
284 : : int* orient = NULL );
285 : : ErrorCode get_lid_inci_child( EntityType type, int deg, int lfid, int leid, std::vector< int >& child_ids,
286 : : std::vector< int >& child_lvids );
287 : :
288 : : // Permutation matrices
289 : : struct pmat
290 : : {
291 : : short int num_comb; // Number of combinations
292 : : int comb[MAX_HE][MAX_HE]; // Combinations
293 : : int lemap[MAX_HE][MAX_HE]; // Local edge map
294 : : int orient[MAX_HE]; // Orientation
295 : : int porder2[MAX_HE][MAX_HE]; // Permuted order degree 2
296 : : int porder3[MAX_HE][MAX_HE]; // Permuted order degree 3
297 : : };
298 : :
299 : : static const pmat permutation[2];
300 : :
301 : : // Print functions
302 : : ErrorCode print_maps_1D( int level );
303 : : ErrorCode print_maps_2D( int level, EntityType type );
304 : : ErrorCode print_maps_3D( int level, EntityType type );
305 : :
306 : : // Coordinates
307 : : ErrorCode compute_coordinates( int cur_level, int deg, EntityType type, EntityHandle* vbuffer, int vtotal,
308 : : double* corner_coords, std::vector< int >& vflag, int nverts_prev );
309 : :
310 : : // Update the ahf maps
311 : :
312 : : ErrorCode update_local_ahf( int deg, EntityType type, EntityHandle* vbuffer, EntityHandle* ent_buffer, int etotal );
313 : :
314 : : ErrorCode update_local_ahf( int deg, EntityType type, int pat_id, EntityHandle* vbuffer, EntityHandle* ent_buffer,
315 : : int etotal );
316 : :
317 : : // ErrorCode update_global_ahf(EntityType type, int cur_level, int deg);
318 : :
319 : : // ErrorCode update_global_ahf(int cur_level, int deg, std::vector<int> &pattern_ids);
320 : :
321 : : ErrorCode update_global_ahf( EntityType type, int cur_level, int deg, std::vector< int >* pattern_ids = NULL );
322 : :
323 : : ErrorCode update_global_ahf_1D( int cur_level, int deg );
324 : :
325 : : ErrorCode update_global_ahf_1D_sub( int cur_level, int deg );
326 : :
327 : : ErrorCode update_ahf_1D( int cur_level );
328 : :
329 : : ErrorCode update_global_ahf_2D( int cur_level, int deg );
330 : :
331 : : ErrorCode update_global_ahf_2D_sub( int cur_level, int deg );
332 : :
333 : : ErrorCode update_global_ahf_3D( int cur_level, int deg, std::vector< int >* pattern_ids = NULL );
334 : :
335 : : // ErrorCode update_global_ahf_3D(int cur_level, int deg);
336 : :
337 : : // ErrorCode update_global_ahf_3D(int cur_level, int deg, std::vector<int> &pattern_ids);
338 : :
339 : : /** Boundary extraction functions
340 : : * Given a vertex at a certain level, it returns a boolean value true if it lies on the domain
341 : : * boundary. Note: This is a specialization of the NestedRefine::is_entity_on_boundary function
342 : : * and applies only to vertex queries. \param entity
343 : : */
344 : : bool is_vertex_on_boundary( const EntityHandle& entity );
345 : : bool is_edge_on_boundary( const EntityHandle& entity );
346 : : bool is_face_on_boundary( const EntityHandle& entity );
347 : : bool is_cell_on_boundary( const EntityHandle& entity );
348 : :
349 : : // ErrorCode find_skin_faces(EntityHandle set, int level, int nskinF);
350 : :
351 : : /** Parallel communication routines
352 : : * We implement two strategies to resolve the shared entities of the newly created entities.
353 : : * The first strategy is to use the existing parallel merge capability which essentially uses
354 : : * a coordinate-based matching of vertices and subsequently the entity handles through
355 : : * their connectivities. The second strategy is an optimized and a new algorithm. It uses
356 : : * the existing shared information from the coarse entities and propagates the parallel
357 : : * information appropriately.
358 : : */
359 : :
360 : : // Send/Recv Buffer storage
361 : : /* struct pbuffer{
362 : : int rank;
363 : : std::vector<int> msgsize;
364 : : std::vector<EntityHandle> localBuff;
365 : : std::vector<EntityHandle> remoteBuff;
366 : : };
367 : :
368 : : pbuffer commBuffers[MAX_PROCS];
369 : :
370 : : //Parallel tag values
371 : : struct parinfo{
372 : : std::multimap<EntityHandle, int> remoteps;
373 : : std::multimap<EntityHandle, EntityHandle> remotehs;
374 : : };
375 : :
376 : : parinfo parallelInfo[MAX_LEVELS];*/
377 : :
378 : : #ifdef MOAB_HAVE_MPI
379 : :
380 : : ErrorCode resolve_shared_ents_parmerge( int level, EntityHandle levelset );
381 : : ErrorCode resolve_shared_ents_opt( EntityHandle* hm_set, int num_levels );
382 : :
383 : : ErrorCode collect_shared_entities_by_dimension( Range sharedEnts, Range& allEnts );
384 : : ErrorCode collect_FList( int to_proc, Range faces, std::vector< EntityHandle >& FList,
385 : : std::vector< EntityHandle >& RList );
386 : : ErrorCode collect_EList( int to_proc, Range edges, std::vector< EntityHandle >& EList,
387 : : std::vector< EntityHandle >& RList );
388 : : ErrorCode collect_VList( int to_proc, Range verts, std::vector< EntityHandle >& VList,
389 : : std::vector< EntityHandle >& RList );
390 : :
391 : : ErrorCode decipher_remote_handles( std::vector< int >& sharedprocs, std::vector< std::vector< int > >& auxinfo,
392 : : std::vector< std::vector< EntityHandle > >& localbuffers,
393 : : std::vector< std::vector< EntityHandle > >& remotebuffers,
394 : : std::multimap< EntityHandle, int >& remProcs,
395 : : std::multimap< EntityHandle, EntityHandle >& remHandles );
396 : :
397 : : ErrorCode decipher_remote_handles_face( int shared_proc, int numfaces, std::vector< EntityHandle >& localFaceList,
398 : : std::vector< EntityHandle >& remFaceList,
399 : : std::multimap< EntityHandle, int >& remProcs,
400 : : std::multimap< EntityHandle, EntityHandle >& remHandles );
401 : :
402 : : ErrorCode decipher_remote_handles_edge( int shared_proc, int numedges, std::vector< EntityHandle >& localEdgeList,
403 : : std::vector< EntityHandle >& remEdgeList,
404 : : std::multimap< EntityHandle, int >& remProcs,
405 : : std::multimap< EntityHandle, EntityHandle >& remHandles );
406 : :
407 : : ErrorCode decipher_remote_handles_vertex( int shared_proc, int numverts,
408 : : std::vector< EntityHandle >& localVertexList,
409 : : std::vector< EntityHandle >& remVertexList,
410 : : std::multimap< EntityHandle, int >& remProcs,
411 : : std::multimap< EntityHandle, EntityHandle >& remHandles );
412 : :
413 : : ErrorCode update_parallel_tags( std::multimap< EntityHandle, int >& remProcs,
414 : : std::multimap< EntityHandle, EntityHandle >& remHandles );
415 : :
416 : : ErrorCode get_data_from_buff( int dim, int type, int level, int entityidx, int nentities,
417 : : std::vector< EntityHandle >& buffer, std::vector< EntityHandle >& data );
418 : :
419 : : bool check_for_parallelinfo( EntityHandle entity, int proc, std::multimap< EntityHandle, int >& remProcs );
420 : :
421 : : ErrorCode check_for_parallelinfo( EntityHandle entity, int proc,
422 : : std::multimap< EntityHandle, EntityHandle >& remHandles,
423 : : std::multimap< EntityHandle, int >& remProcs, EntityHandle& rhandle );
424 : :
425 : : #endif
426 : : };
427 : :
428 : : } // namespace moab
429 : : #endif
|
