1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
/*! \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