 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
#include <MeshGeneration.hpp>
Collaboration diagram for moab::MeshGeneration:Classes | |
| struct | BrickOpts |
Public Member Functions | |
| MeshGeneration (Interface *mbi, EntityHandle rset=0) | |
| virtual | ~MeshGeneration () |
| ErrorCode | BrickInstance (BrickOpts &opts) |
Private Attributes | |
| Interface * | mb |
| EntityHandle | cset |
Definition at line 18 of file MeshGeneration.hpp.
| moab::MeshGeneration::MeshGeneration | ( | Interface * | mbi, |
| EntityHandle | rset = 0 |
||
| ) |
Definition at line 30 of file MeshGeneration.cpp.
References moab::ParallelComm::get_pcomm(), and mb.
: mb( impl ), #ifdef MOAB_HAVE_MPI pc( comm ), #endif cset( rset ) { // ErrorCode error; #ifdef MOAB_HAVE_MPI // Get the Parallel Comm instance to prepare all new sets to work in parallel // in case the user did not provide any arguments if( !comm ) pc = moab::ParallelComm::get_pcomm( mb, 0 ); #endif }
| moab::MeshGeneration::~MeshGeneration | ( | ) | [virtual] |
Definition at line 50 of file MeshGeneration.cpp.
{}
Definition at line 52 of file MeshGeneration.cpp.
References moab::MeshGeneration::BrickOpts::A, moab::Interface::add_entities(), moab::MeshGeneration::BrickOpts::adjEnts, moab::MeshGeneration::BrickOpts::B, moab::Range::begin(), moab::MeshGeneration::BrickOpts::blockSize, moab::MeshGeneration::BrickOpts::C, CC, moab::Interface::create_meshset(), cset, moab::E, moab::Range::end(), ErrorCode, moab::F, moab::Interface::get_adjacencies(), moab::ReadUtilIface::get_element_connect(), moab::Interface::get_entities_by_dimension(), moab::ReadUtilIface::get_node_coords(), moab::MeshGeneration::BrickOpts::GL, iface, moab::Range::insert(), moab::MeshGeneration::BrickOpts::K, moab::MeshGeneration::BrickOpts::keep_skins, moab::MeshGeneration::BrickOpts::M, mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TAG_SPARSE, MB_TYPE_INTEGER, MB_TYPE_OPAQUE, MBHEX, MBTET, moab::ParallelMergeMesh::merge(), moab::Range::merge(), moab::MergeMesh::merge_entities(), moab::MergeMesh::merge_using_integer_tag(), MESHSET_SET, moab::MeshGeneration::BrickOpts::N, moab::MeshGeneration::BrickOpts::newMergeMethod, moab::MeshGeneration::BrickOpts::parmerge, moab::Range::print(), moab::MeshGeneration::BrickOpts::quadratic, moab::Interface::query_interface(), moab::Interface::remove_entities(), size, moab::Range::size(), moab::Interface::tag_delete(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), moab::MeshGeneration::BrickOpts::tetra, moab::Interface::UNION, moab::ReadUtilIface::update_adjacencies(), moab::MeshGeneration::BrickOpts::xsize, moab::MeshGeneration::BrickOpts::ysize, and moab::MeshGeneration::BrickOpts::zsize.
Referenced by main().
{
int A = opts.A, B = opts.B, C = opts.C, M = opts.M, N = opts.N, K = opts.K;
int blockSize = opts.blockSize;
double xsize = opts.xsize, ysize = opts.ysize, zsize = opts.zsize; // The size of the region
bool newMergeMethod = opts.newMergeMethod;
bool quadratic = opts.quadratic;
bool tetra = opts.tetra;
bool adjEnts = opts.adjEnts;
bool parmerge = opts.parmerge;
#ifdef MOAB_HAVE_MPI
int GL = opts.GL; // number of ghost layers
bool keep_skins = opts.keep_skins;
#endif
int rank = 0, size = 1;
#ifndef _MSC_VER /* windows */
clock_t tt = clock();
#endif
#ifdef MOAB_HAVE_MPI
rank = pc->rank();
size = pc->size();
#endif
if( M * N * K != size )
{
if( 0 == rank ) std::cout << "M*N*K = " << M * N * K << " != size = " << size << "\n";
return MB_FAILURE;
}
// Determine m, n, k for processor rank
int m, n, k;
k = rank / ( M * N );
int leftover = rank % ( M * N );
n = leftover / M;
m = leftover % M;
// Used for nodes increments
int q = ( quadratic ) ? 2 : 1;
// Used for element increments
int factor = ( tetra ) ? 6 : 1;
double dx = xsize / ( A * M * blockSize * q ); // Distance between 2 nodes in x direction
double dy = ysize / ( B * N * blockSize * q ); // Distance between 2 nodes in y direction
double dz = zsize / ( C * K * blockSize * q ); // Distance between 2 nodes in z direction
int NX = ( q * M * A * blockSize + 1 );
int NY = ( q * N * B * blockSize + 1 );
int nex = M * A * blockSize; // Number of elements in x direction, used for global id on element
int ney = N * B * blockSize; // Number of elements in y direction ...
// int NZ = (K * C * blockSize + 1); // Not used
int blockSize1 = q * blockSize + 1; // Used for vertices
long num_total_verts = (long)NX * NY * ( K * C * blockSize + 1 );
if( 0 == rank )
{
std::cout << "Generate mesh on " << size << " processors \n";
std::cout << "Total number of vertices: " << num_total_verts << "\n";
}
// int xstride = 1;
int ystride = blockSize1;
int zstride = blockSize1 * blockSize1;
// Generate the block at (a, b, c); it will represent a partition, it will get a partition tag
ReadUtilIface* iface;
ErrorCode rval = mb->query_interface( iface );MB_CHK_SET_ERR( rval, "Can't get reader interface" );
Tag global_id_tag;
rval = mb->tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, global_id_tag );MB_CHK_SET_ERR( rval, "Can't get global id tag" );
// set global ids
Tag new_id_tag;
if( !parmerge )
{
rval =
mb->tag_get_handle( "HANDLEID", sizeof( long ), MB_TYPE_OPAQUE, new_id_tag, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_SET_ERR( rval, "Can't get handle id tag" );
}
Tag part_tag;
int dum_id = -1;
rval =
mb->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, part_tag, MB_TAG_CREAT | MB_TAG_SPARSE, &dum_id );MB_CHK_SET_ERR( rval, "Can't get parallel partition tag" );
Range wsets; // write only part sets
Range localVerts;
Range all3dcells;
for( int a = 0; a < A; a++ )
{
for( int b = 0; b < B; b++ )
{
for( int c = 0; c < C; c++ )
{
// We will generate (q*block + 1)^3 vertices, and block^3 hexas; q is 1 for linear,
// 2 for quadratic the global id of the vertices will come from m, n, k, a, b, c x
// will vary from m*A*q*block + a*q*block to m*A*q*block + (a+1)*q*block etc;
int num_nodes = blockSize1 * blockSize1 * blockSize1;
vector< double* > arrays;
EntityHandle startv;
rval = iface->get_node_coords( 3, num_nodes, 0, startv, arrays );MB_CHK_SET_ERR( rval, "Can't get node coords" );
// Will start with the lower corner:
int x = m * A * q * blockSize + a * q * blockSize;
int y = n * B * q * blockSize + b * q * blockSize;
int z = k * C * q * blockSize + c * q * blockSize;
int ix = 0;
vector< int > gids( num_nodes );
vector< long > lgids( num_nodes );
Range verts( startv, startv + num_nodes - 1 );
for( int kk = 0; kk < blockSize1; kk++ )
{
for( int jj = 0; jj < blockSize1; jj++ )
{
for( int ii = 0; ii < blockSize1; ii++ )
{
arrays[0][ix] = ( x + ii ) * dx;
arrays[1][ix] = ( y + jj ) * dy;
arrays[2][ix] = ( z + kk ) * dz;
gids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + ( z + kk ) * ( NX * NY );
if( !parmerge )
lgids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + (long)( z + kk ) * ( NX * NY );
// Set int tags, some nice values?
ix++;
}
}
}
rval = mb->tag_set_data( global_id_tag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids to vertices" );
if( !parmerge )
{
rval = mb->tag_set_data( new_id_tag, verts, &lgids[0] );MB_CHK_SET_ERR( rval, "Can't set the new handle id tags" );
}
localVerts.merge( verts );
int num_hexas = blockSize * blockSize * blockSize;
int num_el = num_hexas * factor;
EntityHandle starte; // Connectivity
EntityHandle* conn;
int num_v_per_elem = 8;
if( quadratic )
{
num_v_per_elem = 27;
rval = iface->get_element_connect( num_el, 27, MBHEX, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
}
else if( tetra )
{
num_v_per_elem = 4;
rval = iface->get_element_connect( num_el, 4, MBTET, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
}
else
{
rval = iface->get_element_connect( num_el, 8, MBHEX, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
}
Range cells( starte, starte + num_el - 1 ); // Should be elements
// Fill cells
ix = 0;
// Identify the elements at the lower corner, for their global ids
int xe = m * A * blockSize + a * blockSize;
int ye = n * B * blockSize + b * blockSize;
int ze = k * C * blockSize + c * blockSize;
gids.resize( num_el );
lgids.resize( num_el );
int ie = 0; // Index now in the elements, for global ids
for( int kk = 0; kk < blockSize; kk++ )
{
for( int jj = 0; jj < blockSize; jj++ )
{
for( int ii = 0; ii < blockSize; ii++ )
{
EntityHandle corner = startv + q * ii + q * jj * ystride + q * kk * zstride;
// These could overflow for large numbers
gids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + ( ze + kk ) * ( nex * ney ) ) *
factor; // 6 more for tetra
lgids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + (long)( ze + kk ) * ( nex * ney ) ) *
factor; // 6 more for tetra
// EntityHandle eh = starte + ie;
ie++;
if( quadratic )
{
// 4 ----- 19 ----- 7
// . | . |
// 16 25 18 |
// . | . |
// 5 ----- 17 ----- 6 |
// | 12 | 23 15
// | | |
// | 20 | 26 | 22 |
// | | |
// 13 21 | 14 |
// | 0 ----- 11 ----- 3
// | . | .
// | 8 24 | 10
// | . | .
// 1 ----- 9 ----- 2
//
conn[ix] = corner;
conn[ix + 1] = corner + 2;
conn[ix + 2] = corner + 2 + 2 * ystride;
conn[ix + 3] = corner + 2 * ystride;
conn[ix + 4] = corner + 2 * zstride;
conn[ix + 5] = corner + 2 + 2 * zstride;
conn[ix + 6] = corner + 2 + 2 * ystride + 2 * zstride;
conn[ix + 7] = corner + 2 * ystride + 2 * zstride;
conn[ix + 8] = corner + 1; // 0-1
conn[ix + 9] = corner + 2 + ystride; // 1-2
conn[ix + 10] = corner + 1 + 2 * ystride; // 2-3
conn[ix + 11] = corner + ystride; // 3-0
conn[ix + 12] = corner + zstride; // 0-4
conn[ix + 13] = corner + 2 + zstride; // 1-5
conn[ix + 14] = corner + 2 + 2 * ystride + zstride; // 2-6
conn[ix + 15] = corner + 2 * ystride + zstride; // 3-7
conn[ix + 16] = corner + 1 + 2 * zstride; // 4-5
conn[ix + 17] = corner + 2 + ystride + 2 * zstride; // 5-6
conn[ix + 18] = corner + 1 + 2 * ystride + 2 * zstride; // 6-7
conn[ix + 19] = corner + ystride + 2 * zstride; // 4-7
conn[ix + 20] = corner + 1 + zstride; // 0154
conn[ix + 21] = corner + 2 + ystride + zstride; // 1265
conn[ix + 22] = corner + 1 + 2 * ystride + zstride; // 2376
conn[ix + 23] = corner + ystride + zstride; // 0374
conn[ix + 24] = corner + 1 + ystride; // 0123
conn[ix + 25] = corner + 1 + ystride + 2 * zstride; // 4567
conn[ix + 26] = corner + 1 + ystride + zstride; // center
ix += 27;
}
else if( tetra )
{
// E H
// F G
//
// A D
// B C
EntityHandle AA = corner;
EntityHandle BB = corner + 1;
EntityHandle CC = corner + 1 + ystride;
EntityHandle D = corner + ystride;
EntityHandle E = corner + zstride;
EntityHandle F = corner + 1 + zstride;
EntityHandle G = corner + 1 + ystride + zstride;
EntityHandle H = corner + ystride + zstride;
// tet EDHG
conn[ix] = E;
conn[ix + 1] = D;
conn[ix + 2] = H;
conn[ix + 3] = G;
// tet ABCF
conn[ix + 4] = AA;
conn[ix + 5] = BB;
conn[ix + 6] = CC;
conn[ix + 7] = F;
// tet ADEF
conn[ix + 8] = AA;
conn[ix + 9] = D;
conn[ix + 10] = E;
conn[ix + 11] = F;
// tet CGDF
conn[ix + 12] = CC;
conn[ix + 13] = G;
conn[ix + 14] = D;
conn[ix + 15] = F;
// tet ACDF
conn[ix + 16] = AA;
conn[ix + 17] = CC;
conn[ix + 18] = D;
conn[ix + 19] = F;
// tet DGEF
conn[ix + 20] = D;
conn[ix + 21] = G;
conn[ix + 22] = E;
conn[ix + 23] = F;
ix += 24;
for( int ff = 0; ff < factor - 1; ff++ )
{
gids[ie] = gids[ie - 1] + 1; // 6 more for tetra
// eh = starte + ie;
ie++;
}
}
else
{ // Linear hex
conn[ix] = corner;
conn[ix + 1] = corner + 1;
conn[ix + 2] = corner + 1 + ystride;
conn[ix + 3] = corner + ystride;
conn[ix + 4] = corner + zstride;
conn[ix + 5] = corner + 1 + zstride;
conn[ix + 6] = corner + 1 + ystride + zstride;
conn[ix + 7] = corner + ystride + zstride;
ix += 8;
}
}
}
}
EntityHandle part_set;
rval = mb->create_meshset( MESHSET_SET, part_set );MB_CHK_SET_ERR( rval, "Can't create mesh set" );
rval = mb->add_entities( part_set, cells );MB_CHK_SET_ERR( rval, "Can't add entities to set" );
all3dcells.merge( cells );
// update adjacencies now, because some elements are new;
rval = iface->update_adjacencies( starte, num_el, num_v_per_elem, conn );MB_CHK_SET_ERR( rval, "Can't update adjacencies" );
// If needed, add all edges and faces
if( adjEnts )
{
// Generate all adj entities dimension 1 and 2 (edges and faces/ tri or qua)
Range edges, faces;
rval = mb->get_adjacencies( cells, 1, true, edges, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get edges" );
rval = mb->get_adjacencies( cells, 2, true, faces, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get faces" );
// rval = mb->add_entities(part_set, edges);MB_CHK_SET_ERR(rval, "Can't add
// edges to partition set"); rval = mb->add_entities(part_set,
// faces);MB_CHK_SET_ERR(rval, "Can't add faces to partition set");
}
rval = mb->tag_set_data( global_id_tag, cells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids to elements" );
if( !parmerge )
{
rval = mb->tag_set_data( new_id_tag, cells, &lgids[0] );MB_CHK_SET_ERR( rval, "Can't set new ids to elements" );
}
int part_num = a + m * A + ( b + n * B ) * ( M * A ) + ( c + k * C ) * ( M * A * N * B );
rval = mb->tag_set_data( part_tag, &part_set, 1, &part_num );MB_CHK_SET_ERR( rval, "Can't set part tag on set" );
wsets.insert( part_set );
}
}
}
mb->add_entities( cset, all3dcells );
rval = mb->add_entities( cset, wsets );MB_CHK_SET_ERR( rval, "Can't add entity sets" );
#ifdef MOAB_HAVE_MPI
pc->partition_sets() = wsets;
#endif
/*
// Before merge locally
rval = mb->write_file("test0.h5m", 0, ";;PARALLEL=WRITE_PART");MB_CHK_SET_ERR(rval, "Can't write
in parallel, before merging");
*/
// After the mesh is generated on each proc, merge the vertices
MergeMesh mm( mb );
// rval = mb->get_entities_by_dimension(0, 3, all3dcells);MB_CHK_SET_ERR(rval, "Can't get all 3d
// cells elements");
if( 0 == rank )
{
#ifndef _MSC_VER /* windows */
std::cout << "generate local mesh: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
tt = clock();
#endif
std::cout << "number of elements on rank 0: " << all3dcells.size() << endl;
std::cout << "Total number of elements " << all3dcells.size() * size << endl;
std::cout << "Element type: " << ( tetra ? "MBTET" : "MBHEX" )
<< " order:" << ( quadratic ? "quadratic" : "linear" ) << endl;
}
if( A * B * C != 1 )
{ // Merge needed
if( newMergeMethod )
{
rval = mm.merge_using_integer_tag( localVerts, global_id_tag );MB_CHK_SET_ERR( rval, "Can't merge" );
}
else
{
rval = mm.merge_entities( all3dcells, 0.0001 );MB_CHK_SET_ERR( rval, "Can't merge" );
}
#ifndef _MSC_VER /* windows */
if( 0 == rank )
{
std::cout << "merge locally: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
tt = clock();
}
#endif
}
// if adjEnts, add now to each set
if( adjEnts )
{
for( Range::iterator wsit = wsets.begin(); wsit != wsets.end(); ++wsit )
{
EntityHandle ws = *wsit; // write set
Range cells, edges, faces;
rval = mb->get_entities_by_dimension( ws, 3, cells );MB_CHK_SET_ERR( rval, "Can't get cells" );
rval = mb->get_adjacencies( cells, 1, false, edges, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get edges" );
rval = mb->get_adjacencies( cells, 2, false, faces, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get faces" );
rval = mb->add_entities( ws, edges );MB_CHK_SET_ERR( rval, "Can't add edges to partition set" );
rval = mb->add_entities( ws, faces );MB_CHK_SET_ERR( rval, "Can't add faces to partition set" );
}
}
#ifdef MOAB_HAVE_MPI
if( size > 1 )
{
// rval = mb->create_meshset(MESHSET_SET, mesh_set);MB_CHK_SET_ERR(rval, "Can't create new
// set");
if( parmerge )
{
ParallelMergeMesh pm( pc, 0.00001 );
rval = pm.merge();MB_CHK_SET_ERR( rval, "Can't resolve shared ents" );
#ifndef _MSC_VER /* windows */
if( 0 == rank )
{
std::cout << "parallel mesh merge: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
tt = clock();
}
#endif
}
else
{
rval = pc->resolve_shared_ents( cset, -1, -1, &new_id_tag );MB_CHK_SET_ERR( rval, "Can't resolve shared ents" );
#ifndef _MSC_VER /* windows */
if( 0 == rank )
{
std::cout << "resolve shared entities: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds"
<< endl;
tt = clock();
}
#endif
}
if( !keep_skins )
{ // Default is to delete the 1- and 2-dimensional entities
// Delete all quads and edges
Range toDelete;
rval = mb->get_entities_by_dimension( cset, 1, toDelete );MB_CHK_SET_ERR( rval, "Can't get edges" );
rval = mb->get_entities_by_dimension( cset, 2, toDelete );MB_CHK_SET_ERR( rval, "Can't get faces" );
rval = pc->delete_entities( toDelete );MB_CHK_SET_ERR( rval, "Can't delete entities" );
rval = mb->remove_entities( cset, toDelete );MB_CHK_SET_ERR( rval, "Can't remove entities from base set" );
if( 0 == rank )
{
std::cout << "delete edges and faces \n";
toDelete.print( std::cout );
#ifndef _MSC_VER /* windows */
std::cout << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
tt = clock();
#endif
}
}
// do some ghosting if required
if( GL > 0 )
{
rval = pc->exchange_ghost_cells( 3, // int ghost_dim
0, // int bridge_dim
GL, // int num_layers
0, // int addl_ents
true );MB_CHK_ERR( rval ); // bool store_remote_handles
#ifndef _MSC_VER /* windows */
if( 0 == rank )
{
std::cout << "exchange " << GL << " ghost layer(s) :" << ( clock() - tt ) / (double)CLOCKS_PER_SEC
<< " seconds" << endl;
tt = clock();
}
#endif
}
}
#endif
if( !parmerge )
{
rval = mb->tag_delete( new_id_tag );MB_CHK_SET_ERR( rval, "Can't delete new ID tag" );
}
return MB_SUCCESS;
}
EntityHandle moab::MeshGeneration::cset [private] |
Definition at line 53 of file MeshGeneration.hpp.
Referenced by BrickInstance().
Interface* moab::MeshGeneration::mb [private] |
Definition at line 49 of file MeshGeneration.hpp.
Referenced by BrickInstance(), and MeshGeneration().
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.