Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
moab::NestedRefine Class Reference

#include <NestedRefine.hpp>

+ Collaboration diagram for moab::NestedRefine:

Classes

struct  codeperf
 
struct  intFEdge
 Helper. More...
 
struct  level_memory
 
struct  pmat
 
struct  refPatterns
 refPatterns More...
 

Public Member Functions

 NestedRefine (Core *impl, ParallelComm *comm=0, EntityHandle rset=0)
 
 ~NestedRefine ()
 
ErrorCode initialize ()
 
ErrorCode generate_mesh_hierarchy (int num_level, int *level_degrees, std::vector< EntityHandle > &level_sets, bool optimize=false)
 Generate a mesh hierarchy. More...
 
ErrorCode get_connectivity (EntityHandle ent, int level, std::vector< EntityHandle > &conn)
 Given an entity and its level, return its connectivity. More...
 
ErrorCode get_coordinates (EntityHandle *verts, int num_verts, int level, double *coords)
 Given a vector of vertices and their level, return its coordinates. More...
 
ErrorCode get_adjacencies (const EntityHandle source_entity, const unsigned int target_dimension, std::vector< EntityHandle > &target_entities)
 Get the adjacencies associated with an entity. More...
 
ErrorCode child_to_parent (EntityHandle child, int child_level, int parent_level, EntityHandle *parent)
 
ErrorCode parent_to_child (EntityHandle parent, int parent_level, int child_level, std::vector< EntityHandle > &children)
 
ErrorCode vertex_to_entities_up (EntityHandle vertex, int vert_level, int parent_level, std::vector< EntityHandle > &incident_entities)
 
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)
 
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)
 

Public Attributes

codeperf timeall
 

Protected Member Functions

int get_index_from_degree (int degree)
 
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])
 
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)
 
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)
 
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)
 
ErrorCode print_maps_1D (int level)
 
ErrorCode print_maps_2D (int level, EntityType type)
 
ErrorCode print_maps_3D (int level, EntityType type)
 
ErrorCode compute_coordinates (int cur_level, int deg, EntityType type, EntityHandle *vbuffer, int vtotal, double *corner_coords, std::vector< int > &vflag, int nverts_prev)
 
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, 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)
 
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)
 

Protected Attributes

CorembImpl
 
ParallelCommpcomm
 
HalfFacetRepahf
 
CpuTimertm
 
EntityHandle _rset
 
Range _inverts
 
Range _inedges
 
Range _infaces
 
Range _incells
 
EntityType elementype
 
int meshdim
 
int nlevels
 
int level_dsequence [MAX_LEVELS]
 
std::map< int, int > deg_index
 
bool hasghost
 
level_memory level_mesh [MAX_LEVELS]
 

Static Protected Attributes

static const refPatterns refTemplates [9][MAX_DEGREE]
 refPatterns More...
 
static const intFEdge intFacEdg [2][2]
 
static const pmat permutation [2]
 

Detailed Description

Definition at line 34 of file NestedRefine.hpp.

Constructor & Destructor Documentation

◆ NestedRefine()

moab::NestedRefine::NestedRefine ( Core impl,
ParallelComm comm = 0,
EntityHandle  rset = 0 
)

Definition at line 24 of file NestedRefine.cpp.

25  : mbImpl( impl ), pcomm( comm ), _rset( rset )
26 {
28  assert( NULL != impl );
29 
30 #ifdef MOAB_HAVE_MPI
31  // Get the Parallel Comm instance to prepare all new sets to work in parallel
32  // in case the user did not provide any arguments
33  if( !comm ) pcomm = moab::ParallelComm::get_pcomm( mbImpl, 0 );
34 #endif
35  error = initialize();
36  if( error != MB_SUCCESS )
37  {
38  std::cout << "Error initializing NestedRefine\n" << std::endl;
39  exit( 1 );
40  }
41 }

References moab::error(), ErrorCode, moab::ParallelComm::get_pcomm(), initialize(), MB_SUCCESS, mbImpl, and pcomm.

◆ ~NestedRefine()

moab::NestedRefine::~NestedRefine ( )

Definition at line 43 of file NestedRefine.cpp.

44 {
45 #ifdef MOAB_HAVE_AHF
46  ahf = NULL;
47 #else
48  delete ahf;
49 #endif
50  delete tm;
51 }

References ahf, and tm.

Member Function Documentation

◆ child_to_parent()

ErrorCode moab::NestedRefine::child_to_parent ( EntityHandle  child,
int  child_level,
int  parent_level,
EntityHandle parent 
)

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.

Parameters
childEntityHandle of the entity whose parent is requested
child_levelMesh level where the child exists
parent_levelMesh level from which parent is requested
parentPointer to the parent in the requested parent_level

Definition at line 234 of file NestedRefine.cpp.

235 {
236  assert( ( child_level > 0 ) && ( child_level > parent_level ) );
237  EntityType type = mbImpl->type_from_handle( child );
238  assert( type != MBVERTEX );
239 
240  int child_index;
241  if( type == MBEDGE )
242  child_index = child - level_mesh[child_level - 1].start_edge;
243  else if( type == MBTRI || type == MBQUAD )
244  child_index = child - level_mesh[child_level - 1].start_face;
245  else if( type == MBTET || type == MBHEX )
246  child_index = child - level_mesh[child_level - 1].start_cell;
247  else
248  MB_SET_ERR( MB_FAILURE, "Requesting parent for unsupported entity type" );
249 
250  int parent_index;
251  int l = child_level - parent_level;
252  for( int i = 0; i < l; i++ )
253  {
254  int d = get_index_from_degree( level_dsequence[child_level - i - 1] );
255  int nch = refTemplates[type - 1][d].total_new_ents;
256  child_index = child_index / nch;
257  }
258  parent_index = child_index;
259 
260  if( type == MBEDGE )
261  {
262  if( parent_level > 0 )
263  *parent = level_mesh[parent_level - 1].start_edge + parent_index;
264  else
265  *parent = _inedges[parent_index];
266  }
267  else if( type == MBTRI || type == MBQUAD )
268  {
269  if( parent_level > 0 )
270  *parent = level_mesh[parent_level - 1].start_face + parent_index;
271  else
272  *parent = _infaces[parent_index];
273  }
274  else if( type == MBTET || type == MBHEX )
275  {
276  if( parent_level > 0 )
277  *parent = level_mesh[parent_level - 1].start_cell + parent_index;
278  else
279  *parent = _incells[parent_index];
280  }
281 
282  return MB_SUCCESS;
283 }

References _incells, _inedges, _infaces, child, get_index_from_degree(), level_dsequence, level_mesh, MB_SET_ERR, MB_SUCCESS, MBEDGE, MBHEX, mbImpl, MBQUAD, MBTET, MBTRI, MBVERTEX, refTemplates, moab::NestedRefine::level_memory::start_cell, moab::NestedRefine::level_memory::start_edge, moab::NestedRefine::level_memory::start_face, moab::NestedRefine::refPatterns::total_new_ents, and moab::Core::type_from_handle().

Referenced by vertex_to_entities_up().

◆ compute_coordinates()

ErrorCode moab::NestedRefine::compute_coordinates ( int  cur_level,
int  deg,
EntityType  type,
EntityHandle vbuffer,
int  vtotal,
double *  corner_coords,
std::vector< int > &  vflag,
int  nverts_prev 
)
protected

Definition at line 1730 of file NestedRefine.cpp.

1738 {
1739  EntityHandle vstart = level_mesh[cur_level].start_vertex;
1740  int d = get_index_from_degree( deg );
1741 
1742  if( type == MBTRI )
1743  {
1744  double xi, eta, N[3];
1745  int findex = mbImpl->type_from_handle( *( _infaces.begin() ) ) - 1;
1746 
1747  for( int i = 3; i < vtotal; i++ )
1748  {
1749  if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
1750 
1751  xi = refTemplates[findex][d].vert_nat_coord[i - 3][0];
1752  eta = refTemplates[findex][d].vert_nat_coord[i - 3][1];
1753  N[0] = 1 - xi - eta;
1754  N[1] = xi;
1755  N[2] = eta;
1756 
1757  double x = 0, y = 0, z = 0;
1758  for( int j = 0; j < 3; j++ )
1759  {
1760  x += N[j] * corner_coords[3 * j];
1761  y += N[j] * corner_coords[3 * j + 1];
1762  z += N[j] * corner_coords[3 * j + 2];
1763  }
1764 
1765  level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
1766  level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
1767  level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
1768  vflag[vbuffer[i] - vstart - nverts_prev] = 1;
1769  }
1770  }
1771  else if( type == MBQUAD )
1772  {
1773  double xi, eta, N[4];
1774  int findex = mbImpl->type_from_handle( *( _infaces.begin() ) ) - 1;
1775 
1776  for( int i = 4; i < vtotal; i++ )
1777  {
1778  if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
1779 
1780  xi = refTemplates[findex][d].vert_nat_coord[i - 4][0];
1781  eta = refTemplates[findex][d].vert_nat_coord[i - 4][1];
1782  N[0] = ( 1 - xi ) * ( 1 - eta ) / 4;
1783  N[1] = ( 1 + xi ) * ( 1 - eta ) / 4;
1784  N[2] = ( 1 + xi ) * ( 1 + eta ) / 4, N[3] = ( 1 - xi ) * ( 1 + eta ) / 4;
1785 
1786  double x = 0, y = 0, z = 0;
1787  for( int j = 0; j < 4; j++ )
1788  {
1789  x += N[j] * corner_coords[3 * j];
1790  y += N[j] * corner_coords[3 * j + 1];
1791  z += N[j] * corner_coords[3 * j + 2];
1792  }
1793 
1794  level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
1795  level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
1796  level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
1797  vflag[vbuffer[i] - vstart - nverts_prev] = 1;
1798  }
1799  }
1800  else if( type == MBTET )
1801  {
1802  double xi, eta, mu, N[4];
1803  int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
1804 
1805  for( int i = 4; i < vtotal; i++ )
1806  {
1807  if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
1808 
1809  xi = refTemplates[cindex][d].vert_nat_coord[i - 4][0];
1810  eta = refTemplates[cindex][d].vert_nat_coord[i - 4][1];
1811  mu = refTemplates[cindex][d].vert_nat_coord[i - 4][2];
1812 
1813  N[0] = 1 - xi - eta - mu;
1814  N[1] = xi;
1815  N[2] = eta, N[3] = mu;
1816 
1817  double x = 0, y = 0, z = 0;
1818  for( int j = 0; j < 4; j++ )
1819  {
1820  x += N[j] * corner_coords[3 * j];
1821  y += N[j] * corner_coords[3 * j + 1];
1822  z += N[j] * corner_coords[3 * j + 2];
1823  }
1824 
1825  level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
1826  level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
1827  level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
1828  vflag[vbuffer[i] - vstart - nverts_prev] = 1;
1829  }
1830  }
1831  else if( type == MBPRISM )
1832  {
1833  double xi, eta, mu, N[6];
1834  int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
1835 
1836  for( int i = 6; i < vtotal; i++ )
1837  {
1838  if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
1839 
1840  xi = refTemplates[cindex][d].vert_nat_coord[i - 6][0];
1841  eta = refTemplates[cindex][d].vert_nat_coord[i - 6][1];
1842  mu = refTemplates[cindex][d].vert_nat_coord[i - 6][2];
1843 
1844  N[0] = ( 1 - xi - eta ) * ( 1 - mu ), N[1] = xi * ( 1 - mu ), N[2] = eta * ( 1 - mu ),
1845  N[3] = ( 1 - xi - eta ) * ( 1 + mu ), N[4] = xi * ( 1 + mu ), N[5] = eta * ( 1 + mu );
1846 
1847  double x = 0, y = 0, z = 0;
1848  for( int j = 0; j < 6; j++ )
1849  {
1850  x += N[j] * corner_coords[3 * j];
1851  y += N[j] * corner_coords[3 * j + 1];
1852  z += N[j] * corner_coords[3 * j + 2];
1853  }
1854 
1855  level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
1856  level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
1857  level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
1858  vflag[vbuffer[i] - vstart - nverts_prev] = 1;
1859  }
1860  }
1861  else if( type == MBHEX )
1862  {
1863  double xi, eta, mu, N[8];
1864  double d1, d2, d3, s1, s2, s3;
1865  int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
1866 
1867  for( int i = 8; i < vtotal; i++ )
1868  {
1869 
1870  if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
1871 
1872  xi = refTemplates[cindex][d].vert_nat_coord[i - 8][0];
1873  eta = refTemplates[cindex][d].vert_nat_coord[i - 8][1];
1874  mu = refTemplates[cindex][d].vert_nat_coord[i - 8][2];
1875 
1876  d1 = 1 - xi;
1877  d2 = 1 - eta;
1878  d3 = 1 - mu;
1879  s1 = 1 + xi;
1880  s2 = 1 + eta;
1881  s3 = 1 + mu;
1882  N[0] = ( d1 * d2 * d3 ) / 8;
1883  N[1] = ( s1 * d2 * d3 ) / 8;
1884  N[2] = ( s1 * s2 * d3 ) / 8;
1885  N[3] = ( d1 * s2 * d3 ) / 8;
1886  N[4] = ( d1 * d2 * s3 ) / 8;
1887  N[5] = ( s1 * d2 * s3 ) / 8;
1888  N[6] = ( s1 * s2 * s3 ) / 8;
1889  N[7] = ( d1 * s2 * s3 ) / 8;
1890 
1891  double x = 0, y = 0, z = 0;
1892  for( int j = 0; j < 8; j++ )
1893  {
1894  x += N[j] * corner_coords[3 * j];
1895  y += N[j] * corner_coords[3 * j + 1];
1896  z += N[j] * corner_coords[3 * j + 2];
1897  }
1898 
1899  level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
1900  level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
1901  level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
1902 
1903  vflag[vbuffer[i] - vstart - nverts_prev] = 1;
1904  }
1905  }
1906  return MB_SUCCESS;
1907 }

References _incells, _infaces, moab::Range::begin(), moab::NestedRefine::level_memory::coordinates, get_index_from_degree(), level_mesh, MB_SUCCESS, MBHEX, mbImpl, MBPRISM, MBQUAD, MBTET, MBTRI, refTemplates, moab::NestedRefine::level_memory::start_vertex, moab::Core::type_from_handle(), and moab::NestedRefine::refPatterns::vert_nat_coord.

Referenced by construct_hm_2D(), subdivide_cells(), and subdivide_tets().

◆ construct_hm_1D() [1/2]

ErrorCode moab::NestedRefine::construct_hm_1D ( int  cur_level,
int  deg 
)
protected

Definition at line 804 of file NestedRefine.cpp.

805 {
807  int nverts_prev, nents_prev;
808  if( cur_level )
809  {
810  nverts_prev = level_mesh[cur_level - 1].num_verts;
811  nents_prev = level_mesh[cur_level - 1].num_edges;
812  }
813  else
814  {
815  nverts_prev = _inverts.size();
816  nents_prev = _inedges.size();
817  }
818 
819  int d = get_index_from_degree( deg );
820  int vtotal = 2 + refTemplates[0][d].total_new_verts;
821  std::vector< EntityHandle > vbuffer( vtotal );
822 
823  std::vector< EntityHandle > conn;
824  int count_nents = 0;
825  int count_verts = nverts_prev;
826 
827  // Step 1: Create the subentities via refinement of the previous mesh
828  for( int eid = 0; eid < nents_prev; eid++ )
829  {
830  conn.clear();
831 
832  // EntityHandle of the working edge
834  if( cur_level )
835  edge = level_mesh[cur_level - 1].start_edge + eid;
836  else
837  edge = _inedges[eid]; // Makes the assumption initial mesh is contiguous in memory
838 
839  error = get_connectivity( edge, cur_level, conn );MB_CHK_ERR( error );
840 
841  // Add the vertex handles to vbuffer for the current level for the working edge
842 
843  // Since the old vertices are copied first, their local indices do not change as new levels
844  // are added. Clearly the local indices of the new vertices introduced in the current level
845  // is still the same when the old vertices are copied. Thus, there is no need to explicitly
846  // store another map between the old and duplicates in the subsequent levels. The second
847  // part in the following sum is the local index in the previous level.
848 
849  // Add the corners to the vbuffer first.
850 
851  for( int i = 0; i < (int)conn.size(); i++ )
852  {
853  if( cur_level )
854  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
855  else
856  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
857  }
858 
859  // Adding rest of the entityhandles to working buffer for vertices.
860  int num_new_verts = refTemplates[0][d].total_new_verts;
861 
862  for( int i = 0; i < num_new_verts; i++ )
863  {
864  vbuffer[i + 2] = level_mesh[cur_level].start_vertex + count_verts;
865  count_verts += 1;
866  }
867 
868  // Use the template to obtain the subentities
869  int id1, id2;
870  int etotal = refTemplates[0][d].total_new_ents;
871  std::vector< EntityHandle > ent_buffer( etotal );
872 
873  for( int i = 0; i < etotal; i++ )
874  {
875  id1 = refTemplates[0][d].ents_conn[i][0];
876  id2 = refTemplates[0][d].ents_conn[i][1];
877  level_mesh[cur_level].edge_conn[2 * ( count_nents )] = vbuffer[id1];
878  level_mesh[cur_level].edge_conn[2 * ( count_nents ) + 1] = vbuffer[id2];
879  ent_buffer[i] = level_mesh[cur_level].start_edge + count_nents;
880  count_nents += 1;
881  };
882 
883  error = update_local_ahf( deg, MBEDGE, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
884 
885  // Compute the coordinates of the new vertices: Linear interpolation
886  int idx;
887  double xi;
888  for( int i = 0; i < num_new_verts; i++ )
889  {
890  xi = refTemplates[0][d].vert_nat_coord[i][0];
891  idx = vbuffer[i + 2] - level_mesh[cur_level].start_vertex; // index of new vertex in current level
892  if( cur_level )
893  {
894  id1 = conn[0] - level_mesh[cur_level - 1].start_vertex; // index of old end vertices in current level
895  id2 = conn[1] - level_mesh[cur_level - 1].start_vertex;
896  }
897  else
898  {
899  id1 = _inverts.index( conn[0] );
900  id2 = _inverts.index( conn[1] );
901  }
902 
903  level_mesh[cur_level].coordinates[0][idx] =
904  ( 1 - xi ) * level_mesh[cur_level].coordinates[0][id1] + xi * level_mesh[cur_level].coordinates[0][id2];
905  level_mesh[cur_level].coordinates[1][idx] =
906  ( 1 - xi ) * level_mesh[cur_level].coordinates[1][id1] + xi * level_mesh[cur_level].coordinates[1][id2];
907  level_mesh[cur_level].coordinates[2][idx] =
908  ( 1 - xi ) * level_mesh[cur_level].coordinates[2][id1] + xi * level_mesh[cur_level].coordinates[2][id2];
909  }
910  }
911 
912  error = update_global_ahf( MBEDGE, cur_level, deg );MB_CHK_ERR( error );
913 
914  return MB_SUCCESS;
915 }

References _inedges, _inverts, moab::Range::begin(), moab::NestedRefine::level_memory::coordinates, moab::NestedRefine::level_memory::edge_conn, moab::NestedRefine::refPatterns::ents_conn, moab::error(), ErrorCode, get_connectivity(), get_index_from_degree(), moab::Range::index(), level_mesh, MB_CHK_ERR, MB_SUCCESS, MBEDGE, moab::NestedRefine::level_memory::num_edges, moab::NestedRefine::level_memory::num_verts, refTemplates, moab::Range::size(), moab::NestedRefine::level_memory::start_edge, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, moab::NestedRefine::refPatterns::total_new_verts, update_global_ahf(), update_local_ahf(), and moab::NestedRefine::refPatterns::vert_nat_coord.

Referenced by construct_hm_2D(), construct_hm_entities(), subdivide_cells(), and subdivide_tets().

◆ construct_hm_1D() [2/2]

ErrorCode moab::NestedRefine::construct_hm_1D ( int  cur_level,
int  deg,
EntityType  type,
std::vector< EntityHandle > &  trackverts 
)
protected

Definition at line 917 of file NestedRefine.cpp.

921 {
923 
924  int nedges_prev;
925  if( cur_level )
926  nedges_prev = level_mesh[cur_level - 1].num_edges;
927  else
928  nedges_prev = _inedges.size();
929 
930  int d = get_index_from_degree( deg );
931  int nve = refTemplates[0][d].nv_edge;
932  int vtotal = 2 + refTemplates[0][d].total_new_verts;
933  int etotal = refTemplates[0][d].total_new_ents;
934  int ne = 0, dim = 0, index = 0;
935  if( type == MBTRI || type == MBQUAD )
936  {
937  index = type - 2;
938  ne = ahf->lConnMap2D[index].num_verts_in_face;
939  dim = 2;
940  }
941  else if( type == MBTET || type == MBHEX )
942  {
943  index = ahf->get_index_in_lmap( *( _incells.begin() ) );
944  ne = ahf->lConnMap3D[index].num_edges_in_cell;
945  dim = 3;
946  }
947 
948  std::vector< EntityHandle > vbuffer( vtotal );
949  std::vector< EntityHandle > ent_buffer( etotal );
950 
951  std::vector< EntityHandle > adjents, econn, fconn;
952  std::vector< int > leids;
953  int count_nents = 0;
954 
955  // Loop over all the edges and gather the vertices to be used for refinement
956  for( int eid = 0; eid < nedges_prev; eid++ )
957  {
958  adjents.clear();
959  leids.clear();
960  econn.clear();
961  fconn.clear();
962  for( int i = 0; i < vtotal; i++ )
963  vbuffer[i] = 0;
964  for( int i = 0; i < etotal; i++ )
965  ent_buffer[i] = 0;
966 
968  if( cur_level )
969  edge = level_mesh[cur_level - 1].start_edge + eid;
970  else
971  edge = _inedges[eid];
972 
973  error = get_connectivity( edge, cur_level, econn );MB_CHK_ERR( error );
974 
975  for( int i = 0; i < (int)econn.size(); i++ )
976  {
977  if( cur_level )
978  vbuffer[i] = level_mesh[cur_level].start_vertex + ( econn[i] - level_mesh[cur_level - 1].start_vertex );
979  else
980  vbuffer[i] = level_mesh[cur_level].start_vertex + ( econn[i] - *_inverts.begin() );
981  }
982 
983  int fid = -1, lid = -1, idx1 = -1, idx2 = -1;
984 
985  if( dim == 2 )
986  {
987  error = ahf->get_up_adjacencies_2d( edge, adjents, &leids );MB_CHK_ERR( error );
988  if( cur_level )
989  fid = adjents[0] - level_mesh[cur_level - 1].start_face;
990  else
991  fid = _infaces.index( adjents[0] );
992 
993  lid = leids[0];
994  idx1 = lid;
995  idx2 = ahf->lConnMap2D[index].next[lid];
996  }
997  else if( dim == 3 )
998  {
999  error = ahf->get_up_adjacencies_edg_3d( edge, adjents, &leids );MB_CHK_ERR( error );
1000  if( cur_level )
1001  fid = adjents[0] - level_mesh[cur_level - 1].start_cell;
1002  else
1003  fid = _incells.index( adjents[0] );
1004 
1005  lid = leids[0];
1006  idx1 = ahf->lConnMap3D[index].e2v[lid][0];
1007  idx2 = ahf->lConnMap3D[index].e2v[lid][1];
1008  }
1009 
1010  error = get_connectivity( adjents[0], cur_level, fconn );MB_CHK_ERR( error );
1011 
1012  bool orient = false;
1013  if( ( fconn[idx1] == econn[0] ) && ( fconn[idx2] == econn[1] ) ) orient = true;
1014 
1015  if( orient )
1016  {
1017  for( int j = 0; j < nve; j++ )
1018  vbuffer[j + 2] = trackverts[fid * ne * nve + nve * lid + j];
1019  }
1020  else
1021  {
1022  for( int j = 0; j < nve; j++ )
1023  vbuffer[( nve - j - 1 ) + 2] = trackverts[fid * ne * nve + nve * lid + j];
1024  }
1025 
1026  // Use the template to obtain the subentities
1027  int id1, id2;
1028 
1029  for( int i = 0; i < etotal; i++ )
1030  {
1031  id1 = refTemplates[0][d].ents_conn[i][0];
1032  id2 = refTemplates[0][d].ents_conn[i][1];
1033  level_mesh[cur_level].edge_conn[2 * ( count_nents )] = vbuffer[id1];
1034  level_mesh[cur_level].edge_conn[2 * ( count_nents ) + 1] = vbuffer[id2];
1035  ent_buffer[i] = level_mesh[cur_level].start_edge + count_nents;
1036 
1037  count_nents += 1;
1038  };
1039 
1040  error = update_local_ahf( deg, MBEDGE, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
1041  }
1042 
1043  error = update_global_ahf_1D_sub( cur_level, deg );MB_CHK_ERR( error );
1044 
1045  return MB_SUCCESS;
1046 }

References _incells, _inedges, _infaces, _inverts, ahf, moab::Range::begin(), dim, moab::HalfFacetRep::LocalMaps3D::e2v, moab::NestedRefine::level_memory::edge_conn, moab::NestedRefine::refPatterns::ents_conn, moab::error(), ErrorCode, get_connectivity(), get_index_from_degree(), moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::get_up_adjacencies_2d(), moab::HalfFacetRep::get_up_adjacencies_edg_3d(), moab::Range::index(), moab::HalfFacetRep::lConnMap2D, moab::HalfFacetRep::lConnMap3D, level_mesh, MB_CHK_ERR, MB_SUCCESS, MBEDGE, MBHEX, MBQUAD, MBTET, MBTRI, moab::HalfFacetRep::LocalMaps2D::next, moab::NestedRefine::level_memory::num_edges, moab::HalfFacetRep::LocalMaps3D::num_edges_in_cell, moab::HalfFacetRep::LocalMaps2D::num_verts_in_face, moab::NestedRefine::refPatterns::nv_edge, refTemplates, moab::Range::size(), moab::NestedRefine::level_memory::start_cell, moab::NestedRefine::level_memory::start_edge, moab::NestedRefine::level_memory::start_face, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, moab::NestedRefine::refPatterns::total_new_verts, update_global_ahf_1D_sub(), and update_local_ahf().

◆ construct_hm_2D() [1/2]

ErrorCode moab::NestedRefine::construct_hm_2D ( int  cur_level,
int  deg 
)
protected

Definition at line 1048 of file NestedRefine.cpp.

1049 {
1050  ErrorCode error;
1051  int nverts_prev, nents_prev;
1052  if( cur_level )
1053  {
1054  nverts_prev = level_mesh[cur_level - 1].num_verts;
1055  nents_prev = level_mesh[cur_level - 1].num_faces;
1056  }
1057  else
1058  {
1059  nverts_prev = _inverts.size();
1060  nents_prev = _infaces.size();
1061  }
1062 
1063  // Create some book-keeping arrays over the old mesh to avoid introducing duplicate vertices and
1064  // calculating vertices more than once.
1065  EntityType ftype = mbImpl->type_from_handle( *_infaces.begin() );
1066  int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
1067  int findex = ftype - 1;
1068 
1069  int d = get_index_from_degree( deg );
1070  int tnv = refTemplates[findex][d].total_new_verts;
1071  int vtotal = nepf + tnv;
1072  std::vector< EntityHandle > vbuffer( vtotal );
1073  int etotal = refTemplates[findex][d].total_new_ents;
1074  std::vector< EntityHandle > ent_buffer( etotal );
1075 
1076  int nve = refTemplates[findex][d].nv_edge;
1077  std::vector< EntityHandle > trackvertsF( nents_prev * nepf * nve, 0 );
1078  int cur_nverts = level_mesh[cur_level].num_verts;
1079  std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
1080 
1081  int count_nverts = nverts_prev;
1082  int count_nents = 0;
1083  std::vector< EntityHandle > conn, cur_conn;
1084 
1085  // Step 1: Create the subentities via refinement of the previous mesh
1086  for( int fid = 0; fid < nents_prev; fid++ )
1087  {
1088  conn.clear();
1089  cur_conn.clear();
1090  for( int i = 0; i < vtotal; i++ )
1091  vbuffer[i] = 0;
1092  for( int i = 0; i < etotal; i++ )
1093  ent_buffer[i] = 0;
1094 
1095  // EntityHandle of the working face
1097  if( cur_level )
1098  face = level_mesh[cur_level - 1].start_face + fid;
1099  else
1100  face = _infaces[fid];
1101 
1102  error = get_connectivity( face, cur_level, conn );MB_CHK_ERR( error );
1103 
1104  // Step 1: Add vertices from the current level for the working face that will be used for
1105  // subdivision.
1106  // Add the corners to vbuffer
1107  for( int i = 0; i < (int)conn.size(); i++ )
1108  {
1109  if( cur_level )
1110  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
1111  else
1112  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
1113 
1114  cur_conn.push_back( vbuffer[i] );
1115  }
1116 
1117  // Gather vertices already added to tracking array due to refinement of the sibling faces
1118 
1119  for( int i = 0; i < nepf; i++ )
1120  {
1121  for( int j = 0; j < nve; j++ )
1122  {
1123  int id = refTemplates[findex][d].vert_on_edges[i][j];
1124  vbuffer[id] = trackvertsF[fid * nve * nepf + nve * i + j];
1125  }
1126  }
1127 
1128  // Add the remaining vertex handles to vbuffer for the current level for the working face
1129  for( int i = 0; i < tnv; i++ )
1130  {
1131  if( !vbuffer[i + nepf] )
1132  {
1133  vbuffer[i + nepf] = level_mesh[cur_level].start_vertex + count_nverts;
1134  count_nverts += 1;
1135  }
1136  }
1137 
1138  // Step 2: Create the subentities using the template and the vbuffer
1139  int idx;
1140  for( int i = 0; i < etotal; i++ )
1141  {
1142  for( int k = 0; k < nepf; k++ )
1143  {
1144  idx = refTemplates[findex][d].ents_conn[i][k];
1145  level_mesh[cur_level].face_conn[nepf * count_nents + k] = vbuffer[idx];
1146  }
1147  ent_buffer[i] = level_mesh[cur_level].start_face + count_nents;
1148  count_nents += 1;
1149  }
1150 
1151  // Step 3: Update the local AHF maps
1152  error = update_local_ahf( deg, ftype, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
1153 
1154  // Step 4: Add the new vertices to the tracking array
1155  int id;
1156 
1157  for( int i = 0; i < nepf; i++ )
1158  {
1159  // Add the vertices to trackvertsF for fid
1160  for( int j = 0; j < nve; j++ )
1161  {
1162  id = refTemplates[findex][d].vert_on_edges[i][j];
1163  trackvertsF[fid * nepf * nve + nve * i + j] = vbuffer[id];
1164  }
1165 
1166  std::vector< EntityHandle > sibfids;
1167  std::vector< int > sibleids;
1168  std::vector< int > siborient;
1169 
1170  // Add the vertices to trackvertsF for siblings of fid, if any.
1171  error = ahf->get_up_adjacencies_2d( face, i, false, sibfids, &sibleids, &siborient );MB_CHK_ERR( error );
1172 
1173  if( !sibfids.size() ) continue;
1174 
1175  for( int s = 0; s < (int)sibfids.size(); s++ )
1176  {
1177  int sibid;
1178  if( cur_level )
1179  sibid = sibfids[s] - level_mesh[cur_level - 1].start_face;
1180  else
1181  sibid = sibfids[s] - *_infaces.begin();
1182 
1183  if( siborient[s] ) // Same half-edge direction as the current half-edge
1184  {
1185  for( int j = 0; j < nve; j++ )
1186  {
1187  id = refTemplates[findex][d].vert_on_edges[i][j];
1188  trackvertsF[sibid * nepf * nve + nve * sibleids[s] + j] = vbuffer[id];
1189  }
1190  }
1191  else
1192  {
1193  for( int j = 0; j < nve; j++ )
1194  {
1195  id = refTemplates[findex][d].vert_on_edges[i][nve - j - 1];
1196  trackvertsF[sibid * nepf * nve + nve * sibleids[s] + j] = vbuffer[id];
1197  }
1198  }
1199  }
1200  }
1201 
1202  // Step 5: Compute the coordinates of the new vertices, avoids computing more than once via
1203  // the flag_verts array.
1204  std::vector< double > corner_coords( nepf * 3 );
1205  error = get_coordinates( &cur_conn[0], nepf, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
1206 
1207  error = compute_coordinates( cur_level, deg, ftype, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
1208  nverts_prev );MB_CHK_ERR( error );
1209  }
1210 
1211  // Step 6: Update the global maps
1212  error = update_global_ahf( ftype, cur_level, deg );MB_CHK_ERR( error );
1213 
1214  // Step 7: If edges exists, refine them.
1215  if( !_inedges.empty() )
1216  {
1217  error = construct_hm_1D( cur_level, deg, ftype, trackvertsF );MB_CHK_ERR( error );
1218  }
1219 
1220  return MB_SUCCESS;
1221 }

References _inedges, _infaces, _inverts, ahf, moab::Range::begin(), compute_coordinates(), construct_hm_1D(), moab::Range::empty(), moab::NestedRefine::refPatterns::ents_conn, moab::error(), ErrorCode, moab::NestedRefine::level_memory::face_conn, get_connectivity(), get_coordinates(), get_index_from_degree(), moab::HalfFacetRep::get_up_adjacencies_2d(), moab::HalfFacetRep::lConnMap2D, level_mesh, MB_CHK_ERR, MB_SUCCESS, mbImpl, moab::NestedRefine::level_memory::num_faces, moab::NestedRefine::level_memory::num_verts, moab::HalfFacetRep::LocalMaps2D::num_verts_in_face, moab::NestedRefine::refPatterns::nv_edge, refTemplates, moab::Range::size(), moab::NestedRefine::level_memory::start_face, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, moab::NestedRefine::refPatterns::total_new_verts, moab::Core::type_from_handle(), update_global_ahf(), update_local_ahf(), and moab::NestedRefine::refPatterns::vert_on_edges.

Referenced by construct_hm_entities(), subdivide_cells(), and subdivide_tets().

◆ construct_hm_2D() [2/2]

ErrorCode moab::NestedRefine::construct_hm_2D ( int  cur_level,
int  deg,
EntityType  type,
std::vector< EntityHandle > &  trackvertsC_edg,
std::vector< EntityHandle > &  trackvertsF 
)
protected

Definition at line 1223 of file NestedRefine.cpp.

1228 {
1229  ErrorCode error;
1230 
1231  EntityType ftype = MBTRI;
1232  if( type == MBHEX ) ftype = MBQUAD;
1233 
1234  int d = get_index_from_degree( deg );
1235  int findex = ftype - 1;
1236  int cidx = ahf->get_index_in_lmap( *( _incells.begin() ) );
1237 
1238  int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
1239  int nepc = ahf->lConnMap3D[cidx].num_edges_in_cell;
1240  int nfpc = ahf->lConnMap3D[cidx].num_faces_in_cell;
1241 
1242  int tnv = refTemplates[findex][d].total_new_verts;
1243  int nve = refTemplates[findex][d].nv_edge;
1244  int nvf = refTemplates[findex][d].nv_face;
1245  int vtotal = nepf + tnv;
1246  int etotal = refTemplates[findex][d].total_new_ents;
1247 
1248  std::vector< EntityHandle > vbuffer( vtotal );
1249  std::vector< EntityHandle > ent_buffer( etotal );
1250 
1251  std::vector< EntityHandle > adjents, fconn, cconn;
1252  std::vector< int > leids;
1253  int count_nents = 0;
1254 
1255  int nents_prev, ecount;
1256  if( cur_level )
1257  {
1258  nents_prev = level_mesh[cur_level - 1].num_faces;
1259  ecount = level_mesh[cur_level - 1].num_edges * refTemplates[MBEDGE - 1][d].total_new_ents;
1260  ;
1261  }
1262  else
1263  {
1264  nents_prev = _infaces.size();
1265  ecount = _inedges.size() * refTemplates[MBEDGE - 1][d].total_new_ents;
1266  ;
1267  }
1268 
1269  // Step 1: Create the subentities via refinement of the previous mesh
1270  for( int it = 0; it < nents_prev; it++ )
1271  {
1272  fconn.clear();
1273  cconn.clear();
1274  adjents.clear();
1275  leids.clear();
1276  for( int i = 0; i < vtotal; i++ )
1277  vbuffer[i] = 0;
1278  for( int i = 0; i < etotal; i++ )
1279  ent_buffer[i] = 0;
1280 
1281  // EntityHandle of the working face
1283  if( cur_level )
1284  face = level_mesh[cur_level - 1].start_face + it;
1285  else
1286  face = _infaces[it];
1287 
1288  error = get_connectivity( face, cur_level, fconn );MB_CHK_ERR( error );
1289 
1290  // Add the new handles for old connectivity in the buffer
1291  for( int i = 0; i < (int)fconn.size(); i++ )
1292  {
1293  if( cur_level )
1294  vbuffer[i] = level_mesh[cur_level].start_vertex + ( fconn[i] - level_mesh[cur_level - 1].start_vertex );
1295  else
1296  vbuffer[i] = level_mesh[cur_level].start_vertex + ( fconn[i] - *_inverts.begin() );
1297  }
1298 
1299  // Add handles for vertices on edges and faces from the already refined cell
1300  int fid, lid;
1301  error = ahf->get_up_adjacencies_face_3d( face, adjents, &leids );MB_CHK_ERR( error );
1302 
1303  if( cur_level )
1304  fid = adjents[0] - level_mesh[cur_level - 1].start_cell;
1305  else
1306  fid = _incells.index( adjents[0] );
1307 
1308  lid = leids[0];
1309 
1310  error = get_connectivity( adjents[0], cur_level, cconn );MB_CHK_ERR( error );
1311 
1312  // Find the orientation w.r.t the half-face and then add vertices properly.
1313  std::vector< EntityHandle > fac_conn( nepf );
1314  std::vector< EntityHandle > lfac_conn( nepf );
1315  for( int j = 0; j < nepf; j++ )
1316  {
1317  fac_conn[j] = fconn[j];
1318  int id = ahf->lConnMap3D[cidx].hf2v[lid][j];
1319  lfac_conn[j] = cconn[id];
1320  }
1321 
1322  std::vector< int > le_idx, indices;
1323 
1324  error = reorder_indices( deg, &fac_conn[0], &lfac_conn[0], nepf, le_idx, indices );MB_CHK_ERR( error );
1325 
1326  // Add the existing vertices on edges of the already refined cell to the vbuffer
1327  for( int j = 0; j < nepf; j++ )
1328  {
1329  int id = le_idx[j]; // Corresponding local edge
1330  int idx = ahf->lConnMap3D[cidx].f2leid[lid][id]; // Local edge in the cell
1331 
1332  // Get the orientation of the local edge of the face wrt the corresponding local edge in
1333  // the cell
1334  bool eorient = false;
1335  int fnext = ahf->lConnMap2D[ftype - 2].next[j];
1336  int idx1 = ahf->lConnMap3D[cidx].e2v[idx][0];
1337  int idx2 = ahf->lConnMap3D[cidx].e2v[idx][1];
1338  if( ( fconn[j] == cconn[idx1] ) && ( fconn[fnext] == cconn[idx2] ) ) eorient = true;
1339 
1340  if( eorient )
1341  {
1342  for( int k = 0; k < nve; k++ )
1343  {
1344  int ind = refTemplates[findex][d].vert_on_edges[j][k];
1345  vbuffer[ind] = trackvertsE[fid * nepc * nve + nve * idx + k];
1346  }
1347  }
1348  else
1349  {
1350  for( int k = 0; k < nve; k++ )
1351  {
1352  int ind = refTemplates[findex][d].vert_on_edges[j][nve - k - 1];
1353  vbuffer[ind] = trackvertsE[fid * nepc * nve + nve * idx + k];
1354  }
1355  }
1356  }
1357 
1358  // Add the existing vertices on the face of the refine cell to vbuffer
1359  if( nvf )
1360  {
1361  for( int k = 0; k < nvf; k++ )
1362  {
1363  int ind = refTemplates[findex][d].vert_on_faces[0][k];
1364  vbuffer[ind] = trackvertsF[fid * nfpc * nvf + nvf * lid + indices[k] - 1];
1365  }
1366  }
1367 
1368  // Create the subentities using the template and the vbuffer
1369  for( int i = 0; i < etotal; i++ )
1370  {
1371  for( int k = 0; k < nepf; k++ )
1372  {
1373  int idx = refTemplates[findex][d].ents_conn[i][k];
1374  level_mesh[cur_level].face_conn[nepf * count_nents + k] = vbuffer[idx];
1375  }
1376  ent_buffer[i] = level_mesh[cur_level].start_face + count_nents;
1377  count_nents += 1;
1378  }
1379 
1380  error = update_local_ahf( deg, ftype, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
1381 
1382  // Create the interior edges
1383  int id1, id2;
1384 
1385  int ne = intFacEdg[ftype - 2][d].nie;
1386  for( int i = 0; i < ne; i++ )
1387  {
1388  id1 = intFacEdg[ftype - 2][d].ieconn[i][0];
1389  id2 = intFacEdg[ftype - 2][d].ieconn[i][1];
1390  level_mesh[cur_level].edge_conn[2 * ( ecount )] = vbuffer[id1];
1391  level_mesh[cur_level].edge_conn[2 * ( ecount ) + 1] = vbuffer[id2];
1392  ecount += 1;
1393  }
1394  }
1395 
1396  // Step 6: Update the global maps
1397  error = update_global_ahf_2D_sub( cur_level, deg );MB_CHK_ERR( error );
1398 
1399  // Step 7: Update the hf-maps for the edges
1400  error = update_ahf_1D( cur_level );MB_CHK_ERR( error );
1401 
1402  return MB_SUCCESS;
1403 }

References _incells, _inedges, _infaces, _inverts, ahf, moab::Range::begin(), moab::HalfFacetRep::LocalMaps3D::e2v, moab::NestedRefine::level_memory::edge_conn, moab::NestedRefine::refPatterns::ents_conn, moab::error(), ErrorCode, moab::HalfFacetRep::LocalMaps3D::f2leid, moab::NestedRefine::level_memory::face_conn, get_connectivity(), get_index_from_degree(), moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::get_up_adjacencies_face_3d(), moab::HalfFacetRep::LocalMaps3D::hf2v, moab::NestedRefine::intFEdge::ieconn, moab::Range::index(), intFacEdg, moab::HalfFacetRep::lConnMap2D, moab::HalfFacetRep::lConnMap3D, level_mesh, MB_CHK_ERR, MB_SUCCESS, MBEDGE, MBHEX, MBQUAD, MBTRI, moab::HalfFacetRep::LocalMaps2D::next, moab::NestedRefine::intFEdge::nie, moab::NestedRefine::level_memory::num_edges, moab::HalfFacetRep::LocalMaps3D::num_edges_in_cell, moab::NestedRefine::level_memory::num_faces, moab::HalfFacetRep::LocalMaps3D::num_faces_in_cell, moab::HalfFacetRep::LocalMaps2D::num_verts_in_face, moab::NestedRefine::refPatterns::nv_edge, moab::NestedRefine::refPatterns::nv_face, refTemplates, reorder_indices(), moab::Range::size(), moab::NestedRefine::level_memory::start_cell, moab::NestedRefine::level_memory::start_face, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, moab::NestedRefine::refPatterns::total_new_verts, update_ahf_1D(), update_global_ahf_2D_sub(), update_local_ahf(), moab::NestedRefine::refPatterns::vert_on_edges, and moab::NestedRefine::refPatterns::vert_on_faces.

◆ construct_hm_3D()

ErrorCode moab::NestedRefine::construct_hm_3D ( int  cur_level,
int  deg 
)
protected

Definition at line 1405 of file NestedRefine.cpp.

1406 {
1407  ErrorCode error;
1408  EntityType type = mbImpl->type_from_handle( *( _incells.begin() ) );
1409  if( type == MBTET )
1410  {
1411  error = subdivide_tets( cur_level, deg );MB_CHK_ERR( error );
1412  }
1413  else
1414  {
1415  error = subdivide_cells( type, cur_level, deg );MB_CHK_ERR( error );
1416  }
1417 
1418  return MB_SUCCESS;
1419 }

References _incells, moab::Range::begin(), moab::error(), ErrorCode, MB_CHK_ERR, MB_SUCCESS, mbImpl, MBTET, subdivide_cells(), subdivide_tets(), and moab::Core::type_from_handle().

Referenced by construct_hm_entities().

◆ construct_hm_entities()

ErrorCode moab::NestedRefine::construct_hm_entities ( int  cur_level,
int  deg 
)
protected

Definition at line 783 of file NestedRefine.cpp.

784 {
786 
787  // Generate mesh for current level by refining previous level.
788  if( ahf->thismeshtype == CURVE )
789  {
790  error = construct_hm_1D( cur_level, deg );MB_CHK_ERR( error );
791  }
792  else if( ahf->thismeshtype == SURFACE || ahf->thismeshtype == SURFACE_MIXED )
793  {
794  error = construct_hm_2D( cur_level, deg );MB_CHK_ERR( error );
795  }
796  else
797  {
798  error = construct_hm_3D( cur_level, deg );MB_CHK_ERR( error );
799  }
800 
801  return MB_SUCCESS;
802 }

References ahf, construct_hm_1D(), construct_hm_2D(), construct_hm_3D(), moab::CURVE, moab::error(), ErrorCode, MB_CHK_ERR, MB_SUCCESS, moab::SURFACE, moab::SURFACE_MIXED, and moab::HalfFacetRep::thismeshtype.

Referenced by generate_hm().

◆ copy_vertices_from_prev_level()

ErrorCode moab::NestedRefine::copy_vertices_from_prev_level ( int  cur_level)
protected

Definition at line 4165 of file NestedRefine.cpp.

4166 {
4167  ErrorCode error;
4168 
4169  if( cur_level )
4170  {
4171  int nverts_prev = level_mesh[cur_level - 1].num_verts;
4172  for( int i = 0; i < nverts_prev; i++ )
4173  {
4174  level_mesh[cur_level].coordinates[0][i] = level_mesh[cur_level - 1].coordinates[0][i];
4175  level_mesh[cur_level].coordinates[1][i] = level_mesh[cur_level - 1].coordinates[1][i];
4176  level_mesh[cur_level].coordinates[2][i] = level_mesh[cur_level - 1].coordinates[2][i];
4177  }
4178  }
4179  else // Copy the vertices from the input mesh
4180  {
4181  int nverts_in = _inverts.size();
4182  std::vector< double > vcoords( 3 * nverts_in );
4183  error = mbImpl->get_coords( _inverts, &vcoords[0] );MB_CHK_ERR( error );
4184 
4185  for( int i = 0; i < nverts_in; i++ )
4186  {
4187  level_mesh[cur_level].coordinates[0][i] = vcoords[3 * i];
4188  level_mesh[cur_level].coordinates[1][i] = vcoords[3 * i + 1];
4189  level_mesh[cur_level].coordinates[2][i] = vcoords[3 * i + 2];
4190  }
4191  }
4192  return MB_SUCCESS;
4193  // To add: Map from old vertices to new duplicates: NOT NEEDED
4194 }

References _inverts, moab::NestedRefine::level_memory::coordinates, moab::error(), ErrorCode, moab::Core::get_coords(), level_mesh, MB_CHK_ERR, MB_SUCCESS, mbImpl, moab::NestedRefine::level_memory::num_verts, and moab::Range::size().

Referenced by generate_hm().

◆ count_subentities()

ErrorCode moab::NestedRefine::count_subentities ( EntityHandle  set,
int  cur_level,
int *  nedges,
int *  nfaces 
)
protected

Definition at line 4520 of file NestedRefine.cpp.

4521 {
4522  ErrorCode error;
4523 
4524  if( cur_level >= 0 )
4525  {
4526  Range edges, faces, cells;
4527 
4528  error = mbImpl->get_entities_by_dimension( set, 1, edges );MB_CHK_ERR( error );
4529 
4530  error = mbImpl->get_entities_by_dimension( set, 2, faces );MB_CHK_ERR( error );
4531 
4532  error = mbImpl->get_entities_by_dimension( set, 3, cells );MB_CHK_ERR( error );
4533 
4534  error = ahf->count_subentities( edges, faces, cells, nedges, nfaces );MB_CHK_ERR( error );
4535  }
4536  else
4537  {
4538  error = ahf->count_subentities( _inedges, _infaces, _incells, nedges, nfaces );MB_CHK_ERR( error );
4539  }
4540 
4541  return MB_SUCCESS;
4542 }

References _incells, _inedges, _infaces, ahf, moab::HalfFacetRep::count_subentities(), moab::error(), ErrorCode, moab::Core::get_entities_by_dimension(), MB_CHK_ERR, MB_SUCCESS, and mbImpl.

Referenced by estimate_hm_storage().

◆ create_hm_storage_single_level()

ErrorCode moab::NestedRefine::create_hm_storage_single_level ( EntityHandle set,
int  cur_level,
int  estL[4] 
)
protected

Definition at line 624 of file NestedRefine.cpp.

625 {
626  // Obtain chunks of memory for the current level. Add them to a particular meshset.
627  EntityHandle set_handle;
628  ErrorCode error = mbImpl->create_meshset( MESHSET_SET, set_handle );MB_CHK_SET_ERR( error, "Cannot create mesh for the current level" );
629  *set = set_handle;
630 
631  ReadUtilIface* read_iface;
632  error = mbImpl->query_interface( read_iface );MB_CHK_ERR( error );
633 
634  // Vertices
635  error = read_iface->get_node_coords( 3, estL[0], 0, level_mesh[cur_level].start_vertex,
636  level_mesh[cur_level].coordinates );MB_CHK_ERR( error );
637  level_mesh[cur_level].num_verts = estL[0];
638 
639  Range newverts( level_mesh[cur_level].start_vertex, level_mesh[cur_level].start_vertex + estL[0] - 1 );
640  error = mbImpl->add_entities( *set, newverts );MB_CHK_ERR( error );
641  level_mesh[cur_level].verts = newverts;
642 
643  Tag gidtag;
645  error = read_iface->assign_ids( gidtag, newverts, level_mesh[cur_level].start_vertex );MB_CHK_ERR( error );
646 
647  // Edges
648  if( estL[1] )
649  {
650  error = read_iface->get_element_connect( estL[1], 2, MBEDGE, 0, level_mesh[cur_level].start_edge,
651  level_mesh[cur_level].edge_conn );MB_CHK_ERR( error );
652  level_mesh[cur_level].num_edges = estL[1];
653 
654  Range newedges( level_mesh[cur_level].start_edge, level_mesh[cur_level].start_edge + estL[1] - 1 );
655  error = mbImpl->add_entities( *set, newedges );MB_CHK_ERR( error );
656  level_mesh[cur_level].edges = newedges;
657  }
658  else
659  level_mesh[cur_level].num_edges = 0;
660 
661  // Faces
662  if( estL[2] )
663  {
664  EntityType type = mbImpl->type_from_handle( *( _infaces.begin() ) );
665  int nvpf = ahf->lConnMap2D[type - 2].num_verts_in_face;
666  error = read_iface->get_element_connect( estL[2], nvpf, type, 0, level_mesh[cur_level].start_face,
667  level_mesh[cur_level].face_conn );MB_CHK_ERR( error );
668  level_mesh[cur_level].num_faces = estL[2];
669 
670  Range newfaces( level_mesh[cur_level].start_face, level_mesh[cur_level].start_face + estL[2] - 1 );
671  error = mbImpl->add_entities( *set, newfaces );MB_CHK_ERR( error );
672  level_mesh[cur_level].faces = newfaces;
673  }
674  else
675  level_mesh[cur_level].num_faces = 0;
676 
677  // Cells
678  if( estL[3] )
679  {
680  EntityType type = mbImpl->type_from_handle( *( _incells.begin() ) );
681  int index = ahf->get_index_in_lmap( *_incells.begin() );
682  int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
683  error = read_iface->get_element_connect( estL[3], nvpc, type, 0, level_mesh[cur_level].start_cell,
684  level_mesh[cur_level].cell_conn );MB_CHK_ERR( error );
685  level_mesh[cur_level].num_cells = estL[3];
686 
687  Range newcells( level_mesh[cur_level].start_cell, level_mesh[cur_level].start_cell + estL[3] - 1 );
688  error = mbImpl->add_entities( *set, newcells );MB_CHK_ERR( error );
689  level_mesh[cur_level].cells = newcells;
690  }
691  else
692  level_mesh[cur_level].num_cells = 0;
693 
694  // Resize the ahf maps
695  error = ahf->resize_hf_maps( level_mesh[cur_level].start_vertex, level_mesh[cur_level].num_verts,
696  level_mesh[cur_level].start_edge, level_mesh[cur_level].num_edges,
697  level_mesh[cur_level].start_face, level_mesh[cur_level].num_faces,
698  level_mesh[cur_level].start_cell, level_mesh[cur_level].num_cells );MB_CHK_ERR( error );
699 
701 
702  // If the mesh type changes, then update the member variable in ahf to use the applicable
703  // adjacency matrix
704  MESHTYPE nwmesh = ahf->get_mesh_type( level_mesh[cur_level].num_verts, level_mesh[cur_level].num_edges,
705  level_mesh[cur_level].num_faces, level_mesh[cur_level].num_cells );MB_CHK_ERR( error );
706  if( ahf->thismeshtype != nwmesh ) ahf->thismeshtype = nwmesh;
707 
708  return MB_SUCCESS;
709 }

References _incells, _infaces, moab::Core::add_entities(), ahf, moab::ReadUtilIface::assign_ids(), moab::Range::begin(), moab::NestedRefine::level_memory::cells, moab::Core::create_meshset(), moab::NestedRefine::level_memory::edges, moab::error(), ErrorCode, moab::NestedRefine::level_memory::faces, moab::ReadUtilIface::get_element_connect(), moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::get_mesh_type(), moab::ReadUtilIface::get_node_coords(), GLOBAL_ID_TAG_NAME, moab::HalfFacetRep::lConnMap2D, moab::HalfFacetRep::lConnMap3D, level_mesh, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MBEDGE, mbImpl, MESHSET_SET, moab::NestedRefine::level_memory::num_cells, moab::NestedRefine::level_memory::num_edges, moab::NestedRefine::level_memory::num_faces, moab::NestedRefine::level_memory::num_verts, moab::HalfFacetRep::LocalMaps3D::num_verts_in_cell, moab::HalfFacetRep::LocalMaps2D::num_verts_in_face, moab::Interface::query_interface(), moab::HalfFacetRep::resize_hf_maps(), moab::Core::tag_get_handle(), moab::HalfFacetRep::thismeshtype, moab::Core::type_from_handle(), moab::HalfFacetRep::update_entity_ranges(), and moab::NestedRefine::level_memory::verts.

Referenced by generate_hm().

◆ estimate_hm_storage()

ErrorCode moab::NestedRefine::estimate_hm_storage ( EntityHandle  set,
int  level_degree,
int  cur_level,
int  hmest[4] 
)
protected

Definition at line 564 of file NestedRefine.cpp.

565 {
567 
568  // Obtain the size of input mesh.
569  int nverts_prev, nedges_prev, nfaces_prev, ncells_prev;
570  if( cur_level )
571  {
572  nverts_prev = level_mesh[cur_level - 1].num_verts;
573  nedges_prev = level_mesh[cur_level - 1].num_edges;
574  nfaces_prev = level_mesh[cur_level - 1].num_faces;
575  ncells_prev = level_mesh[cur_level - 1].num_cells;
576  }
577  else
578  {
579  nverts_prev = _inverts.size();
580  nedges_prev = _inedges.size();
581  nfaces_prev = _infaces.size();
582  ncells_prev = _incells.size();
583  }
584 
585  // Estimate mesh size of current level mesh.
586  int nedges = 0, nfaces = 0;
587  error = count_subentities( set, cur_level - 1, &nedges, &nfaces );MB_CHK_ERR( error );
588 
589  int d = get_index_from_degree( level_degree );
590  int nverts = refTemplates[MBEDGE - 1][d].nv_edge * nedges;
591  hmest[0] = nverts_prev + nverts;
592  hmest[1] = nedges_prev * refTemplates[MBEDGE - 1][d].total_new_ents;
593  hmest[2] = 0;
594  hmest[3] = 0;
595 
596  int findex, cindex;
597  if( nfaces_prev != 0 )
598  {
599  EntityHandle start_face;
600  if( cur_level )
601  start_face = level_mesh[cur_level - 1].start_face;
602  else
603  start_face = *_infaces.begin();
604  findex = mbImpl->type_from_handle( start_face ) - 1;
605  hmest[2] = nfaces_prev * refTemplates[findex][d].total_new_ents;
606 
607  if( meshdim == 2 ) hmest[0] += refTemplates[findex][d].nv_face * nfaces_prev;
608 
609  if( meshdim == 3 ) hmest[1] += nfaces_prev * intFacEdg[findex - 1][d].nie;
610  }
611 
612  if( ncells_prev != 0 )
613  {
614  cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
615  hmest[3] = ncells_prev * refTemplates[cindex][d].total_new_ents;
616 
617  hmest[0] += refTemplates[cindex][d].nv_face * nfaces;
618  hmest[0] += refTemplates[cindex][d].nv_cell * ncells_prev;
619  }
620 
621  return MB_SUCCESS;
622 }

References _incells, _inedges, _infaces, _inverts, moab::Range::begin(), count_subentities(), moab::error(), ErrorCode, get_index_from_degree(), intFacEdg, level_mesh, MB_CHK_ERR, MB_SUCCESS, MBEDGE, mbImpl, meshdim, moab::NestedRefine::intFEdge::nie, moab::NestedRefine::level_memory::num_cells, moab::NestedRefine::level_memory::num_edges, moab::NestedRefine::level_memory::num_faces, moab::NestedRefine::level_memory::num_verts, moab::NestedRefine::refPatterns::nv_cell, moab::NestedRefine::refPatterns::nv_edge, moab::NestedRefine::refPatterns::nv_face, refTemplates, moab::Range::size(), moab::NestedRefine::level_memory::start_face, moab::NestedRefine::refPatterns::total_new_ents, and moab::Core::type_from_handle().

Referenced by generate_hm().

◆ exchange_ghosts()

ErrorCode moab::NestedRefine::exchange_ghosts ( std::vector< EntityHandle > &  lsets,
int  num_glayers 
)

Definition at line 454 of file NestedRefine.cpp.

455 {
457 
458  if( hasghost ) return MB_SUCCESS;
459 
460  hasghost = true;
461 #ifdef MOAB_HAVE_MPI
462  error = pcomm->exchange_ghost_cells( meshdim, 0, num_glayers, 0, true, false );MB_CHK_ERR( error );
463  {
464  Range empty_range;
466  // error = pcomm->assign_global_ids(lsets[i], 0, 1, false, true, false);MB_CHK_ERR(error);
467  }
468 #else
469  MB_SET_ERR( MB_FAILURE, "Requesting ghost layers for a serial mesh" );
470 #endif
471 
472  Range* lverts = new Range[lsets.size()];
473  Range* lents = new Range[lsets.size()];
474  for( size_t i = 0; i < lsets.size(); i++ )
475  {
476  error = mbImpl->get_entities_by_dimension( lsets[i], meshdim, lents[i] );MB_CHK_ERR( error );
477  error = mbImpl->get_connectivity( lents[i], lverts[i] );MB_CHK_ERR( error );
478 
479  for( int gl = 0; gl < num_glayers; gl++ )
480  {
481  error = mbImpl->get_adjacencies( lverts[i], meshdim, false, lents[i], Interface::UNION );MB_CHK_ERR( error );
482  error = mbImpl->get_connectivity( lents[i], lverts[i] );MB_CHK_ERR( error );
483  }
484  }
485  for( size_t i = 0; i < lsets.size(); i++ )
486  {
487  error = mbImpl->add_entities( lsets[i], lverts[i] );MB_CHK_ERR( error );
488  error = mbImpl->add_entities( lsets[i], lents[i] );MB_CHK_ERR( error );
489  }
490 
491  delete[] lverts;
492  delete[] lents;
493  return MB_SUCCESS;
494 }

References moab::Core::add_entities(), moab::error(), ErrorCode, moab::ParallelComm::exchange_ghost_cells(), moab::ParallelComm::exchange_tags(), moab::Core::get_adjacencies(), moab::Core::get_connectivity(), moab::Core::get_entities_by_dimension(), GLOBAL_ID_TAG_NAME, hasghost, MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, mbImpl, meshdim, pcomm, and moab::Interface::UNION.

Referenced by main().

◆ find_shortest_diagonal_octahedron()

int moab::NestedRefine::find_shortest_diagonal_octahedron ( int  cur_level,
int  deg,
EntityHandle vbuffer 
)
protected

Definition at line 4580 of file NestedRefine.cpp.

4581 {
4582  ErrorCode error;
4583  double coords[18];
4584  error = get_octahedron_corner_coords( cur_level, deg, vbuffer, coords );
4585  if( error != MB_SUCCESS ) MB_SET_ERR( MB_FAILURE, "Error in obtaining octahedron corner coordinates" );
4586 
4587  int diag_map[6] = { 1, 3, 2, 4, 5, 0 };
4588  double length = std::numeric_limits< double >::max();
4589 
4590  int diag = 0;
4591  double x, y, z;
4592  x = y = z = 0;
4593 
4594  for( int d = 0; d < 3; d++ )
4595  {
4596  int id1 = diag_map[2 * d];
4597  int id2 = diag_map[2 * d + 1];
4598  x = coords[3 * id1] - coords[3 * id2];
4599  y = coords[3 * id1 + 1] - coords[3 * id2 + 1];
4600  z = coords[3 * id1 + 2] - coords[3 * id2 + 2];
4601  double dlen = sqrt( x * x + y * y + z * z );
4602  if( dlen < length )
4603  {
4604  length = dlen;
4605  diag = d + 1;
4606  }
4607  }
4608 
4609  return diag;
4610 }

References moab::error(), ErrorCode, get_octahedron_corner_coords(), length(), MB_SET_ERR, and MB_SUCCESS.

Referenced by subdivide_tets().

◆ generate_hm()

ErrorCode moab::NestedRefine::generate_hm ( int *  level_degrees,
int  num_level,
EntityHandle hm_set,
bool  optimize 
)
protected

Definition at line 715 of file NestedRefine.cpp.

716 {
718 
719  Tag gidtag;
721 
722  nlevels = num_level;
723 
724  timeall.tm_total = 0;
725  timeall.tm_refine = 0;
726  timeall.tm_resolve = 0;
727 
728  for( int l = 0; l < num_level; l++ )
729  {
730  double tstart;
731  tstart = tm->time_elapsed();
732 
733  // Estimate storage
734  int hmest[4] = { 0, 0, 0, 0 };
735  EntityHandle set;
736  if( l )
737  set = hm_set[l - 1];
738  else
739  set = _rset;
740  error = estimate_hm_storage( set, level_degrees[l], l, hmest );MB_CHK_ERR( error );
741 
742  // Create arrays for storing the current level
743  error = create_hm_storage_single_level( &hm_set[l], l, hmest );MB_CHK_ERR( error );
744 
745  // Copy the old vertices along with their coordinates
747 
748  // Create the new entities and new vertices
749  error = construct_hm_entities( l, level_degrees[l] );MB_CHK_ERR( error );
750 
751  timeall.tm_refine += tm->time_elapsed() - tstart;
752 
753  // Go into parallel communication
754  if( !optimize )
755  {
756 #ifdef MOAB_HAVE_MPI
757  if( pcomm && ( pcomm->size() > 1 ) )
758  {
759  double tpstart = tm->time_elapsed();
760  error = resolve_shared_ents_parmerge( l, hm_set[l] );MB_CHK_ERR( error );
761  timeall.tm_resolve += tm->time_elapsed() - tpstart;
762  }
763 #endif
764  }
765  }
766 
767  if( optimize )
768  {
769 #ifdef MOAB_HAVE_MPI
770  if( pcomm && ( pcomm->size() > 1 ) )
771  {
772  double tpstart = tm->time_elapsed();
773  error = resolve_shared_ents_opt( hm_set, nlevels );MB_CHK_ERR( error );
774  timeall.tm_resolve = tm->time_elapsed() - tpstart;
775  }
776 #endif
777  }
779 
780  return MB_SUCCESS;
781 }

References _rset, construct_hm_entities(), copy_vertices_from_prev_level(), create_hm_storage_single_level(), moab::error(), ErrorCode, estimate_hm_storage(), GLOBAL_ID_TAG_NAME, MB_CHK_ERR, MB_SUCCESS, mbImpl, nlevels, pcomm, moab::ParallelComm::size(), moab::Core::tag_get_handle(), moab::CpuTimer::time_elapsed(), timeall, tm, moab::NestedRefine::codeperf::tm_refine, moab::NestedRefine::codeperf::tm_resolve, and moab::NestedRefine::codeperf::tm_total.

Referenced by generate_mesh_hierarchy().

◆ generate_mesh_hierarchy()

ErrorCode moab::NestedRefine::generate_mesh_hierarchy ( int  num_level,
int *  level_degrees,
std::vector< EntityHandle > &  level_sets,
bool  optimize = false 
)

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.

Parameters
level_degreesInteger array storing the degrees used in each level.
num_levelThe total number of levels in the hierarchy.
hm_setEntityHandle STL vector that returns the handles of the sets created for each mesh level.

Definition at line 114 of file NestedRefine.cpp.

118 {
119  assert( num_level > 0 );
120  nlevels = num_level;
121 
123  std::vector< moab::EntityHandle > hmsets( num_level );
124 
125  if( meshdim <= 2 )
126  {
127  for( int i = 0; i < num_level; i++ )
128  {
129  assert( ( level_degrees[i] == 2 ) || ( level_degrees[i] == 3 ) || ( level_degrees[i] == 5 ) );
130  level_dsequence[i] = level_degrees[i];
131  }
132  }
133  else
134  {
135  for( int i = 0; i < num_level; i++ )
136  {
137  assert( ( level_degrees[i] == 2 ) || ( level_degrees[i] == 3 ) );
138  level_dsequence[i] = level_degrees[i];
139  }
140  }
141 
142  error = generate_hm( level_degrees, num_level, &hmsets[0], optimize );MB_CHK_ERR( error );
143 
144  // copy the entity handles
145  level_sets.resize( num_level + 1 );
146  level_sets[0] = _rset;
147  for( int i = 0; i < num_level; i++ )
148  level_sets[i + 1] = hmsets[i];
149 
150  return MB_SUCCESS;
151 }

References _rset, moab::error(), ErrorCode, generate_hm(), level_dsequence, MB_CHK_ERR, MB_SUCCESS, meshdim, and nlevels.

Referenced by main().

◆ get_adjacencies()

ErrorCode moab::NestedRefine::get_adjacencies ( const EntityHandle  source_entity,
const unsigned int  target_dimension,
std::vector< EntityHandle > &  target_entities 
)

Get the adjacencies associated with an entity.

Given an entity of dimension d, gather all the adjacent D dimensional entities where D >, = , < d .

Parameters
source_entityEntityHandle to which adjacent entities have to be found.
target_dimensionInt Dimension of the desired adjacent entities.
target_entitiesVector in which the adjacent EntityHandle are returned.

Definition at line 223 of file NestedRefine.cpp.

227 {
229  error = ahf->get_adjacencies( source_entity, target_dimension, target_entities );MB_CHK_ERR( error );
230 
231  return MB_SUCCESS;
232 }

References ahf, moab::error(), ErrorCode, moab::HalfFacetRep::get_adjacencies(), MB_CHK_ERR, and MB_SUCCESS.

◆ get_connectivity()

ErrorCode moab::NestedRefine::get_connectivity ( EntityHandle  ent,
int  level,
std::vector< EntityHandle > &  conn 
)

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.

Parameters
entEntityHandle of the entity
levelInteger level of the entity for which connectivity is requested
connstd::vector returning the connectivity of the entity

Definition at line 153 of file NestedRefine.cpp.

154 {
156  EntityType type = mbImpl->type_from_handle( ent );
157  EntityHandle start_ent;
158  if( !conn.empty() ) conn.clear();
159  if( level > 0 )
160  {
161  if( type == MBEDGE )
162  {
163  conn.reserve( 2 );
164  start_ent = level_mesh[level - 1].start_edge;
165  EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
166  conn.push_back( level_mesh[level - 1].edge_conn[2 * offset] );
167  conn.push_back( level_mesh[level - 1].edge_conn[2 * offset + 1] );
168  }
169  else if( type == MBTRI || type == MBQUAD )
170  {
171  int num_corners = ahf->lConnMap2D[type - 2].num_verts_in_face;
172  conn.reserve( num_corners );
173  start_ent = level_mesh[level - 1].start_face;
174  EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
175 
176  for( int i = 0; i < num_corners; i++ )
177  conn.push_back( level_mesh[level - 1].face_conn[num_corners * offset + i] );
178  }
179  else if( type == MBTET || type == MBHEX )
180  {
181  int index = ahf->get_index_in_lmap( *_incells.begin() );
182  int num_corners = ahf->lConnMap3D[index].num_verts_in_cell;
183  conn.reserve( num_corners );
184  start_ent = level_mesh[level - 1].start_cell;
185  EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
186  for( int i = 0; i < num_corners; i++ )
187  conn.push_back( level_mesh[level - 1].cell_conn[num_corners * offset + i] );
188  }
189  else
190  MB_SET_ERR( MB_FAILURE, "Requesting connectivity for an unsupported entity type" );
191  }
192  else
193  {
194  error = mbImpl->get_connectivity( &ent, 1, conn );MB_CHK_ERR( error );
195  }
196 
197  return MB_SUCCESS;
198 }

References _incells, ahf, moab::Range::begin(), moab::NestedRefine::level_memory::cell_conn, moab::error(), ErrorCode, moab::NestedRefine::level_memory::face_conn, moab::Core::get_connectivity(), moab::HalfFacetRep::get_index_in_lmap(), moab::ID_FROM_HANDLE(), moab::HalfFacetRep::lConnMap2D, moab::HalfFacetRep::lConnMap3D, level_mesh, MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, MBEDGE, MBHEX, mbImpl, MBQUAD, MBTET, MBTRI, moab::HalfFacetRep::LocalMaps3D::num_verts_in_cell, moab::HalfFacetRep::LocalMaps2D::num_verts_in_face, moab::NestedRefine::level_memory::start_cell, moab::NestedRefine::level_memory::start_edge, moab::NestedRefine::level_memory::start_face, and moab::Core::type_from_handle().

Referenced by construct_hm_1D(), construct_hm_2D(), get_local_vid(), reorder_indices(), subdivide_cells(), subdivide_tets(), update_global_ahf_1D_sub(), update_global_ahf_2D(), and update_global_ahf_2D_sub().

◆ get_coordinates()

ErrorCode moab::NestedRefine::get_coordinates ( EntityHandle verts,
int  num_verts,
int  level,
double *  coords 
)

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.

Parameters
vertsstd::vector of the entity handles of the vertices
num_vertsThe number of vertices
levelInteger level of the entity for which connectivity is requested
coordsdouble pointer returning the coordinates of the vertices

Definition at line 200 of file NestedRefine.cpp.

201 {
202  if( level > 0 )
203  {
204  EntityID vstart = ID_FROM_HANDLE( level_mesh[level - 1].start_vertex );
205  for( int i = 0; i < num_verts; i++ )
206  {
207  const EntityHandle& vid = verts[i];
208  EntityID offset = ID_FROM_HANDLE( vid ) - vstart;
209  coords[3 * i] = level_mesh[level - 1].coordinates[0][offset];
210  coords[3 * i + 1] = level_mesh[level - 1].coordinates[1][offset];
211  coords[3 * i + 2] = level_mesh[level - 1].coordinates[2][offset];
212  }
213  }
214  else
215  {
217  error = mbImpl->get_coords( verts, num_verts, coords );MB_CHK_ERR( error );
218  }
219 
220  return MB_SUCCESS;
221 }

References moab::NestedRefine::level_memory::coordinates, moab::error(), ErrorCode, moab::Core::get_coords(), moab::ID_FROM_HANDLE(), level_mesh, MB_CHK_ERR, MB_SUCCESS, and mbImpl.

Referenced by construct_hm_2D(), subdivide_cells(), and subdivide_tets().

◆ get_index_from_degree()

◆ get_lid_inci_child()

ErrorCode moab::NestedRefine::get_lid_inci_child ( EntityType  type,
int  deg,
int  lfid,
int  leid,
std::vector< int > &  child_ids,
std::vector< int > &  child_lvids 
)
protected

Definition at line 3985 of file NestedRefine.cpp.

3991 {
3992  int index = ahf->get_index_in_lmap( *_incells.begin() );
3993  int d = get_index_from_degree( deg );
3994 
3995  // int lv0 = ahf->lConnMap3D[index].e2v[leid][0];
3996  // int lv1 = ahf->lConnMap3D[index].e2v[leid][1];
3997  int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
3998 
3999  int nv = refTemplates[type - 1][d].nv_edge;
4000  int nch = refTemplates[type - 1][d].ents_on_pent[lfid][0];
4001 
4002  for( int i = 0; i < nch; i++ )
4003  {
4004  int id = refTemplates[type - 1][d].ents_on_pent[lfid][i + 1] - 1;
4005  for( int j = 0; j < nvpc; j++ )
4006  {
4007  int lv = refTemplates[type - 1][d].ents_conn[id][j];
4008  for( int k = 0; k < nv; k++ )
4009  {
4010  if( lv == refTemplates[type - 1][d].vert_on_edges[leid][k] )
4011  {
4012  child_ids.push_back( id );
4013  child_lvids.push_back( j );
4014  }
4015  }
4016  }
4017  }
4018 
4019  return MB_SUCCESS;
4020 }

References _incells, ahf, moab::Range::begin(), moab::NestedRefine::refPatterns::ents_conn, moab::NestedRefine::refPatterns::ents_on_pent, get_index_from_degree(), moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::lConnMap3D, MB_SUCCESS, moab::HalfFacetRep::LocalMaps3D::num_verts_in_cell, moab::NestedRefine::refPatterns::nv_edge, and refTemplates.

◆ get_local_vid()

int moab::NestedRefine::get_local_vid ( EntityHandle  vid,
EntityHandle  ent,
int  level 
)
protected

Definition at line 4612 of file NestedRefine.cpp.

4613 {
4614  ErrorCode error;
4615  // Given a vertex, find its local id in the given entity
4616  std::vector< EntityHandle > conn;
4617 
4618  error = get_connectivity( ent, level + 1, conn );
4619  if( error != MB_SUCCESS ) MB_SET_ERR( MB_FAILURE, "Error in getting connectivity of the requested entity" );
4620 
4621  int lid = -1;
4622  for( int i = 0; i < (int)conn.size(); i++ )
4623  {
4624  if( conn[i] == vid )
4625  {
4626  lid = i;
4627  break;
4628  }
4629  }
4630  if( lid < 0 ) MB_SET_ERR( MB_FAILURE, "Error in getting local vertex id in the given entity" );
4631  return lid;
4632 }

References moab::error(), ErrorCode, get_connectivity(), MB_SET_ERR, and MB_SUCCESS.

Referenced by update_global_ahf_1D(), update_global_ahf_1D_sub(), update_global_ahf_2D(), update_global_ahf_2D_sub(), and update_global_ahf_3D().

◆ get_octahedron_corner_coords()

ErrorCode moab::NestedRefine::get_octahedron_corner_coords ( int  cur_level,
int  deg,
EntityHandle vbuffer,
double *  ocoords 
)
protected

Definition at line 4544 of file NestedRefine.cpp.

4545 {
4546  int lid[6] = { 0, 0, 0, 0, 0, 0 };
4547 
4548  if( deg == 2 )
4549  {
4550  lid[0] = 5;
4551  lid[1] = 8;
4552  lid[2] = 9;
4553  lid[3] = 6;
4554  lid[4] = 4;
4555  lid[5] = 7;
4556  }
4557  else if( deg == 3 )
4558  {
4559  lid[0] = 19;
4560  lid[1] = 16;
4561  lid[2] = 18;
4562  lid[3] = 9;
4563  lid[4] = 4;
4564  lid[5] = 10;
4565  }
4566 
4567  EntityHandle vstart = level_mesh[cur_level].start_vertex;
4568 
4569  for( int i = 0; i < 6; i++ )
4570  {
4571  EntityHandle vid = vbuffer[lid[i]];
4572  ocoords[3 * i] = level_mesh[cur_level].coordinates[0][vid - vstart];
4573  ocoords[3 * i + 1] = level_mesh[cur_level].coordinates[1][vid - vstart];
4574  ocoords[3 * i + 2] = level_mesh[cur_level].coordinates[2][vid - vstart];
4575  }
4576 
4577  return MB_SUCCESS;
4578 }

References moab::NestedRefine::level_memory::coordinates, level_mesh, MB_SUCCESS, and moab::NestedRefine::level_memory::start_vertex.

Referenced by find_shortest_diagonal_octahedron().

◆ get_vertex_duplicates()

ErrorCode moab::NestedRefine::get_vertex_duplicates ( EntityHandle  vertex,
int  level,
EntityHandle dupvertex 
)

Definition at line 425 of file NestedRefine.cpp.

426 {
427  if( ( vertex - *_inverts.begin() ) > _inverts.size() )
428  MB_SET_ERR( MB_FAILURE, "Requesting duplicates for non-coarse vertices" );
429 
430  dupvertex = level_mesh[level - 1].start_vertex + ( vertex - *_inverts.begin() );
431 
432  return MB_SUCCESS;
433 }

References _inverts, moab::Range::begin(), level_mesh, MB_SET_ERR, MB_SUCCESS, moab::Range::size(), and moab::NestedRefine::level_memory::start_vertex.

◆ initialize()

ErrorCode moab::NestedRefine::initialize ( )

Definition at line 53 of file NestedRefine.cpp.

54 {
56 
57  tm = new CpuTimer();
58  if( !tm ) return MB_MEMORY_ALLOCATION_FAILED;
59 
60 #ifdef MOAB_HAVE_AHF
61  ahf = mbImpl->a_half_facet_rep();
62 #else
63  ahf = new HalfFacetRep( mbImpl, pcomm, _rset, true );
64  if( !ahf ) return MB_MEMORY_ALLOCATION_FAILED;
65 #endif
66 
67  // Check for mixed entity type
68  bool chk_mixed = ahf->check_mixed_entity_type();
69  if( chk_mixed ) MB_SET_ERR( MB_NOT_IMPLEMENTED, "Encountered a mesh with mixed entity types" );
70 
73 
74  // Check for supported entity type
75  if( !_incells.empty() )
76  {
77  EntityType type = mbImpl->type_from_handle( _incells[0] );
78  if( type != MBTET && type != MBHEX )
79  MB_SET_ERR( MB_FAILURE, "Not supported 3D entity types: MBPRISM, MBPYRAMID, MBKNIFE, MBPOLYHEDRON" );
80 
81  meshdim = 3;
82  elementype = type;
83  }
84  else if( !_infaces.empty() )
85  {
86  EntityType type = mbImpl->type_from_handle( _infaces[0] );
87  if( type == MBPOLYGON ) MB_SET_ERR( MB_FAILURE, "Not supported 2D entity type: POLYGON" );
88 
89  meshdim = 2;
90  elementype = type;
91  }
92  else if( !_inedges.empty() )
93  {
94  meshdim = 1;
96  }
97  else
98  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Encountered a mixed-dimensional or invalid mesh" );
99 
100  // Initialize std::map to get indices of degrees.
101  deg_index[2] = 0;
102  deg_index[3] = 1;
103  deg_index[5] = 2;
104 
105  // Set ghost flag to false
106  hasghost = false;
107  return MB_SUCCESS;
108 }

References _incells, _inedges, _infaces, _inverts, _rset, ahf, moab::HalfFacetRep::check_mixed_entity_type(), deg_index, elementype, moab::Range::empty(), moab::error(), ErrorCode, moab::HalfFacetRep::get_entity_ranges(), hasghost, moab::HalfFacetRep::initialize(), MB_CHK_ERR, MB_MEMORY_ALLOCATION_FAILED, MB_NOT_IMPLEMENTED, MB_SET_ERR, MB_SUCCESS, MBEDGE, MBHEX, mbImpl, MBPOLYGON, MBTET, meshdim, pcomm, tm, and moab::Core::type_from_handle().

Referenced by NestedRefine().

◆ is_cell_on_boundary()

bool moab::NestedRefine::is_cell_on_boundary ( const EntityHandle entity)
protected

Definition at line 4138 of file NestedRefine.cpp.

4139 {
4140  if( meshdim != 3 )
4141  MB_SET_ERR( MB_FAILURE, "Requesting boundary information for a cell entity type on a curve or surface mesh" );
4142 
4143  bool is_border = false;
4144  int index = ahf->get_index_in_lmap( *_incells.begin() );
4145  int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
4146  EntityHandle sibents[6];
4147  int siblids[6];
4148 
4149  ErrorCode error = ahf->get_sibling_map( elementype, entity, &sibents[0], &siblids[0], nfpc );MB_CHK_ERR( error );
4150 
4151  for( int i = 0; i < nfpc; i++ )
4152  {
4153  if( sibents[i] == 0 )
4154  {
4155  is_border = true;
4156  break;
4157  }
4158  }
4159  return is_border;
4160 }

References _incells, ahf, moab::Range::begin(), elementype, moab::error(), ErrorCode, moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::get_sibling_map(), moab::HalfFacetRep::lConnMap3D, MB_CHK_ERR, MB_SET_ERR, meshdim, and moab::HalfFacetRep::LocalMaps3D::num_faces_in_cell.

Referenced by is_entity_on_boundary().

◆ is_edge_on_boundary()

bool moab::NestedRefine::is_edge_on_boundary ( const EntityHandle entity)
protected

Definition at line 4053 of file NestedRefine.cpp.

4054 {
4055  ErrorCode error;
4056  bool is_border = false;
4057  if( meshdim == 1 ) // The edge has a vertex on the boundary in the curve mesh
4058  {
4059  EntityHandle sibents[2];
4060  int siblids[2];
4061  error = ahf->get_sibling_map( MBEDGE, entity, &sibents[0], &siblids[0], 2 );MB_CHK_ERR( error );
4062  for( int i = 0; i < 2; i++ )
4063  {
4064  if( sibents[i] == 0 )
4065  {
4066  is_border = true;
4067  break;
4068  }
4069  }
4070  }
4071  else if( meshdim == 2 ) // The edge is on the boundary of the 2d mesh
4072  {
4073  std::vector< EntityHandle > adjents;
4074  error = ahf->get_up_adjacencies_2d( entity, adjents );MB_CHK_ERR( error );
4075  if( adjents.size() == 1 ) is_border = true;
4076  }
4077  else if( meshdim == 3 ) // The edge is part of a face on the boundary of the 3d mesh
4078  {
4079  std::vector< EntityHandle > adjents;
4080  std::vector< int > leids;
4081  error = ahf->get_up_adjacencies_edg_3d( entity, adjents, &leids );MB_CHK_ERR( error );
4082  assert( !adjents.empty() );
4083 
4084  int index = ahf->get_index_in_lmap( adjents[0] );
4085  int nhf = ahf->lConnMap3D[index].num_faces_in_cell;
4086 
4087  for( int i = 0; i < (int)adjents.size(); i++ )
4088  {
4089  EntityHandle sibents[6];
4090  int siblids[6];
4091  error = ahf->get_sibling_map( elementype, adjents[0], &sibents[0], &siblids[0], nhf );MB_CHK_ERR( error );
4092  for( int k = 0; k < 2; k++ )
4093  {
4094  int hf = ahf->lConnMap3D[index].e2hf[leids[0]][k];
4095  if( sibents[hf] == 0 )
4096  {
4097  is_border = true;
4098  break;
4099  }
4100  }
4101  }
4102  }
4103  return is_border;
4104 }

References ahf, moab::HalfFacetRep::LocalMaps3D::e2hf, elementype, moab::error(), ErrorCode, moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::get_sibling_map(), moab::HalfFacetRep::get_up_adjacencies_2d(), moab::HalfFacetRep::get_up_adjacencies_edg_3d(), moab::HalfFacetRep::lConnMap3D, MB_CHK_ERR, MBEDGE, meshdim, and moab::HalfFacetRep::LocalMaps3D::num_faces_in_cell.

Referenced by is_entity_on_boundary().

◆ is_entity_on_boundary()

bool moab::NestedRefine::is_entity_on_boundary ( const EntityHandle entity)

Given an entity at a certain level, it returns a boolean value true if it lies on the domain boundary.

Parameters
entity

Definition at line 435 of file NestedRefine.cpp.

436 {
437  bool is_border = false;
438  EntityType type = mbImpl->type_from_handle( entity );
439 
440  if( type == MBVERTEX )
441  is_border = is_vertex_on_boundary( entity );
442  else if( type == MBEDGE )
443  is_border = is_edge_on_boundary( entity );
444  else if( type == MBTRI || type == MBQUAD )
445  is_border = is_face_on_boundary( entity );
446  else if( type == MBTET || type == MBHEX )
447  is_border = is_cell_on_boundary( entity );
448  else
449  MB_SET_ERR( MB_FAILURE, "Requesting boundary information for unsupported entity type" );
450 
451  return is_border;
452 }

References is_cell_on_boundary(), is_edge_on_boundary(), is_face_on_boundary(), is_vertex_on_boundary(), MB_SET_ERR, MBEDGE, MBHEX, mbImpl, MBQUAD, MBTET, MBTRI, MBVERTEX, and moab::Core::type_from_handle().

◆ is_face_on_boundary()

bool moab::NestedRefine::is_face_on_boundary ( const EntityHandle entity)
protected

Definition at line 4106 of file NestedRefine.cpp.

4107 {
4108  ErrorCode error;
4109  bool is_border = false;
4110 
4111  if( meshdim == 1 )
4112  MB_SET_ERR( MB_FAILURE, "Requesting boundary information for a face entity type on a curve mesh" );
4113  else if( meshdim == 2 ) // The face has a local edge on the boundary of the 2d mesh
4114  {
4115  EntityHandle sibents[4];
4116  int siblids[4];
4117  int nepf = ahf->lConnMap2D[elementype - 2].num_verts_in_face;
4118  error = ahf->get_sibling_map( elementype, entity, &sibents[0], &siblids[0], nepf );MB_CHK_ERR( error );
4119 
4120  for( int i = 0; i < nepf; i++ )
4121  {
4122  if( sibents[i] == 0 )
4123  {
4124  is_border = true;
4125  break;
4126  }
4127  }
4128  }
4129  else if( meshdim == 3 ) // The face lies on the boundary of the 3d mesh
4130  {
4131  std::vector< EntityHandle > adjents;
4132  error = ahf->get_up_adjacencies_face_3d( entity, adjents );MB_CHK_ERR( error );
4133  if( adjents.size() == 1 ) is_border = true;
4134  }
4135  return is_border;
4136 }

References ahf, elementype, moab::error(), ErrorCode, moab::HalfFacetRep::get_sibling_map(), moab::HalfFacetRep::get_up_adjacencies_face_3d(), moab::HalfFacetRep::lConnMap2D, MB_CHK_ERR, MB_SET_ERR, meshdim, and moab::HalfFacetRep::LocalMaps2D::num_verts_in_face.

Referenced by is_entity_on_boundary().

◆ is_vertex_on_boundary()

bool moab::NestedRefine::is_vertex_on_boundary ( const EntityHandle entity)
protected

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.

Parameters
entity

Definition at line 4026 of file NestedRefine.cpp.

4027 {
4028  ErrorCode error;
4029  EntityHandle sibents[27];
4030  int siblids[27];
4031  std::vector< EntityHandle > ent;
4032  std::vector< int > lid;
4033 
4034  int nhf;
4035  if( elementype == MBEDGE )
4036  nhf = 2;
4037  else if( ( elementype == MBTRI ) || ( elementype == MBQUAD ) )
4039  else if( ( elementype == MBTET ) || ( elementype == MBHEX ) )
4040  {
4041  int idx = ahf->get_index_in_lmap( *_incells.begin() );
4042  nhf = ahf->lConnMap3D[idx].num_faces_in_cell;
4043  }
4044  else
4045  MB_SET_ERR( MB_FAILURE, "Requesting vertex boundary information for an unsupported entity type" );
4046 
4048  error = ahf->get_sibling_map( elementype, ent[0], &sibents[0], &siblids[0], nhf );MB_CHK_ERR( error );
4049 
4050  return ( sibents[lid[0]] == 0 );
4051 }

References _incells, ahf, moab::Range::begin(), elementype, moab::error(), ErrorCode, moab::HalfFacetRep::get_incident_map(), moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::get_sibling_map(), moab::HalfFacetRep::lConnMap2D, moab::HalfFacetRep::lConnMap3D, MB_CHK_ERR, MB_SET_ERR, MBEDGE, MBHEX, MBQUAD, MBTET, MBTRI, moab::HalfFacetRep::LocalMaps3D::num_faces_in_cell, and moab::HalfFacetRep::LocalMaps2D::num_verts_in_face.

Referenced by is_entity_on_boundary().

◆ parent_to_child()

ErrorCode moab::NestedRefine::parent_to_child ( EntityHandle  parent,
int  parent_level,
int  child_level,
std::vector< EntityHandle > &  children 
)

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.

Parameters
parentEntityHandle of the entity whose children in subsequent level is requested
parent_levelMesh level where the parent exists
child_levelMesh level from which its children are requested
childrenVector containing all childrens from the requested child_level

Definition at line 285 of file NestedRefine.cpp.

289 {
290  assert( ( child_level > 0 ) && ( child_level > parent_level ) );
291  EntityType type = mbImpl->type_from_handle( parent );
292  assert( type != MBVERTEX );
293 
294  int parent_index;
295  if( type == MBEDGE )
296  {
297  if( parent_level > 0 )
298  parent_index = parent - level_mesh[parent_level - 1].start_edge;
299  else
300  parent_index = _inedges.index( parent );
301  }
302  else if( type == MBTRI || type == MBQUAD )
303  {
304  if( parent_level > 0 )
305  parent_index = parent - level_mesh[parent_level - 1].start_face;
306  else
307  parent_index = _infaces.index( parent );
308  }
309  else if( type == MBTET || type == MBHEX )
310  {
311  if( parent_level > 0 )
312  parent_index = parent - level_mesh[parent_level - 1].start_cell;
313  else
314  parent_index = _incells.index( parent );
315  }
316  else
317  MB_SET_ERR( MB_FAILURE, "Requesting children for unsupported entity type" );
318 
319  int start, end;
320  start = end = parent_index;
321  for( int i = parent_level; i < child_level; i++ )
322  {
324  int nch = refTemplates[type - 1][d].total_new_ents;
325  start = start * nch;
326  end = end * nch + nch - 1;
327  }
328 
329  int num_child = end - start;
330  children.reserve( num_child );
331 
332  for( int i = start; i <= end; i++ )
333  {
335  if( type == MBEDGE )
336  child = level_mesh[child_level - 1].start_edge + i;
337  else if( type == MBTRI || type == MBQUAD )
338  child = level_mesh[child_level - 1].start_face + i;
339  else if( type == MBTET || type == MBHEX )
340  child = level_mesh[child_level - 1].start_cell + i;
341 
342  children.push_back( child );
343  }
344 
345  return MB_SUCCESS;
346 }

References _incells, _inedges, _infaces, child, children, get_index_from_degree(), moab::Range::index(), level_dsequence, level_mesh, MB_SET_ERR, MB_SUCCESS, MBEDGE, MBHEX, mbImpl, MBQUAD, MBTET, MBTRI, MBVERTEX, refTemplates, moab::NestedRefine::level_memory::start_cell, moab::NestedRefine::level_memory::start_edge, moab::NestedRefine::level_memory::start_face, moab::NestedRefine::refPatterns::total_new_ents, and moab::Core::type_from_handle().

Referenced by update_special_tags(), and vertex_to_entities_down().

◆ print_maps_1D()

ErrorCode moab::NestedRefine::print_maps_1D ( int  level)
protected

◆ print_maps_2D()

ErrorCode moab::NestedRefine::print_maps_2D ( int  level,
EntityType  type 
)
protected

◆ print_maps_3D()

ErrorCode moab::NestedRefine::print_maps_3D ( int  level,
EntityType  type 
)
protected

◆ reorder_indices() [1/4]

ErrorCode moab::NestedRefine::reorder_indices ( EntityHandle face1_conn,
EntityHandle face2_conn,
int  nvF,
int *  conn_map,
int &  comb,
int *  orient = NULL 
)
protected

Definition at line 4477 of file NestedRefine.cpp.

4483 {
4484  // Given connectivities of two faces and a degree, get the permuted indices of the children
4485  // faces w.r.t first face.
4486 
4487  // Step 1: First find the combination
4488  int nco = permutation[nvF - 3].num_comb;
4489  int c = 0;
4490  for( int i = 0; i < nco; i++ )
4491  {
4492  int count = 0;
4493  for( int j = 0; j < nvF; j++ )
4494  {
4495  int id = permutation[nvF - 3].comb[i][j];
4496  if( face1_conn[j] == face2_conn[id] ) count += 1;
4497  }
4498 
4499  if( count == nvF )
4500  {
4501  c = i;
4502  break;
4503  }
4504  }
4505 
4506  if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
4507 
4508  comb = c;
4509 
4510  if( orient ) orient[0] = permutation[nvF - 3].orient[c];
4511 
4512  for( int j = 0; j < nvF; j++ )
4513  {
4514  conn_map[j] = permutation[nvF - 3].comb[c][j];
4515  }
4516 
4517  return MB_SUCCESS;
4518 }

References moab::NestedRefine::pmat::comb, MB_SET_ERR, MB_SUCCESS, moab::NestedRefine::pmat::num_comb, moab::NestedRefine::pmat::orient, and permutation.

◆ reorder_indices() [2/4]

ErrorCode moab::NestedRefine::reorder_indices ( int  cur_level,
int  deg,
EntityHandle  cell,
int  lfid,
EntityHandle  sib_cell,
int  sib_lfid,
int  index,
int *  id_sib 
)
protected

Definition at line 4322 of file NestedRefine.cpp.

4330 {
4331  // Reorders the indices of either vertices or children cell local ids to match with order of the
4332  // given cell and a local face. index = 0 : vertices,
4333  // = 1 : face
4334 
4335  assert( deg == 2 || deg == 3 );
4336 
4337  ErrorCode error;
4338  int idx = ahf->get_index_in_lmap( *_incells.begin() );
4339  int nvF = ahf->lConnMap3D[idx].hf2v_num[lfid];
4340  int nco = permutation[nvF - 3].num_comb;
4341 
4342  if( !index && ( ( nvF == 3 && deg == 3 ) || ( nvF == 4 && deg == 2 ) ) )
4343  {
4344  id_sib[0] = 1;
4345  }
4346  else
4347  {
4348  // Get connectivity of the cell and its sibling cell
4349  std::vector< EntityHandle > conn, sib_conn;
4350  error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
4351 
4352  error = get_connectivity( sib_cell, cur_level, sib_conn );MB_CHK_ERR( error );
4353 
4354  // Get the connectivity of the local face in the cell and its sibling
4355  std::vector< EntityHandle > lface( nvF );
4356  std::vector< EntityHandle > lface_sib( nvF );
4357  for( int i = 0; i < nvF; i++ )
4358  {
4359  int id = ahf->lConnMap3D[idx].hf2v[lfid][i];
4360  lface[i] = conn[id];
4361 
4362  id = ahf->lConnMap3D[idx].hf2v[sib_lfid][i];
4363  lface_sib[i] = sib_conn[id];
4364  }
4365 
4366  // Find the combination
4367  int c = 0;
4368  for( int i = 0; i < nco; i++ )
4369  {
4370  int count = 0;
4371  for( int j = 0; j < nvF; j++ )
4372  {
4373  int id = permutation[nvF - 3].comb[i][j];
4374  if( lface[j] == lface_sib[id] ) count += 1;
4375  }
4376 
4377  if( count == nvF )
4378  {
4379  c = i;
4380  break;
4381  }
4382  }
4383 
4384  if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
4385 
4386  // Get the ordered indices
4387  if( ( ( !index ) && ( nvF == 4 ) && ( deg == 3 ) ) || ( deg == 2 ) )
4388  {
4389  for( int i = 0; i < 4; i++ )
4390  id_sib[i] = permutation[nvF - 3].porder2[c][i];
4391  }
4392  else
4393  {
4394  for( int i = 0; i < 9; i++ )
4395  id_sib[i] = permutation[nvF - 3].porder3[c][i];
4396  }
4397  }
4398 
4399  return MB_SUCCESS;
4400 }

References _incells, ahf, moab::Range::begin(), moab::NestedRefine::pmat::comb, moab::error(), ErrorCode, get_connectivity(), moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::LocalMaps3D::hf2v, moab::HalfFacetRep::LocalMaps3D::hf2v_num, moab::HalfFacetRep::lConnMap3D, MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, moab::NestedRefine::pmat::num_comb, and permutation.

Referenced by construct_hm_2D(), update_global_ahf_3D(), and update_tracking_verts().

◆ reorder_indices() [3/4]

ErrorCode moab::NestedRefine::reorder_indices ( int  deg,
EntityHandle face1_conn,
EntityHandle face2_conn,
int  nvF,
std::vector< int > &  lemap,
std::vector< int > &  vidx,
int *  leorient = NULL 
)
protected

Definition at line 4402 of file NestedRefine.cpp.

4409 {
4410  // Given the connectivities of two faces, get the permuted indices w.r.t first face.
4411  // Step 1: First find the orientation
4412  int nco = permutation[nvF - 3].num_comb;
4413  int c = 0;
4414  for( int i = 0; i < nco; i++ )
4415  {
4416  int count = 0;
4417  for( int j = 0; j < nvF; j++ )
4418  {
4419  int id = permutation[nvF - 3].comb[i][j];
4420  if( face1_conn[j] == face2_conn[id] ) count += 1;
4421  }
4422 
4423  if( count == nvF )
4424  {
4425  c = i;
4426  break;
4427  }
4428  }
4429 
4430  if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
4431 
4432  // Add the corresponding local edges
4433  lemap.reserve( nvF );
4434  for( int i = 0; i < nvF; i++ )
4435  {
4436  lemap.push_back( permutation[nvF - 3].lemap[c][i] );
4437  }
4438  if( leorient ) leorient[0] = permutation[nvF - 3].orient[c];
4439 
4440  if( nvF == 3 && deg == 2 ) return MB_SUCCESS;
4441 
4442  if( ( nvF == 3 && deg == 3 ) || ( nvF == 4 && deg == 2 ) )
4443  {
4444  vidx.push_back( 1 );
4445  }
4446  else if( nvF == 4 && deg == 3 )
4447  {
4448  for( int i = 0; i < 4; i++ )
4449  vidx.push_back( permutation[nvF - 3].porder2[c][i] );
4450  }
4451 
4452  return MB_SUCCESS;
4453 }

References moab::NestedRefine::pmat::comb, MB_SET_ERR, MB_SUCCESS, moab::NestedRefine::pmat::num_comb, moab::NestedRefine::pmat::orient, permutation, and moab::NestedRefine::pmat::porder2.

◆ reorder_indices() [4/4]

ErrorCode moab::NestedRefine::reorder_indices ( int  deg,
int  nvF,
int  comb,
int *  childfid_map 
)
protected

Definition at line 4455 of file NestedRefine.cpp.

4456 {
4457  // Given connectivities of two faces and a degree, get the permuted indices of the children
4458  // faces w.r.t first face.
4459 
4460  assert( deg == 2 || deg == 3 );
4461 
4462  // Get the ordered indices
4463  if( deg == 2 )
4464  {
4465  for( int i = 0; i < 4; i++ )
4466  childfid_map[i] = permutation[nvF - 3].porder2[comb][i];
4467  }
4468  else
4469  {
4470  for( int i = 0; i < 9; i++ )
4471  childfid_map[i] = permutation[nvF - 3].porder3[comb][i];
4472  }
4473 
4474  return MB_SUCCESS;
4475 }

References MB_SUCCESS, and permutation.

◆ subdivide_cells()

ErrorCode moab::NestedRefine::subdivide_cells ( EntityType  type,
int  cur_level,
int  deg 
)
protected

Definition at line 1421 of file NestedRefine.cpp.

1422 {
1423  ErrorCode error;
1424  int nverts_prev, nents_prev;
1425  if( cur_level )
1426  {
1427  nverts_prev = level_mesh[cur_level - 1].num_verts;
1428  nents_prev = level_mesh[cur_level - 1].num_cells;
1429  }
1430  else
1431  {
1432  nverts_prev = _inverts.size();
1433  nents_prev = _incells.size();
1434  }
1435 
1436  // Create some book-keeping arrays over the parent mesh to avoid introducing duplicate vertices
1437  int cindex = type - 1;
1438  int d = get_index_from_degree( deg );
1439  int ne = refTemplates[cindex][d].nv_edge;
1440  int nvf = refTemplates[cindex][d].nv_face;
1441  int nvtotal = refTemplates[cindex][d].total_new_verts;
1442 
1443  int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
1444  int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
1445  int nepc = ahf->lConnMap3D[index].num_edges_in_cell;
1446  int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
1447 
1448  int vtotal = nvpc + nvtotal;
1449  std::vector< EntityHandle > vbuffer( vtotal );
1450 
1451  std::vector< EntityHandle > trackvertsC_edg( nepc * ne * nents_prev, 0 );
1452  std::vector< EntityHandle > trackvertsC_face( nfpc * nvf * nents_prev, 0 );
1453 
1454  int cur_nverts = level_mesh[cur_level].num_verts;
1455  std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
1456 
1457  int count_nverts = nverts_prev;
1458  int count_ents = 0;
1459  std::vector< EntityHandle > conn, cur_conn;
1460 
1461  // Step 1: Create the subentities via refinement of the previous mesh
1462  for( int cid = 0; cid < nents_prev; cid++ )
1463  {
1464  conn.clear();
1465  cur_conn.clear();
1466  for( int i = 0; i < vtotal; i++ )
1467  vbuffer[i] = 0;
1468 
1469  // EntityHandle of the working cell
1470  EntityHandle cell;
1471  if( cur_level )
1472  cell = level_mesh[cur_level - 1].start_cell + cid;
1473  else
1474  cell = _incells[cid];
1475 
1476  error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
1477 
1478  // Step 1: Add vertices from the current level for the working face that will be used for
1479  // subdivision.
1480  // Add the corners to vbuffer
1481  for( int i = 0; i < (int)conn.size(); i++ )
1482  {
1483  if( cur_level )
1484  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
1485  else
1486  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
1487 
1488  cur_conn.push_back( vbuffer[i] );
1489  }
1490 
1491  // Gather vertices already added to tracking array due to refinement of the sibling cells
1492  for( int i = 0; i < nepc; i++ )
1493  {
1494  for( int j = 0; j < ne; j++ )
1495  {
1496  int idx = refTemplates[cindex][d].vert_on_edges[i][j];
1497  vbuffer[idx] = trackvertsC_edg[cid * nepc * ne + ne * i + j];
1498  }
1499  }
1500 
1501  // Add remaining new vertex handles
1502  for( int i = 0; i < nfpc; i++ )
1503  {
1504  for( int j = 0; j < nvf; j++ )
1505  {
1506  int idx = refTemplates[cindex][d].vert_on_faces[i][j];
1507  vbuffer[idx] = trackvertsC_face[cid * nfpc * nvf + nvf * i + j];
1508  }
1509  }
1510 
1511  // Add the remaining vertex handles to vbuffer for the current level for the working cell
1512  for( int i = 0; i < nvtotal; i++ )
1513  {
1514  if( !vbuffer[i + nvpc] )
1515  {
1516  vbuffer[i + nvpc] = level_mesh[cur_level].start_vertex + count_nverts;
1517  count_nverts += 1;
1518  }
1519  }
1520 
1521  // Step 2: Use the template to obtain the subentities. The coordinates and local ahf maps
1522  // are also constructed. Connectivity of the children
1523  int etotal = refTemplates[type - 1][d].total_new_ents;
1524  std::vector< EntityHandle > ent_buffer( etotal );
1525 
1526  for( int i = 0; i < etotal; i++ )
1527  {
1528  for( int k = 0; k < nvpc; k++ )
1529  {
1530  int idx = refTemplates[type - 1][d].ents_conn[i][k];
1531  level_mesh[cur_level].cell_conn[nvpc * count_ents + k] = vbuffer[idx];
1532  }
1533  ent_buffer[i] = level_mesh[cur_level].start_cell + count_ents;
1534  count_ents += 1;
1535  }
1536 
1537  // Step 3: Update local ahf maps
1538  error = update_local_ahf( deg, type, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
1539 
1540  // Step 4: Update tracking information
1541  error = update_tracking_verts( cell, cur_level, deg, trackvertsC_edg, trackvertsC_face, &vbuffer[0] );MB_CHK_ERR( error );
1542 
1543  // Step 5: Coordinates of the new vertices
1544  std::vector< double > corner_coords( nvpc * 3 );
1545  error = get_coordinates( &cur_conn[0], nvpc, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
1546 
1547  error = compute_coordinates( cur_level, deg, type, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
1548  nverts_prev );MB_CHK_ERR( error );
1549  }
1550 
1551  // error = ahf->print_tags(3);
1552 
1553  // Step 6: Update the global maps
1554  error = update_global_ahf( type, cur_level, deg );MB_CHK_ERR( error );
1555 
1556  // Step 7: If edges exists, refine them as well.
1557  if( level_mesh[cur_level].num_edges != 0 )
1558  {
1559  error = construct_hm_1D( cur_level, deg, type, trackvertsC_edg );MB_CHK_ERR( error );
1560  }
1561 
1562  // Step 8: If faces exists, refine them as well.
1563  if( !_infaces.empty() )
1564  {
1565  error = construct_hm_2D( cur_level, deg, type, trackvertsC_edg, trackvertsC_face );MB_CHK_ERR( error );
1566  }
1567 
1568  // error = ahf->print_tags(3);
1569 
1570  return MB_SUCCESS;
1571 }

References _incells, _infaces, _inverts, ahf, moab::Range::begin(), moab::NestedRefine::level_memory::cell_conn, compute_coordinates(), construct_hm_1D(), construct_hm_2D(), moab::Range::empty(), moab::NestedRefine::refPatterns::ents_conn, moab::error(), ErrorCode, get_connectivity(), get_coordinates(), get_index_from_degree(), moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::lConnMap3D, level_mesh, MB_CHK_ERR, MB_SUCCESS, moab::NestedRefine::level_memory::num_cells, moab::HalfFacetRep::LocalMaps3D::num_edges_in_cell, moab::HalfFacetRep::LocalMaps3D::num_faces_in_cell, moab::NestedRefine::level_memory::num_verts, moab::HalfFacetRep::LocalMaps3D::num_verts_in_cell, moab::NestedRefine::refPatterns::nv_edge, moab::NestedRefine::refPatterns::nv_face, refTemplates, moab::Range::size(), moab::NestedRefine::level_memory::start_cell, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, moab::NestedRefine::refPatterns::total_new_verts, update_global_ahf(), update_local_ahf(), update_tracking_verts(), moab::NestedRefine::refPatterns::vert_on_edges, and moab::NestedRefine::refPatterns::vert_on_faces.

Referenced by construct_hm_3D().

◆ subdivide_tets()

ErrorCode moab::NestedRefine::subdivide_tets ( int  cur_level,
int  deg 
)
protected

Definition at line 1573 of file NestedRefine.cpp.

1574 {
1575  ErrorCode error;
1576  int nverts_prev, nents_prev;
1577  if( cur_level )
1578  {
1579  nverts_prev = level_mesh[cur_level - 1].num_verts;
1580  nents_prev = level_mesh[cur_level - 1].num_cells;
1581  }
1582  else
1583  {
1584  nverts_prev = _inverts.size();
1585  nents_prev = _incells.size();
1586  }
1587 
1588  EntityType type = MBTET;
1589  int cindex = type - 1;
1590  int d = get_index_from_degree( deg );
1591  int ne = refTemplates[cindex][d].nv_edge;
1592  int nvf = refTemplates[cindex][d].nv_face;
1593  int nvtotal = refTemplates[cindex][d].total_new_verts;
1594 
1595  int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
1596  int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
1597  int nepc = ahf->lConnMap3D[index].num_edges_in_cell;
1598  int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
1599 
1600  // Create vertex buffer
1601  int vtotal = nvpc + nvtotal;
1602  std::vector< EntityHandle > vbuffer( vtotal );
1603 
1604  // Create book-keeping arrays over the parent mesh to avoid introducing duplicate vertices
1605  std::vector< EntityHandle > trackvertsC_edg( nepc * ne * nents_prev, 0 );
1606  std::vector< EntityHandle > trackvertsC_face( nfpc * nvf * nents_prev, 0 );
1607 
1608  int cur_nverts = level_mesh[cur_level].num_verts;
1609  std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
1610  std::vector< int > cell_patterns( nents_prev, 0 );
1611 
1612  int count_nverts = nverts_prev;
1613  int count_ents = 0;
1614  std::vector< EntityHandle > conn, cur_conn;
1615 
1616  // Step 1: Create the subentities via refinement of the previous mesh
1617  for( int cid = 0; cid < nents_prev; cid++ )
1618  {
1619  conn.clear();
1620  cur_conn.clear();
1621  for( int i = 0; i < vtotal; i++ )
1622  vbuffer[i] = 0;
1623 
1624  // EntityHandle of the working cell
1625  EntityHandle cell;
1626  if( cur_level )
1627  cell = level_mesh[cur_level - 1].start_cell + cid;
1628  else
1629  cell = _incells[cid];
1630 
1631  error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
1632 
1633  // Step 1: Add vertices from the current level for the working face that will be used for
1634  // subdivision.
1635  // Add the corners to vbuffer
1636  for( int i = 0; i < (int)conn.size(); i++ )
1637  {
1638  if( cur_level )
1639  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
1640  else
1641  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
1642 
1643  cur_conn.push_back( vbuffer[i] );
1644  }
1645 
1646  // Gather vertices already added to tracking array due to refinement of the sibling cells
1647  for( int i = 0; i < nepc; i++ )
1648  {
1649  for( int j = 0; j < ne; j++ )
1650  {
1651  int idx = refTemplates[cindex][d].vert_on_edges[i][j];
1652  vbuffer[idx] = trackvertsC_edg[cid * nepc * ne + ne * i + j];
1653  }
1654  }
1655 
1656  // Add remaining new vertex handles
1657  for( int i = 0; i < nfpc; i++ )
1658  {
1659  for( int j = 0; j < nvf; j++ )
1660  {
1661  int idx = refTemplates[cindex][d].vert_on_faces[i][j];
1662  vbuffer[idx] = trackvertsC_face[cid * nfpc * nvf + nvf * i + j];
1663  }
1664  }
1665 
1666  // Add the remaining vertex handles to vbuffer for the current level for the working cell
1667  for( int i = 0; i < nvtotal; i++ )
1668  {
1669  if( !vbuffer[i + nvpc] )
1670  {
1671  vbuffer[i + nvpc] = level_mesh[cur_level].start_vertex + count_nverts;
1672  count_nverts += 1;
1673  }
1674  }
1675 
1676  // Step 2: Coordinates of the new vertices
1677  std::vector< double > corner_coords( nvpc * 3 );
1678  error = get_coordinates( &cur_conn[0], nvpc, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
1679 
1680  error = compute_coordinates( cur_level, deg, type, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
1681  nverts_prev );MB_CHK_ERR( error );
1682 
1683  // Step 3: Choose the tet refine pattern to be used for this tet
1684  int diag = find_shortest_diagonal_octahedron( cur_level, deg, &vbuffer[0] );
1685  int pat_id = diag + 2;
1686  cell_patterns[cid] = pat_id;
1687 
1688  // Step 4: Use the template to obtain the subentities. The coordinates and local ahf maps
1689  // are also constructed. Connectivity of the children
1690  int etotal = refTemplates[pat_id][d].total_new_ents;
1691  std::vector< EntityHandle > ent_buffer( etotal );
1692 
1693  for( int i = 0; i < etotal; i++ )
1694  {
1695  for( int k = 0; k < nvpc; k++ )
1696  {
1697  int idx = refTemplates[pat_id][d].ents_conn[i][k];
1698  level_mesh[cur_level].cell_conn[nvpc * count_ents + k] = vbuffer[idx];
1699  }
1700  ent_buffer[i] = level_mesh[cur_level].start_cell + count_ents;
1701  count_ents += 1;
1702  }
1703 
1704  // Step 5: Update local ahf maps
1705  error = update_local_ahf( deg, MBTET, pat_id, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
1706 
1707  // Step 6: Update tracking information
1708  error = update_tracking_verts( cell, cur_level, deg, trackvertsC_edg, trackvertsC_face, &vbuffer[0] );MB_CHK_ERR( error );
1709  }
1710 
1711  // Step 7: Update the global maps
1712  // error = update_global_ahf(cur_level, deg, cell_patterns); MB_CHK_ERR(error);
1713  error = update_global_ahf( type, cur_level, deg, &cell_patterns );MB_CHK_ERR( error );
1714 
1715  // Step 8: If edges exists, refine them as well.
1716  if( level_mesh[cur_level].num_edges != 0 )
1717  {
1718  error = construct_hm_1D( cur_level, deg, type, trackvertsC_edg );MB_CHK_ERR( error );
1719  }
1720 
1721  // Step 9: If faces exists, refine them as well.
1722  if( !_infaces.empty() )
1723  {
1724  error = construct_hm_2D( cur_level, deg, type, trackvertsC_edg, trackvertsC_face );MB_CHK_ERR( error );
1725  }
1726 
1727  return MB_SUCCESS;
1728 }

References _incells, _infaces, _inverts, ahf, moab::Range::begin(), moab::NestedRefine::level_memory::cell_conn, compute_coordinates(), construct_hm_1D(), construct_hm_2D(), moab::Range::empty(), moab::NestedRefine::refPatterns::ents_conn, moab::error(), ErrorCode, find_shortest_diagonal_octahedron(), get_connectivity(), get_coordinates(), get_index_from_degree(), moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::lConnMap3D, level_mesh, MB_CHK_ERR, MB_SUCCESS, MBTET, moab::NestedRefine::level_memory::num_cells, moab::HalfFacetRep::LocalMaps3D::num_edges_in_cell, moab::HalfFacetRep::LocalMaps3D::num_faces_in_cell, moab::NestedRefine::level_memory::num_verts, moab::HalfFacetRep::LocalMaps3D::num_verts_in_cell, moab::NestedRefine::refPatterns::nv_edge, moab::NestedRefine::refPatterns::nv_face, refTemplates, moab::Range::size(), moab::NestedRefine::level_memory::start_cell, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, moab::NestedRefine::refPatterns::total_new_verts, update_global_ahf(), update_local_ahf(), update_tracking_verts(), moab::NestedRefine::refPatterns::vert_on_edges, and moab::NestedRefine::refPatterns::vert_on_faces.

Referenced by construct_hm_3D().

◆ update_ahf_1D()

ErrorCode moab::NestedRefine::update_ahf_1D ( int  cur_level)
protected

Definition at line 3398 of file NestedRefine.cpp.

3399 {
3400  ErrorCode error;
3401  error = ahf->determine_sibling_halfverts( level_mesh[cur_level].verts, level_mesh[cur_level].edges );MB_CHK_ERR( error );
3402 
3404 
3405  return MB_SUCCESS;
3406 }

References ahf, moab::HalfFacetRep::determine_incident_halfverts(), moab::HalfFacetRep::determine_sibling_halfverts(), moab::error(), ErrorCode, level_mesh, MB_CHK_ERR, and MB_SUCCESS.

Referenced by construct_hm_2D().

◆ update_global_ahf()

ErrorCode moab::NestedRefine::update_global_ahf ( EntityType  type,
int  cur_level,
int  deg,
std::vector< int > *  pattern_ids = NULL 
)
protected

Definition at line 3113 of file NestedRefine.cpp.

3114 {
3115 
3116  ErrorCode error;
3117 
3118  // Get the number of half-facets and number of children of each type
3119  if( type == MBEDGE )
3120  {
3121  assert( pattern_ids == NULL );
3122  error = update_global_ahf_1D( cur_level, deg );MB_CHK_ERR( error );
3123  }
3124  else if( type == MBTRI || type == MBQUAD )
3125  {
3126  assert( pattern_ids == NULL );
3127  error = update_global_ahf_2D( cur_level, deg );MB_CHK_ERR( error );
3128  }
3129  else if( type == MBHEX )
3130  {
3131  assert( pattern_ids == NULL );
3132  error = update_global_ahf_3D( cur_level, deg );MB_CHK_ERR( error );
3133  }
3134  else if( type == MBTET )
3135  {
3136  assert( pattern_ids != NULL );
3137  error = update_global_ahf_3D( cur_level, deg, pattern_ids );MB_CHK_ERR( error );
3138  }
3139  else
3140  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Requesting AHF update for an unsupported mesh entity type" );
3141 
3142  return MB_SUCCESS;
3143 }

References moab::error(), ErrorCode, MB_CHK_ERR, MB_NOT_IMPLEMENTED, MB_SET_ERR, MB_SUCCESS, MBEDGE, MBHEX, MBQUAD, MBTET, MBTRI, update_global_ahf_1D(), update_global_ahf_2D(), and update_global_ahf_3D().

Referenced by construct_hm_1D(), construct_hm_2D(), subdivide_cells(), and subdivide_tets().

◆ update_global_ahf_1D()

ErrorCode moab::NestedRefine::update_global_ahf_1D ( int  cur_level,
int  deg 
)
protected

Definition at line 3153 of file NestedRefine.cpp.

3154 {
3155  ErrorCode error;
3156  int d = get_index_from_degree( deg );
3157  int nhf, nchilds, nverts_prev, nents_prev;
3158  nhf = 2;
3159  nchilds = refTemplates[0][d].total_new_ents;
3160  if( cur_level )
3161  {
3162  nverts_prev = level_mesh[cur_level - 1].num_verts;
3163  nents_prev = level_mesh[cur_level - 1].num_edges;
3164  }
3165  else
3166  {
3167  nverts_prev = _inverts.size();
3168  nents_prev = _inedges.size();
3169  }
3170 
3171  std::vector< EntityHandle > inci_ent, child_ents;
3172  std::vector< int > inci_lid, child_lids;
3173 
3174  // Update the vertex to half-facet maps for duplicate vertices
3175  for( int i = 0; i < nverts_prev; i++ )
3176  {
3177  inci_ent.clear();
3178  inci_lid.clear();
3179  child_ents.clear();
3180  child_lids.clear();
3181 
3182  // Vertex id in the previous mesh and the current one
3183  EntityHandle vid;
3184  if( cur_level )
3185  vid = level_mesh[cur_level - 1].start_vertex + i;
3186  else
3187  vid = _inverts[i];
3188  EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
3189 
3190  // Get the incident half-vert in the previous mesh
3191  error = ahf->get_incident_map( MBEDGE, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
3192 
3193  // Obtain the corresponding incident child in the current mesh
3194  int lvid = get_local_vid( vid, inci_ent[0], cur_level - 1 );
3195  if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
3196  int chid = refTemplates[0][d].v2hf[lvid][0] - 1;
3197 
3198  int pid;
3199  if( cur_level )
3200  pid = inci_ent[0] - level_mesh[cur_level - 1].start_edge;
3201  else
3202  pid = inci_ent[0] - *_inedges.begin();
3203 
3204  int ind = nchilds * pid;
3205 
3206  child_ents.push_back( level_mesh[cur_level].start_edge + ind + chid );
3207  child_lids.push_back( refTemplates[0][d].v2hf[lvid][1] );
3208 
3209  error = ahf->set_incident_map( MBEDGE, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
3210  }
3211 
3212  // Update the sibling half-facet maps across entities
3213  for( int i = 0; i < nents_prev; i++ )
3214  {
3215  EntityHandle ent;
3216  if( cur_level )
3217  ent = level_mesh[cur_level - 1].start_edge + i;
3218  else
3219  ent = _inedges[i];
3220 
3221  std::vector< EntityHandle > sib_entids( nhf );
3222  std::vector< int > sib_lids( nhf );
3223 
3224  error = ahf->get_sibling_map( MBEDGE, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3225 
3226  int id, idx;
3227 
3228  for( int l = 0; l < nhf; l++ )
3229  {
3230  if( !sib_entids[l] ) continue;
3231 
3232  // Find the child incident on the half-facet
3233  id = refTemplates[0][d].ents_on_pent[l][1] - 1;
3234  idx = nchilds * i;
3235  EntityHandle child_ent = level_mesh[cur_level].start_edge + idx + id;
3236  int ch_lid = l;
3237 
3238  // Find the sibling of the child
3239  std::vector< EntityHandle > sib_childs( nhf );
3240  std::vector< int > sib_chlids( nhf );
3241 
3242  error = ahf->get_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
3243 
3244  // If the sibling already exists, dont do anything
3245  if( sib_childs[ch_lid] ) continue;
3246 
3247  // Get the correponding child of the sibling of the current parent
3248  int psib;
3249  if( cur_level )
3250  psib = sib_entids[l] - level_mesh[cur_level - 1].start_edge;
3251  else
3252  psib = sib_entids[l] - *_inedges.begin();
3253 
3254  int plid = sib_lids[l];
3255 
3256  id = refTemplates[0][d].ents_on_pent[plid][1] - 1;
3257  idx = nchilds * psib;
3258 
3259  EntityHandle psib_child = level_mesh[cur_level].start_edge + idx + id;
3260  int psib_chlid = plid;
3261 
3262  // Set the siblings
3263  sib_childs[ch_lid] = psib_child;
3264  sib_chlids[ch_lid] = psib_chlid;
3265 
3266  error = ahf->set_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
3267  }
3268  }
3269 
3270  return MB_SUCCESS;
3271 }

References _inedges, _inverts, ahf, moab::Range::begin(), moab::NestedRefine::refPatterns::ents_on_pent, moab::error(), ErrorCode, moab::HalfFacetRep::get_incident_map(), get_index_from_degree(), get_local_vid(), moab::HalfFacetRep::get_sibling_map(), level_mesh, MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, MBEDGE, moab::NestedRefine::level_memory::num_edges, moab::NestedRefine::level_memory::num_verts, refTemplates, moab::HalfFacetRep::set_incident_map(), moab::HalfFacetRep::set_sibling_map(), moab::Range::size(), moab::NestedRefine::level_memory::start_edge, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, and moab::NestedRefine::refPatterns::v2hf.

Referenced by update_global_ahf().

◆ update_global_ahf_1D_sub()

ErrorCode moab::NestedRefine::update_global_ahf_1D_sub ( int  cur_level,
int  deg 
)
protected

Definition at line 3273 of file NestedRefine.cpp.

3274 {
3275  ErrorCode error;
3276  int d = get_index_from_degree( deg );
3277  int nhf, nchilds, nents_prev;
3278  nhf = 2;
3279  nchilds = refTemplates[0][d].total_new_ents;
3280  if( cur_level )
3281  {
3282  nents_prev = level_mesh[cur_level - 1].num_edges;
3283  }
3284  else
3285  {
3286  nents_prev = _inedges.size();
3287  }
3288 
3289  // Update the sibling half-facet maps across entities
3290 
3291  std::vector< EntityHandle > conn;
3292  for( int i = 0; i < nents_prev; i++ )
3293  {
3294  EntityHandle ent;
3295  if( cur_level )
3296  ent = level_mesh[cur_level - 1].start_edge + i;
3297  else
3298  ent = _inedges[i];
3299 
3300  // Set incident hv maps
3301  conn.clear();
3302  error = get_connectivity( ent, cur_level, conn );MB_CHK_ERR( error );
3303 
3304  std::vector< EntityHandle > inci_ent, child_ents;
3305  std::vector< int > inci_lid, child_lids;
3306  for( int j = 0; j < 2; j++ )
3307  {
3308  inci_ent.clear();
3309  inci_lid.clear();
3310  child_ents.clear();
3311  child_lids.clear();
3312 
3313  // Get the entityhandle of the vertex from previous level in the current level
3314  EntityHandle cur_vid;
3315  if( cur_level )
3316  cur_vid = level_mesh[cur_level].start_vertex + ( conn[j] - level_mesh[cur_level - 1].start_vertex );
3317  else
3318  cur_vid = level_mesh[cur_level].start_vertex + ( conn[j] - *_inverts.begin() );
3319 
3320  // Obtain the incident half-facet. If exists, then no need to assign another
3321  error = ahf->get_incident_map( MBEDGE, cur_vid, inci_ent, inci_lid );MB_CHK_ERR( error );
3322  if( inci_ent[0] != 0 ) continue;
3323 
3324  // Get the incident half-facet on the old vertex
3325  error = ahf->get_incident_map( MBEDGE, conn[j], inci_ent, inci_lid );MB_CHK_ERR( error );
3326 
3327  // Obtain the corresponding incident child in the current mesh
3328  int lvid = get_local_vid( conn[j], inci_ent[0], cur_level - 1 );
3329  if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
3330  int chid = refTemplates[0][d].v2hf[lvid][0] - 1;
3331 
3332  int pid;
3333  if( cur_level )
3334  pid = inci_ent[0] - level_mesh[cur_level - 1].start_edge;
3335  else
3336  pid = inci_ent[0] - *_inedges.begin();
3337 
3338  int ind = nchilds * pid;
3339 
3340  child_ents.push_back( level_mesh[cur_level].start_edge + ind + chid );
3341  child_lids.push_back( refTemplates[0][d].v2hf[lvid][1] );
3342 
3343  error = ahf->set_incident_map( MBEDGE, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
3344  }
3345 
3346  std::vector< EntityHandle > sib_entids( nhf );
3347  std::vector< int > sib_lids( nhf );
3348 
3349  error = ahf->get_sibling_map( MBEDGE, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3350 
3351  int id, idx;
3352 
3353  for( int l = 0; l < nhf; l++ )
3354  {
3355  if( !sib_entids[l] ) continue;
3356 
3357  // Find the child incident on the half-facet
3358  id = refTemplates[0][d].ents_on_pent[l][1] - 1;
3359  idx = nchilds * i;
3360  EntityHandle child_ent = level_mesh[cur_level].start_edge + idx + id;
3361  int ch_lid = l;
3362 
3363  // Find the sibling of the child
3364  std::vector< EntityHandle > sib_childs( nhf );
3365  std::vector< int > sib_chlids( nhf );
3366 
3367  error = ahf->get_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
3368 
3369  // If the sibling already exists, dont do anything
3370  if( sib_childs[ch_lid] ) continue;
3371 
3372  // Get the correponding child of the sibling of the current parent
3373  int psib;
3374  if( cur_level )
3375  psib = sib_entids[l] - level_mesh[cur_level - 1].start_edge;
3376  else
3377  psib = sib_entids[l] - *_inedges.begin();
3378 
3379  int plid = sib_lids[l];
3380 
3381  id = refTemplates[0][d].ents_on_pent[plid][1] - 1;
3382  idx = nchilds * psib;
3383 
3384  EntityHandle psib_child = level_mesh[cur_level].start_edge + idx + id;
3385  int psib_chlid = plid;
3386 
3387  // Set the siblings
3388  sib_childs[ch_lid] = psib_child;
3389  sib_chlids[ch_lid] = psib_chlid;
3390 
3391  error = ahf->set_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
3392  }
3393  }
3394 
3395  return MB_SUCCESS;
3396 }

References _inedges, _inverts, ahf, moab::Range::begin(), moab::Range::clear(), moab::NestedRefine::refPatterns::ents_on_pent, moab::error(), ErrorCode, get_connectivity(), moab::HalfFacetRep::get_incident_map(), get_index_from_degree(), get_local_vid(), moab::HalfFacetRep::get_sibling_map(), level_mesh, MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, MBEDGE, moab::NestedRefine::level_memory::num_edges, refTemplates, moab::HalfFacetRep::set_incident_map(), moab::HalfFacetRep::set_sibling_map(), moab::Range::size(), moab::NestedRefine::level_memory::start_edge, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, and moab::NestedRefine::refPatterns::v2hf.

Referenced by construct_hm_1D().

◆ update_global_ahf_2D()

ErrorCode moab::NestedRefine::update_global_ahf_2D ( int  cur_level,
int  deg 
)
protected

Definition at line 3408 of file NestedRefine.cpp.

3409 {
3410  ErrorCode error;
3411 
3412  EntityType type = mbImpl->type_from_handle( *_infaces.begin() );
3413  int nhf, nchilds, nverts_prev, nents_prev;
3414 
3415  nhf = ahf->lConnMap2D[type - 2].num_verts_in_face;
3416  int d = get_index_from_degree( deg );
3417  nchilds = refTemplates[type - 1][d].total_new_ents;
3418 
3419  if( cur_level )
3420  {
3421  nverts_prev = level_mesh[cur_level - 1].num_verts;
3422  nents_prev = level_mesh[cur_level - 1].num_faces;
3423  }
3424  else
3425  {
3426  nverts_prev = _inverts.size();
3427  nents_prev = _infaces.size();
3428  }
3429 
3430  std::vector< EntityHandle > inci_ent, child_ents;
3431  std::vector< int > inci_lid, child_lids;
3432 
3433  // Update the vertex to half-edge maps for old/duplicate vertices
3434  for( int i = 0; i < nverts_prev; i++ )
3435  {
3436  inci_ent.clear();
3437  inci_lid.clear();
3438  child_ents.clear();
3439  child_lids.clear();
3440 
3441  // Vertex id in the previous mesh
3442  EntityHandle vid;
3443  if( cur_level )
3444  vid = level_mesh[cur_level - 1].start_vertex + i;
3445  else
3446  vid = _inverts[i];
3447  EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
3448 
3449  // Get the incident half-vert in the previous mesh
3450  error = ahf->get_incident_map( type, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
3451 
3452  // Obtain the corresponding incident child in the current mesh
3453  for( int j = 0; j < (int)inci_ent.size(); j++ )
3454  {
3455  int lvid = get_local_vid( vid, inci_ent[j], cur_level - 1 );
3456  if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
3457  int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
3458 
3459  int pid;
3460  if( cur_level )
3461  pid = inci_ent[j] - level_mesh[cur_level - 1].start_face;
3462  else
3463  pid = inci_ent[j] - *_infaces.begin();
3464 
3465  int ind = nchilds * pid;
3466 
3467  child_ents.push_back( level_mesh[cur_level].start_face + ind + chid );
3468  child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
3469  }
3470  error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
3471  }
3472 
3473  EntityHandle fedge[2];
3474 
3475  // Update the sibling half-facet maps across entities
3476  for( int i = 0; i < nents_prev; i++ )
3477  {
3478  EntityHandle ent;
3479  if( cur_level )
3480  ent = level_mesh[cur_level - 1].start_face + i;
3481  else
3482  ent = _infaces[i];
3483 
3484  std::vector< EntityHandle > fid_conn;
3485  error = get_connectivity( ent, cur_level, fid_conn );
3486  if( MB_SUCCESS != error ) return error;
3487 
3488  std::vector< EntityHandle > sib_entids( nhf );
3489  std::vector< int > sib_lids( nhf );
3490 
3491  error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3492 
3493  int id, idx;
3494 
3495  for( int l = 0; l < nhf; l++ )
3496  {
3497  if( !sib_entids[l] ) continue;
3498 
3499  int nidx = ahf->lConnMap2D[type - 2].next[l];
3500  fedge[0] = fid_conn[l];
3501  fedge[1] = fid_conn[nidx];
3502 
3503  EntityHandle sfid = sib_entids[l];
3504  int slid = sib_lids[l];
3505 
3506  std::vector< EntityHandle > conn;
3507  error = get_connectivity( sfid, cur_level, conn );
3508  if( MB_SUCCESS != error ) return error;
3509 
3510  bool orient = true;
3511  nidx = ahf->lConnMap2D[type - 2].next[slid];
3512  if( ( fedge[1] == conn[slid] ) && ( fedge[0] == conn[nidx] ) ) orient = false;
3513 
3514  if( orient ) assert( ( fedge[0] == conn[slid] ) && ( fedge[1] == conn[nidx] ) );
3515 
3516  // Find the childrens incident on the half-facet
3517  int nch = refTemplates[type - 1][d].ents_on_pent[l][0];
3518  idx = nchilds * i;
3519 
3520  // Loop over all the incident childrens
3521  for( int k = 0; k < nch; k++ )
3522  {
3523  id = refTemplates[type - 1][d].ents_on_pent[l][k + 1] - 1;
3524  EntityHandle child_ent = level_mesh[cur_level].start_face + idx + id;
3525  int child_lid = l;
3526 
3527  // Find the sibling of the child
3528  EntityHandle child_sibent;
3529  int child_siblid;
3530  error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
3531 
3532  if( child_sibent != 0 ) continue;
3533 
3534  // Get the correponding child of the sibling of the current parent
3535  int psib;
3536  if( cur_level )
3537  psib = sfid - level_mesh[cur_level - 1].start_face;
3538  else
3539  psib = sfid - *_infaces.begin();
3540 
3541  int plid = slid;
3542 
3543  if( orient )
3544  id = refTemplates[type - 1][d].ents_on_pent[plid][k + 1] - 1;
3545  else
3546  id = refTemplates[type - 1][d].ents_on_pent[plid][nch - k] - 1;
3547 
3548  int sidx = nchilds * psib;
3549 
3550  EntityHandle psib_child = level_mesh[cur_level].start_face + sidx + id;
3551  int psib_chlid = plid;
3552 
3553  // Set the siblings
3554  error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
3555  }
3556  }
3557  }
3558 
3559  return MB_SUCCESS;
3560 }

References _infaces, _inverts, ahf, moab::Range::begin(), moab::NestedRefine::refPatterns::ents_on_pent, moab::error(), ErrorCode, get_connectivity(), moab::HalfFacetRep::get_incident_map(), get_index_from_degree(), get_local_vid(), moab::HalfFacetRep::get_sibling_map(), moab::HalfFacetRep::lConnMap2D, level_mesh, MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, mbImpl, moab::HalfFacetRep::LocalMaps2D::next, moab::NestedRefine::level_memory::num_faces, moab::NestedRefine::level_memory::num_verts, moab::HalfFacetRep::LocalMaps2D::num_verts_in_face, refTemplates, moab::HalfFacetRep::set_incident_map(), moab::HalfFacetRep::set_sibling_map(), moab::Range::size(), moab::NestedRefine::level_memory::start_face, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, moab::Core::type_from_handle(), and moab::NestedRefine::refPatterns::v2hf.

Referenced by update_global_ahf().

◆ update_global_ahf_2D_sub()

ErrorCode moab::NestedRefine::update_global_ahf_2D_sub ( int  cur_level,
int  deg 
)
protected

Definition at line 3562 of file NestedRefine.cpp.

3563 {
3564  ErrorCode error;
3565  int d = get_index_from_degree( deg );
3566  EntityType type = mbImpl->type_from_handle( *_infaces.begin() );
3567  int nhf, nchilds, nents_prev;
3568  nhf = ahf->lConnMap2D[type - 2].num_verts_in_face;
3569  nchilds = refTemplates[type - 1][d].total_new_ents;
3570 
3571  if( cur_level )
3572  nents_prev = level_mesh[cur_level - 1].num_faces;
3573  else
3574  nents_prev = _infaces.size();
3575 
3576  EntityHandle fedge[2];
3577 
3578  // Update the sibling half-facet maps across entities
3579  for( int i = 0; i < nents_prev; i++ )
3580  {
3581  EntityHandle ent;
3582  if( cur_level )
3583  ent = level_mesh[cur_level - 1].start_face + i;
3584  else
3585  ent = _infaces[i];
3586 
3587  std::vector< EntityHandle > fid_conn;
3588  error = get_connectivity( ent, cur_level, fid_conn );
3589  if( MB_SUCCESS != error ) return error;
3590 
3591  std::vector< EntityHandle > inci_ent, child_ents;
3592  std::vector< int > inci_lid, child_lids;
3593 
3594  // Set incident half-edges
3595  for( int j = 0; j < nhf; j++ )
3596  {
3597  inci_ent.clear();
3598  inci_lid.clear();
3599  child_ents.clear();
3600  child_lids.clear();
3601  EntityHandle cur_vid;
3602  if( cur_level )
3603  cur_vid = level_mesh[cur_level].start_vertex + ( fid_conn[j] - level_mesh[cur_level - 1].start_vertex );
3604  else
3605  cur_vid = level_mesh[cur_level].start_vertex + ( fid_conn[j] - *_inverts.begin() );
3606 
3607  // Obtain the incident half-facet. If exists, then no need to assign another
3608  error = ahf->get_incident_map( type, cur_vid, inci_ent, inci_lid );MB_CHK_ERR( error );
3609  if( inci_ent[0] != 0 ) continue;
3610 
3611  // Get the incident half-facet on the old vertex
3612  error = ahf->get_incident_map( type, fid_conn[j], inci_ent, inci_lid );MB_CHK_ERR( error );
3613 
3614  // Obtain the corresponding incident child in the current mesh
3615  for( int k = 0; k < (int)inci_ent.size(); k++ )
3616  {
3617  int lvid = get_local_vid( fid_conn[j], inci_ent[k], cur_level - 1 );
3618  if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
3619  int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
3620 
3621  int pid;
3622  if( cur_level )
3623  pid = inci_ent[k] - level_mesh[cur_level - 1].start_face;
3624  else
3625  pid = inci_ent[k] - *_infaces.begin();
3626 
3627  int ind = nchilds * pid;
3628 
3629  child_ents.push_back( level_mesh[cur_level].start_face + ind + chid );
3630  child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
3631  }
3632 
3633  error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
3634  }
3635 
3636  // Set sibling half-edges
3637  std::vector< EntityHandle > sib_entids( nhf );
3638  std::vector< int > sib_lids( nhf );
3639 
3640  error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3641 
3642  int id, idx;
3643 
3644  for( int l = 0; l < nhf; l++ )
3645  {
3646  if( !sib_entids[l] ) continue;
3647 
3648  int nidx = ahf->lConnMap2D[type - 2].next[l];
3649  fedge[0] = fid_conn[l];
3650  fedge[1] = fid_conn[nidx];
3651 
3652  EntityHandle sfid = sib_entids[l];
3653  int slid = sib_lids[l];
3654 
3655  std::vector< EntityHandle > conn;
3656  error = get_connectivity( sfid, cur_level, conn );MB_CHK_ERR( error );
3657 
3658  assert( (int)conn.size() > nidx && (int)conn.size() > slid );
3659 
3660  bool orient = true;
3661  nidx = ahf->lConnMap2D[type - 2].next[slid];
3662  if( ( fedge[1] == conn[slid] ) && ( fedge[0] == conn[nidx] ) ) orient = false;
3663 
3664  if( orient ) assert( ( fedge[0] == conn[slid] ) && ( fedge[1] == conn[nidx] ) );
3665 
3666  // Find the childrens incident on the half-facet
3667  int nch = refTemplates[type - 1][d].ents_on_pent[l][0];
3668  idx = nchilds * i;
3669 
3670  // Loop over all the incident childrens
3671  for( int k = 0; k < nch; k++ )
3672  {
3673  id = refTemplates[type - 1][d].ents_on_pent[l][k + 1] - 1;
3674  EntityHandle child_ent = level_mesh[cur_level].start_face + idx + id;
3675  int child_lid = l;
3676 
3677  // Find the sibling of the child
3678  EntityHandle child_sibent;
3679  int child_siblid;
3680  error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
3681 
3682  if( child_sibent != 0 ) continue;
3683 
3684  // Get the correponding child of the sibling of the current parent
3685  int psib;
3686  if( cur_level )
3687  psib = sfid - level_mesh[cur_level - 1].start_face;
3688  else
3689  psib = sfid - *_infaces.begin();
3690 
3691  int plid = slid;
3692 
3693  if( orient )
3694  id = refTemplates[type - 1][d].ents_on_pent[plid][k + 1] - 1;
3695  else
3696  id = refTemplates[type - 1][d].ents_on_pent[plid][nch - k] - 1;
3697 
3698  int sidx = nchilds * psib;
3699 
3700  EntityHandle psib_child = level_mesh[cur_level].start_face + sidx + id;
3701  int psib_chlid = plid;
3702 
3703  // Set the siblings
3704  error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
3705  }
3706  }
3707  }
3708 
3709  return MB_SUCCESS;
3710 }

References _infaces, _inverts, ahf, moab::Range::begin(), moab::NestedRefine::refPatterns::ents_on_pent, moab::error(), ErrorCode, get_connectivity(), moab::HalfFacetRep::get_incident_map(), get_index_from_degree(), get_local_vid(), moab::HalfFacetRep::get_sibling_map(), moab::HalfFacetRep::lConnMap2D, level_mesh, MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, mbImpl, moab::HalfFacetRep::LocalMaps2D::next, moab::NestedRefine::level_memory::num_faces, moab::HalfFacetRep::LocalMaps2D::num_verts_in_face, refTemplates, moab::HalfFacetRep::set_incident_map(), moab::HalfFacetRep::set_sibling_map(), moab::Range::size(), moab::NestedRefine::level_memory::start_face, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, moab::Core::type_from_handle(), and moab::NestedRefine::refPatterns::v2hf.

Referenced by construct_hm_2D().

◆ update_global_ahf_3D()

ErrorCode moab::NestedRefine::update_global_ahf_3D ( int  cur_level,
int  deg,
std::vector< int > *  pattern_ids = NULL 
)
protected

Definition at line 3712 of file NestedRefine.cpp.

3713 {
3714  ErrorCode error;
3715  int nvpc, ne, nhf, nchilds, nverts_prev, nents_prev;
3716 
3717  EntityType type = mbImpl->type_from_handle( *_incells.begin() );
3718  int index = ahf->get_index_in_lmap( *_incells.begin() );
3719  int d = get_index_from_degree( deg );
3720 
3721  nhf = ahf->lConnMap3D[index].num_faces_in_cell;
3722  ne = ahf->lConnMap3D[index].num_edges_in_cell;
3723  nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
3724  nchilds = refTemplates[type - 1][d].total_new_ents;
3725 
3726  if( cur_level )
3727  {
3728  nverts_prev = level_mesh[cur_level - 1].num_verts;
3729  nents_prev = level_mesh[cur_level - 1].num_cells;
3730  }
3731  else
3732  {
3733  nverts_prev = _inverts.size();
3734  nents_prev = _incells.size();
3735  }
3736 
3737  std::vector< EntityHandle > inci_ent, child_ents;
3738  std::vector< int > inci_lid, child_lids;
3739 
3740  // Step 1: Update the V2HF maps for old/duplicate vertices
3741  for( int i = 0; i < nverts_prev; i++ )
3742  {
3743  inci_ent.clear();
3744  inci_lid.clear();
3745  child_ents.clear();
3746  child_lids.clear();
3747 
3748  // Vertex id in the previous mesh
3749  EntityHandle vid;
3750  if( cur_level )
3751  vid = level_mesh[cur_level - 1].start_vertex + i;
3752  else
3753  vid = _inverts[i];
3754  EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
3755 
3756  // Get the incident half-vert in the previous mesh
3757  error = ahf->get_incident_map( type, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
3758 
3759  // Obtain the corresponding incident child in the current mesh
3760  for( int j = 0; j < (int)inci_ent.size(); j++ )
3761  {
3762  int lvid = get_local_vid( vid, inci_ent[j], cur_level - 1 );
3763  if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
3764  int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
3765 
3766  int pid;
3767  if( cur_level )
3768  pid = inci_ent[j] - level_mesh[cur_level - 1].start_cell;
3769  else
3770  pid = inci_ent[j] - *_incells.begin();
3771 
3772  int ind = nchilds * pid;
3773 
3774  // EntityHandle child_ent = level_mesh[cur_level].start_cell + ind+chid ;
3775  // int child_lid = refTemplates[type-1][d].v2hf[lvid][1];
3776  child_ents.push_back( level_mesh[cur_level].start_cell + ind + chid );
3777  child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
3778  }
3779 
3780  error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
3781  }
3782 
3783  // error = ahf->determine_incident_halffaces( level_mesh[cur_level].cells);MB_CHK_ERR(error);
3784 
3785  // Step 2: Update SIBHFS maps
3786  for( int i = 0; i < nents_prev; i++ )
3787  {
3788  EntityHandle ent;
3789  if( cur_level )
3790  ent = level_mesh[cur_level - 1].start_cell + i;
3791  else
3792  ent = _incells[i];
3793 
3794  std::vector< EntityHandle > sib_entids( nhf );
3795  std::vector< int > sib_lids( nhf );
3796 
3797  error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3798 
3799  int id, idx;
3800 
3801  for( int l = 0; l < nhf; l++ )
3802  {
3803 
3804  if( !sib_entids[l] ) continue;
3805 
3806  // Get the number of children incident on this half-face
3807  int pat_id;
3808  if( type == MBTET )
3809  pat_id = ( *pattern_ids )[i];
3810  else
3811  pat_id = type - 1;
3812  int nch = refTemplates[pat_id][d].ents_on_pent[l][0];
3813 
3814  // Get the order of children indices incident on this half-face
3815  std::vector< int > id_sib( nch );
3816  for( int k = 0; k < nch; k++ )
3817  id_sib[k] = 0;
3818 
3819  error = reorder_indices( cur_level, deg, ent, l, sib_entids[l], sib_lids[l], 1, &id_sib[0] );MB_CHK_ERR( error );
3820 
3821  // Get the parent index of the sibling cell
3822  int psib;
3823  if( cur_level )
3824  psib = sib_entids[l] - level_mesh[cur_level - 1].start_cell;
3825  else
3826  psib = sib_entids[l] - *_incells.begin();
3827 
3828  int plid = sib_lids[l];
3829  int sidx = nchilds * psib;
3830  int sibpat_id;
3831  if( type == MBTET )
3832  sibpat_id = ( *pattern_ids )[psib];
3833  else
3834  sibpat_id = type - 1;
3835 
3836  // Loop over all the childs incident on the working half-face
3837  idx = nchilds * i;
3838 
3839  for( int k = 0; k < nch; k++ )
3840  {
3841  id = refTemplates[pat_id][d].ents_on_pent[l][k + 1] - 1;
3842  EntityHandle child_ent = level_mesh[cur_level].start_cell + idx + id;
3843  int child_lid = l;
3844 
3845  // Find the sibling of the working child
3846  EntityHandle child_sibent;
3847  int child_siblid;
3848  error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
3849 
3850  if( child_sibent != 0 ) continue;
3851 
3852  // Get the correponding child of the sibling of the current parent
3853  // We have already computed the order the children on incident corresponding to the
3854  // working half-face
3855  id = refTemplates[sibpat_id][d].ents_on_pent[plid][id_sib[k]] - 1;
3856 
3857  EntityHandle psib_child = level_mesh[cur_level].start_cell + sidx + id;
3858  int psib_chlid = plid;
3859 
3860  // Set the siblings of children incident on current half-face
3861  error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
3862 
3863  // Set the sibling of the sibling of the children to the children
3864  error = ahf->set_sibling_map( type, psib_child, psib_chlid, child_ent, child_lid );MB_CHK_ERR( error );
3865  }
3866  }
3867 
3868  // Loop over edges to check if there are any non-manifold edges. If there are then the v2hfs
3869  // map should be updated for the new vertices on it.
3870  const EntityHandle* conn;
3871  error = mbImpl->get_connectivity( ent, conn, nvpc );MB_CHK_ERR( error );
3872 
3873  int nv = refTemplates[type - 1][d].nv_edge; //#verts on each edge
3874  for( int l = 0; l < ne; l++ )
3875  {
3876  id = ahf->lConnMap3D[index].e2v[l][0];
3877  EntityHandle v_start = conn[id];
3878  id = ahf->lConnMap3D[index].e2v[l][1];
3879  EntityHandle v_end = conn[id];
3880 
3881  bool visited = false;
3882 
3883  std::vector< EntityHandle > inci_ent1, inci_ent2;
3884  std::vector< int > inci_lid1, inci_lid2;
3885  error = ahf->get_incident_map( type, v_start, inci_ent1, inci_lid1 );MB_CHK_ERR( error );
3886  error = ahf->get_incident_map( type, v_end, inci_ent2, inci_lid2 );MB_CHK_ERR( error );
3887 
3888  if( inci_ent1.size() > 1 && inci_ent2.size() > 1 )
3889  {
3890  std::vector< EntityHandle > cell_comps;
3891  std::vector< int > leid_comps;
3892 
3893  error = ahf->get_up_adjacencies_edg_3d_comp( ent, l, cell_comps, &leid_comps );MB_CHK_ERR( error );
3894 
3895  int ncomps = cell_comps.size();
3896  std::vector< EntityHandle > edgverts;
3897  std::vector< EntityHandle > compchildents( nv * ncomps );
3898  std::vector< int > compchildlfids( nv * ncomps );
3899 
3900  for( int s = 0; s < nv * ncomps; s++ )
3901  {
3902  compchildents[s] = 0;
3903  compchildlfids[s] = 0;
3904  }
3905 
3906  for( int j = 0; j < (int)cell_comps.size(); j++ )
3907  {
3908  int ind;
3909  if( cur_level )
3910  ind = level_mesh[cur_level - 1].cells.index( cell_comps[j] );
3911  else
3912  ind = _incells.index( cell_comps[j] );
3913 
3914  for( int k = 0; k < nv; k++ )
3915  {
3916  int chid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k] - 1;
3917  int lfid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k + 1];
3918  int lvid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k + 2];
3919 
3920  EntityHandle childcell = level_mesh[cur_level].start_cell + ind * nchilds + chid;
3921 
3922  const EntityHandle* econn;
3923  error = mbImpl->get_connectivity( childcell, econn, nvpc );MB_CHK_ERR( error );
3924 
3925  EntityHandle vert = econn[lvid];
3926 
3927  if( ahf->check_nonmanifold_vertices( type, vert ) )
3928  {
3929  visited = true;
3930  break;
3931  }
3932 
3933  if( edgverts.empty() )
3934  {
3935  edgverts.push_back( vert );
3936  compchildents[0] = childcell;
3937  compchildlfids[0] = lfid;
3938  }
3939  else
3940  {
3941  std::vector< EntityHandle >::iterator it;
3942  it = find( edgverts.begin(), edgverts.end(), vert );
3943  int indx = it - edgverts.begin();
3944 
3945  if( it == edgverts.end() )
3946  {
3947  edgverts.push_back( vert );
3948  compchildents[k * ncomps] = childcell;
3949  compchildlfids[k * ncomps] = lfid;
3950  }
3951  else
3952  {
3953  compchildents[indx * ncomps + j] = childcell;
3954  compchildlfids[indx * ncomps + j] = lfid;
3955  }
3956  }
3957  }
3958  }
3959 
3960  if( visited )
3961  {
3962  break;
3963  }
3964 
3965  // Set the incident half-facet map
3966  for( int k = 0; k < nv; k++ )
3967  {
3968  std::vector< EntityHandle > set_childents;
3969  std::vector< int > set_childlfids;
3970  for( int j = 0; j < ncomps; j++ )
3971  {
3972  set_childents.push_back( compchildents[k * ncomps + j] );
3973  set_childlfids.push_back( compchildlfids[k * ncomps + j] );
3974  }
3975 
3976  error = ahf->set_incident_map( type, edgverts[k], set_childents, set_childlfids );MB_CHK_ERR( error );
3977  }
3978  }
3979  }
3980  }
3981 
3982  return MB_SUCCESS;
3983 }

References _incells, _inverts, ahf, moab::Range::begin(), moab::NestedRefine::level_memory::cells, moab::HalfFacetRep::check_nonmanifold_vertices(), moab::HalfFacetRep::LocalMaps3D::e2v, moab::NestedRefine::refPatterns::ents_on_pent, moab::NestedRefine::refPatterns::ents_on_vedge, moab::error(), ErrorCode, moab::Core::get_connectivity(), moab::HalfFacetRep::get_incident_map(), get_index_from_degree(), moab::HalfFacetRep::get_index_in_lmap(), get_local_vid(), moab::HalfFacetRep::get_sibling_map(), moab::HalfFacetRep::get_up_adjacencies_edg_3d_comp(), moab::Range::index(), moab::HalfFacetRep::lConnMap3D, level_mesh, MB_CHK_ERR, MB_SET_ERR, MB_SUCCESS, mbImpl, MBTET, moab::NestedRefine::level_memory::num_cells, moab::HalfFacetRep::LocalMaps3D::num_edges_in_cell, moab::HalfFacetRep::LocalMaps3D::num_faces_in_cell, moab::NestedRefine::level_memory::num_verts, moab::HalfFacetRep::LocalMaps3D::num_verts_in_cell, moab::NestedRefine::refPatterns::nv_edge, refTemplates, reorder_indices(), moab::HalfFacetRep::set_incident_map(), moab::HalfFacetRep::set_sibling_map(), moab::Range::size(), moab::NestedRefine::level_memory::start_cell, moab::NestedRefine::level_memory::start_vertex, moab::NestedRefine::refPatterns::total_new_ents, moab::Core::type_from_handle(), and moab::NestedRefine::refPatterns::v2hf.

Referenced by update_global_ahf().

◆ update_local_ahf() [1/2]

ErrorCode moab::NestedRefine::update_local_ahf ( int  deg,
EntityType  type,
EntityHandle vbuffer,
EntityHandle ent_buffer,
int  etotal 
)
protected

Definition at line 3100 of file NestedRefine.cpp.

3105 {
3106  ErrorCode error;
3107  assert( type != MBTET );
3108  error = update_local_ahf( deg, type, type - 1, vbuffer, ent_buffer, etotal );MB_CHK_ERR( error );
3109 
3110  return MB_SUCCESS;
3111 }

References moab::error(), ErrorCode, MB_CHK_ERR, MB_SUCCESS, and MBTET.

Referenced by construct_hm_1D(), construct_hm_2D(), subdivide_cells(), and subdivide_tets().

◆ update_local_ahf() [2/2]

ErrorCode moab::NestedRefine::update_local_ahf ( int  deg,
EntityType  type,
int  pat_id,
EntityHandle vbuffer,
EntityHandle ent_buffer,
int  etotal 
)
protected

Definition at line 3005 of file NestedRefine.cpp.

3011 {
3012  ErrorCode error;
3013  int nhf = 0, nv = 0, total_new_verts = 0;
3014  int d = get_index_from_degree( deg );
3015 
3016  // Get the number of half-facets
3017  if( type == MBEDGE )
3018  {
3019  nhf = 2;
3020  nv = 2;
3021  total_new_verts = refTemplates[0][d].total_new_verts;
3022  }
3023  else if( type == MBTRI || type == MBQUAD )
3024  {
3025  nhf = ahf->lConnMap2D[type - 2].num_verts_in_face;
3026  nv = nhf;
3027  total_new_verts = refTemplates[pat_id][d].total_new_verts;
3028  }
3029  else if( type == MBTET || type == MBHEX )
3030  {
3031  int index = ahf->get_index_in_lmap( *_incells.begin() );
3032  nhf = ahf->lConnMap3D[index].num_faces_in_cell;
3033  nv = ahf->lConnMap3D[index].num_verts_in_cell;
3034  total_new_verts = refTemplates[pat_id][d].total_new_verts;
3035  }
3036 
3037  std::vector< EntityHandle > ent;
3038  std::vector< int > lid;
3039 
3040  // Update the vertex to half-facet map
3041  for( int i = 0; i < total_new_verts; i++ )
3042  {
3043  ent.clear();
3044  lid.clear();
3045  EntityHandle vid = vbuffer[i + nv];
3046  error = ahf->get_incident_map( type, vid, ent, lid );MB_CHK_ERR( error );
3047 
3048  if( ent[0] ) continue;
3049 
3050  int id = refTemplates[pat_id][d].v2hf[i + nv][0] - 1;
3051  ent[0] = ent_buffer[id];
3052  lid[0] = refTemplates[pat_id][d].v2hf[i + nv][1];
3053 
3054  error = ahf->set_incident_map( type, vid, ent, lid );MB_CHK_ERR( error );
3055  }
3056 
3057  // Update the sibling half-facet map
3058  for( int i = 0; i < etotal; i++ )
3059  {
3060  std::vector< EntityHandle > sib_entids( nhf );
3061  std::vector< int > sib_lids( nhf );
3062 
3063  error = ahf->get_sibling_map( type, ent_buffer[i], &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3064 
3065  for( int l = 0; l < nhf; l++ )
3066  {
3067  if( sib_entids[l] ) continue;
3068 
3069  // Fill out the sibling values
3070  int id = refTemplates[pat_id][d].ents_opphfs[i][2 * l];
3071  if( id )
3072  {
3073  sib_entids[l] = ent_buffer[id - 1];
3074  sib_lids[l] = refTemplates[pat_id][d].ents_opphfs[i][2 * l + 1];
3075  }
3076  else
3077  {
3078  sib_entids[l] = 0;
3079  sib_lids[l] = 0;
3080  }
3081  }
3082 
3083  error = ahf->set_sibling_map( type, ent_buffer[i], &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3084 
3085  for( int l = 0; l < nhf; l++ )
3086  {
3087  if( sib_entids[l] )
3088  {
3089 
3090  EntityHandle set_entid = ent_buffer[i];
3091  int set_lid = l;
3092 
3093  error = ahf->set_sibling_map( type, sib_entids[l], sib_lids[l], set_entid, set_lid );MB_CHK_ERR( error );
3094  }
3095  }
3096  }
3097  return MB_SUCCESS;
3098 }

References _incells, ahf, moab::Range::begin(), moab::NestedRefine::refPatterns::ents_opphfs, moab::error(), ErrorCode, moab::HalfFacetRep::get_incident_map(), get_index_from_degree(), moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::get_sibling_map(), moab::HalfFacetRep::lConnMap2D, moab::HalfFacetRep::lConnMap3D, MB_CHK_ERR, MB_SUCCESS, MBEDGE, MBHEX, MBQUAD, MBTET, MBTRI, moab::HalfFacetRep::LocalMaps3D::num_faces_in_cell, moab::HalfFacetRep::LocalMaps3D::num_verts_in_cell, moab::HalfFacetRep::LocalMaps2D::num_verts_in_face, refTemplates, moab::HalfFacetRep::set_incident_map(), moab::HalfFacetRep::set_sibling_map(), moab::NestedRefine::refPatterns::total_new_verts, and moab::NestedRefine::refPatterns::v2hf.

◆ update_special_tags()

ErrorCode moab::NestedRefine::update_special_tags ( int  level,
EntityHandle lset 
)

Definition at line 496 of file NestedRefine.cpp.

497 {
498  assert( level > 0 && level < nlevels + 1 );
499 
501  std::vector< Tag > mtags( 3 );
502 
506 
507  for( int i = 0; i < 3; i++ )
508  {
509  // Gather sets of a particular tag
510  Range sets;
511  error = mbImpl->get_entities_by_type_and_tag( _rset, MBENTITYSET, &mtags[i], NULL, 1, sets );MB_CHK_ERR( error );
512 
513  // Loop over all sets, gather entities in each set and add their children at all levels to
514  // the set
515  Range set_ents;
516  Range::iterator set_it;
517  std::vector< EntityHandle > childs;
518 
519  for( set_it = sets.begin(); set_it != sets.end(); ++set_it )
520  {
521  // Get the entities in the set, recursively
522  set_ents.clear();
523  childs.clear();
524  error = mbImpl->get_entities_by_handle( *set_it, set_ents, true );MB_CHK_ERR( error );
525 
526  // Gather child entities at the input level
527  for( Range::iterator sit = set_ents.begin(); sit != set_ents.end(); sit++ )
528  {
529  EntityType type = mbImpl->type_from_handle( *sit );
530  if( type == MBVERTEX )
531  {
532  Range conn;
533  std::vector< EntityHandle > cents;
534  error = vertex_to_entities_down( *sit, 0, level, cents );MB_CHK_ERR( error );
535  error = mbImpl->get_connectivity( &cents[0], (int)cents.size(), conn, true );MB_CHK_ERR( error );
536  childs.insert( childs.end(), cents.begin(), cents.end() );
537  }
538  else
539  {
540  error = parent_to_child( *sit, 0, level, childs );MB_CHK_ERR( error );
541  }
542 
543  std::sort( childs.begin(), childs.end() );
544  childs.erase( std::unique( childs.begin(), childs.end() ), childs.end() );
545 
546  // Add child entities to tagged sets
547  error = mbImpl->add_entities( *set_it, &childs[0], childs.size() );MB_CHK_ERR( error );
548  }
549 
550  // Remove the coarse entities
551  error = mbImpl->remove_entities( *set_it, set_ents );MB_CHK_ERR( error );
552 
553  // Add
554  error = mbImpl->add_entities( lset, &( *set_it ), 1 );MB_CHK_ERR( error );
555  }
556  }
557  return MB_SUCCESS;
558 }

References _rset, moab::Core::add_entities(), moab::Range::begin(), moab::Range::clear(), DIRICHLET_SET_TAG_NAME, moab::Range::end(), moab::error(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_entities_by_handle(), moab::Core::get_entities_by_type_and_tag(), MATERIAL_SET_TAG_NAME, MB_CHK_ERR, MB_SUCCESS, MB_TYPE_INTEGER, MBENTITYSET, mbImpl, MBVERTEX, NEUMANN_SET_TAG_NAME, nlevels, parent_to_child(), moab::Core::remove_entities(), moab::Core::tag_get_handle(), moab::Core::type_from_handle(), and vertex_to_entities_down().

◆ update_tracking_verts()

ErrorCode moab::NestedRefine::update_tracking_verts ( EntityHandle  cid,
int  cur_level,
int  deg,
std::vector< EntityHandle > &  trackvertsC_edg,
std::vector< EntityHandle > &  trackvertsC_face,
EntityHandle vbuffer 
)
protected

Definition at line 4196 of file NestedRefine.cpp.

4202 {
4203  // The vertices in the vbuffer are added to appropriate edges and faces of cells that are
4204  // incident on the working cell.
4205  ErrorCode error;
4206 
4207  EntityHandle cstart_prev;
4208  if( cur_level )
4209  cstart_prev = level_mesh[cur_level - 1].start_cell;
4210  else
4211  cstart_prev = *_incells.begin();
4212 
4213  EntityType cell_type = mbImpl->type_from_handle( cstart_prev );
4214  int cindex = cell_type - 1;
4215  int d = get_index_from_degree( deg );
4216 
4217  int nve = refTemplates[cindex][d].nv_edge;
4218  int nvf = refTemplates[cindex][d].nv_face;
4219 
4220  int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
4221  int nepc = ahf->lConnMap3D[index].num_edges_in_cell;
4222  int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
4223 
4224  // Step 1: Add the vertices on an edge of the working cell to tracking array of incident cells.
4225  for( int i = 0; i < nepc; i++ )
4226  {
4227  // Add the vertices to edges of the current cell
4228  for( int j = 0; j < nve; j++ )
4229  {
4230  int id = refTemplates[cindex][d].vert_on_edges[i][j];
4231  int idx = cid - cstart_prev;
4232  int aid = idx * nve * nepc + nve * i + j;
4233 
4234  if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
4235  }
4236 
4237  // Obtain all the incident cells
4238  std::vector< EntityHandle > inc_cids;
4239  std::vector< int > inc_leids, inc_orient;
4240 
4241  error = ahf->get_up_adjacencies_edg_3d( cid, i, inc_cids, &inc_leids, &inc_orient );MB_CHK_ERR( error );
4242 
4243  if( inc_cids.size() == 1 ) continue;
4244 
4245  // Add the vertices to the edges of the incident cells
4246  for( int k = 0; k < (int)inc_cids.size(); k++ )
4247  {
4248  if( inc_cids[k] == cid ) continue;
4249 
4250  int idx = inc_cids[k] - cstart_prev;
4251 
4252  if( inc_orient[k] ) // Same edge direction as the current edge
4253  {
4254  for( int j = 0; j < nve; j++ )
4255  {
4256  int id = refTemplates[cindex][d].vert_on_edges[i][j];
4257  int aid = idx * nve * nepc + nve * inc_leids[k] + j;
4258 
4259  if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
4260  }
4261  }
4262  else
4263  {
4264  for( int j = 0; j < nve; j++ )
4265  {
4266  int id = refTemplates[cindex][d].vert_on_edges[i][nve - j - 1];
4267  int aid = idx * nve * nepc + nve * inc_leids[k] + j;
4268 
4269  if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
4270  }
4271  }
4272  }
4273  }
4274 
4275  // Step 2: Add the vertices on a face of the working cell to tracking array of incident cells.
4276  if( nvf )
4277  {
4278 
4279  for( int i = 0; i < nfpc; i++ )
4280  {
4281  // Add vertices to the tracking array of vertices on faces for the current cell
4282  std::vector< EntityHandle > face_vbuf( nvf, 0 );
4283  for( int j = 0; j < nvf; j++ )
4284  {
4285  int id = refTemplates[cindex][d].vert_on_faces[i][j];
4286  int idx = cid - cstart_prev;
4287  int aid = idx * nvf * nfpc + nvf * i + j;
4288 
4289  if( !trackvertsC_face[aid] ) trackvertsC_face[aid] = vbuffer[id];
4290 
4291  face_vbuf[j] = vbuffer[id];
4292  }
4293 
4294  // Obtain all the incident cells
4295  std::vector< EntityHandle > sib_cids;
4296  std::vector< int > sib_lfids;
4297  error = ahf->get_up_adjacencies_face_3d( cid, i, sib_cids, &sib_lfids );MB_CHK_ERR( error );
4298 
4299  if( sib_cids.size() == 1 ) continue;
4300 
4301  // Reorder the vertex local ids incident on the half-face
4302  std::vector< int > id_sib( nvf );
4303  for( int k = 0; k < nvf; k++ )
4304  id_sib[k] = 0;
4305 
4306  error = reorder_indices( cur_level, deg, sib_cids[1], sib_lfids[1], cid, i, 0, &id_sib[0] );MB_CHK_ERR( error );
4307 
4308  // Add vertices to the tracking array of vertices on faces for the sibling cell of the
4309  // current cell
4310  for( int j = 0; j < nvf; j++ )
4311  {
4312  int idx = sib_cids[1] - cstart_prev;
4313  int aid = idx * nvf * nfpc + nvf * sib_lfids[1] + j;
4314 
4315  if( !trackvertsC_face[aid] ) trackvertsC_face[aid] = face_vbuf[id_sib[j] - 1];
4316  }
4317  }
4318  }
4319  return MB_SUCCESS;
4320 }

References _incells, ahf, moab::Range::begin(), moab::error(), ErrorCode, get_index_from_degree(), moab::HalfFacetRep::get_index_in_lmap(), moab::HalfFacetRep::get_up_adjacencies_edg_3d(), moab::HalfFacetRep::get_up_adjacencies_face_3d(), moab::HalfFacetRep::lConnMap3D, level_mesh, MB_CHK_ERR, MB_SUCCESS, mbImpl, moab::HalfFacetRep::LocalMaps3D::num_edges_in_cell, moab::HalfFacetRep::LocalMaps3D::num_faces_in_cell, moab::NestedRefine::refPatterns::nv_edge, moab::NestedRefine::refPatterns::nv_face, refTemplates, reorder_indices(), moab::NestedRefine::level_memory::start_cell, moab::Core::type_from_handle(), moab::NestedRefine::refPatterns::vert_on_edges, and moab::NestedRefine::refPatterns::vert_on_faces.

Referenced by subdivide_cells(), and subdivide_tets().

◆ vertex_to_entities_down()

ErrorCode moab::NestedRefine::vertex_to_entities_down ( EntityHandle  vertex,
int  vert_level,
int  child_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

Parameters
vertexEntityHandle of the vertex
vert_levelMesh level of the vertex
child_levelMesh level from which child entities are requested
incident_entitiesVector containing entities from the child level

Definition at line 388 of file NestedRefine.cpp.

392 {
393  assert( vert_level < child_level );
395 
396  // Step 1: Get the incident entities at the current level
397  std::vector< EntityHandle > inents;
398  if( meshdim == 1 )
399  {
401  }
402  else if( meshdim == 2 )
403  {
405  }
406  else if( meshdim == 3 )
407  {
409  }
410 
411  // Step 2: Loop over all the incident entities at the current level and gather their parents
412  std::vector< EntityHandle > childs;
413  for( int i = 0; i < (int)inents.size(); i++ )
414  {
415  childs.clear();
416  EntityHandle ent = inents[i];
417  error = parent_to_child( ent, vert_level, child_level, childs );MB_CHK_ERR( error );
418  for( int j = 0; j < (int)childs.size(); j++ )
419  incident_entities.push_back( childs[j] );
420  }
421 
422  return MB_SUCCESS;
423 }

References ahf, moab::error(), ErrorCode, moab::HalfFacetRep::get_up_adjacencies_1d(), moab::HalfFacetRep::get_up_adjacencies_vert_2d(), moab::HalfFacetRep::get_up_adjacencies_vert_3d(), MB_CHK_ERR, MB_SUCCESS, meshdim, and parent_to_child().

Referenced by update_special_tags().

◆ vertex_to_entities_up()

ErrorCode moab::NestedRefine::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 entities from any previous levels that contains it.

Parameters
vertexEntityHandle of the vertex
vert_levelMesh level of the vertex
parent_levelMesh level from which entities containing vertex is requested
incident_entitiesVector containing entities from the parent level incident on the vertex

Definition at line 348 of file NestedRefine.cpp.

352 {
353  assert( vert_level > parent_level );
355 
356  // Step 1: Get the incident entities at the current level
357  std::vector< EntityHandle > inents;
358  if( meshdim == 1 )
359  {
361  }
362  else if( meshdim == 2 )
363  {
365  }
366  else if( meshdim == 3 )
367  {
369  }
370 
371  // Step 2: Loop over all the incident entities at the current level and gather their parents
372  for( int i = 0; i < (int)inents.size(); i++ )
373  {
374  EntityHandle ent = inents[i];
375  EntityHandle parent;
376  error = child_to_parent( ent, vert_level, parent_level, &parent );MB_CHK_ERR( error );
377  incident_entities.push_back( parent );
378  }
379 
380  // Step 3: Sort and remove duplicates
381  std::sort( incident_entities.begin(), incident_entities.end() );
382  incident_entities.erase( std::unique( incident_entities.begin(), incident_entities.end() ),
383  incident_entities.end() );
384 
385  return MB_SUCCESS;
386 }

References ahf, child_to_parent(), moab::error(), ErrorCode, moab::HalfFacetRep::get_up_adjacencies_1d(), moab::HalfFacetRep::get_up_adjacencies_vert_2d(), moab::HalfFacetRep::get_up_adjacencies_vert_3d(), MB_CHK_ERR, MB_SUCCESS, and meshdim.

Member Data Documentation

◆ _incells

◆ _inedges

◆ _infaces

◆ _inverts

◆ _rset

EntityHandle moab::NestedRefine::_rset
protected

◆ ahf

◆ deg_index

std::map< int, int > moab::NestedRefine::deg_index
protected

Definition at line 180 of file NestedRefine.hpp.

Referenced by get_index_from_degree(), and initialize().

◆ elementype

EntityType moab::NestedRefine::elementype
protected

◆ hasghost

bool moab::NestedRefine::hasghost
protected

Definition at line 181 of file NestedRefine.hpp.

Referenced by exchange_ghosts(), and initialize().

◆ intFacEdg

const NestedRefine::intFEdge moab::NestedRefine::intFacEdg
staticprotected
Initial value:
= {
{ { 3, { { 3, 4 }, { 4, 5 }, { 5, 3 } } },
{ 9, { { 8, 3 }, { 3, 9 }, { 9, 4 }, { 4, 5 }, { 5, 9 }, { 9, 8 }, { 7, 9 }, { 9, 6 }, { 6, 7 } } } },
{ { 4, { { 4, 8 }, { 7, 8 }, { 8, 6 }, { 8, 5 } } },
{ 12,
{ { 4, 12 },
{ 12, 15 },
{ 15, 9 },
{ 5, 13 },
{ 13, 14 },
{ 14, 8 },
{ 11, 12 },
{ 12, 13 },
{ 13, 6 },
{ 10, 15 },
{ 15, 14 },
{ 14, 7 } } } } }

Definition at line 241 of file NestedRefine.hpp.

Referenced by construct_hm_2D(), and estimate_hm_storage().

◆ level_dsequence

int moab::NestedRefine::level_dsequence[MAX_LEVELS]
protected

Definition at line 179 of file NestedRefine.hpp.

Referenced by child_to_parent(), generate_mesh_hierarchy(), and parent_to_child().

◆ level_mesh

◆ mbImpl

◆ meshdim

◆ nlevels

int moab::NestedRefine::nlevels
protected

Definition at line 178 of file NestedRefine.hpp.

Referenced by generate_hm(), generate_mesh_hierarchy(), and update_special_tags().

◆ pcomm

ParallelComm* moab::NestedRefine::pcomm
protected

Definition at line 170 of file NestedRefine.hpp.

Referenced by exchange_ghosts(), generate_hm(), initialize(), and NestedRefine().

◆ permutation

const NestedRefine::pmat moab::NestedRefine::permutation
staticprotected

Definition at line 333 of file NestedRefine.hpp.

Referenced by reorder_indices().

◆ refTemplates

◆ timeall

codeperf moab::NestedRefine::timeall

Definition at line 166 of file NestedRefine.hpp.

Referenced by generate_hm(), and main().

◆ tm

CpuTimer* moab::NestedRefine::tm
protected

Definition at line 172 of file NestedRefine.hpp.

Referenced by generate_hm(), initialize(), and ~NestedRefine().


The documentation for this class was generated from the following files: