54 int A = opts.
A, B = opts.
B, C = opts.
C, M = opts.
M, N = opts.
N, K = opts.
K;
59 bool tetra = opts.
tetra;
68 int rank = 0,
size = 1;
77 if( M * N * K !=
size )
79 if( 0 == rank ) std::cout <<
"M*N*K = " << M * N * K <<
" != size = " <<
size <<
"\n";
86 int leftover = rank % ( M * N );
91 int q = ( quadratic ) ? 2 : 1;
93 int factor = ( tetra ) ? 6 : 1;
95 double dx = xsize / ( A * M * blockSize * q );
96 double dy = ysize / ( B * N * blockSize * q );
97 double dz = zsize / ( C * K * blockSize * q );
99 int NX = ( q * M * A * blockSize + 1 );
100 int NY = ( q * N * B * blockSize + 1 );
101 int nex = M * A * blockSize;
102 int ney = N * B * blockSize;
104 int blockSize1 = q * blockSize + 1;
105 long num_total_verts = (long)NX * NY * ( K * C * blockSize + 1 );
108 std::cout <<
"Generate mesh on " <<
size <<
" processors \n";
109 std::cout <<
"Total number of vertices: " << num_total_verts <<
"\n";
112 int ystride = blockSize1;
114 int zstride = blockSize1 * blockSize1;
138 for(
int a = 0; a < A; a++ )
140 for(
int b = 0; b < B; b++ )
142 for(
int c = 0; c < C; c++ )
147 int num_nodes = blockSize1 * blockSize1 * blockSize1;
149 vector< double* > arrays;
151 rval =
iface->get_node_coords( 3, num_nodes, 0, startv, arrays );
MB_CHK_SET_ERR( rval,
"Can't get node coords" );
154 int x = m * A * q * blockSize + a * q * blockSize;
155 int y = n * B * q * blockSize + b * q * blockSize;
156 int z = k * C * q * blockSize + c * q * blockSize;
158 vector< int > gids( num_nodes );
159 vector< long > lgids( num_nodes );
160 Range verts( startv, startv + num_nodes - 1 );
161 for(
int kk = 0; kk < blockSize1; kk++ )
163 for(
int jj = 0; jj < blockSize1; jj++ )
165 for(
int ii = 0; ii < blockSize1; ii++ )
167 arrays[0][ix] = ( x + ii ) * dx;
168 arrays[1][ix] = ( y + jj ) * dy;
169 arrays[2][ix] = ( z + kk ) * dz;
170 gids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + ( z + kk ) * ( NX * NY );
172 lgids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + (long)( z + kk ) * ( NX * NY );
185 localVerts.
merge( verts );
186 int num_hexas = blockSize * blockSize * blockSize;
187 int num_el = num_hexas * factor;
191 int num_v_per_elem = 8;
195 rval =
iface->get_element_connect( num_el, 27,
MBHEX, 0, starte, conn );
MB_CHK_SET_ERR( rval,
"Can't get element connectivity" );
200 rval =
iface->get_element_connect( num_el, 4,
MBTET, 0, starte, conn );
MB_CHK_SET_ERR( rval,
"Can't get element connectivity" );
204 rval =
iface->get_element_connect( num_el, 8,
MBHEX, 0, starte, conn );
MB_CHK_SET_ERR( rval,
"Can't get element connectivity" );
207 Range cells( starte, starte + num_el - 1 );
211 int xe = m * A * blockSize + a * blockSize;
212 int ye = n * B * blockSize + b * blockSize;
213 int ze = k * C * blockSize + c * blockSize;
214 gids.resize( num_el );
215 lgids.resize( num_el );
217 for(
int kk = 0; kk < blockSize; kk++ )
219 for(
int jj = 0; jj < blockSize; jj++ )
221 for(
int ii = 0; ii < blockSize; ii++ )
223 EntityHandle corner = startv + q * ii + q * jj * ystride + q * kk * zstride;
225 gids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + ( ze + kk ) * ( nex * ney ) ) *
227 lgids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + (long)( ze + kk ) * ( nex * ney ) ) *
251 conn[ix + 1] = corner + 2;
252 conn[ix + 2] = corner + 2 + 2 * ystride;
253 conn[ix + 3] = corner + 2 * ystride;
254 conn[ix + 4] = corner + 2 * zstride;
255 conn[ix + 5] = corner + 2 + 2 * zstride;
256 conn[ix + 6] = corner + 2 + 2 * ystride + 2 * zstride;
257 conn[ix + 7] = corner + 2 * ystride + 2 * zstride;
258 conn[ix + 8] = corner + 1;
259 conn[ix + 9] = corner + 2 + ystride;
260 conn[ix + 10] = corner + 1 + 2 * ystride;
261 conn[ix + 11] = corner + ystride;
262 conn[ix + 12] = corner + zstride;
263 conn[ix + 13] = corner + 2 + zstride;
264 conn[ix + 14] = corner + 2 + 2 * ystride + zstride;
265 conn[ix + 15] = corner + 2 * ystride + zstride;
266 conn[ix + 16] = corner + 1 + 2 * zstride;
267 conn[ix + 17] = corner + 2 + ystride + 2 * zstride;
268 conn[ix + 18] = corner + 1 + 2 * ystride + 2 * zstride;
269 conn[ix + 19] = corner + ystride + 2 * zstride;
270 conn[ix + 20] = corner + 1 + zstride;
271 conn[ix + 21] = corner + 2 + ystride + zstride;
272 conn[ix + 22] = corner + 1 + 2 * ystride + zstride;
273 conn[ix + 23] = corner + ystride + zstride;
274 conn[ix + 24] = corner + 1 + ystride;
275 conn[ix + 25] = corner + 1 + ystride + 2 * zstride;
276 conn[ix + 26] = corner + 1 + ystride + zstride;
331 for(
int ff = 0; ff < factor - 1; ff++ )
333 gids[ie] = gids[ie - 1] + 1;
343 conn[ix + 1] = corner + 1;
344 conn[ix + 2] = corner + 1 + ystride;
345 conn[ix + 3] = corner + ystride;
346 conn[ix + 4] = corner + zstride;
347 conn[ix + 5] = corner + 1 + zstride;
348 conn[ix + 6] = corner + 1 + ystride + zstride;
349 conn[ix + 7] = corner + ystride + zstride;
359 all3dcells.
merge( cells );
361 rval =
iface->update_adjacencies( starte, num_el, num_v_per_elem, conn );
MB_CHK_SET_ERR( rval,
"Can't update adjacencies" );
379 int part_num = a + m * A + ( b + n * B ) * ( M * A ) + ( c + k * C ) * ( M * A * N * B );
389 pc->partition_sets() = wsets;
406 std::cout <<
"generate local mesh: " << ( clock() - tt ) / (
double)CLOCKS_PER_SEC <<
" seconds" << endl;
410 std::cout <<
"number of elements on rank 0: " << all3dcells.
size() << endl;
411 std::cout <<
"Total number of elements " << all3dcells.
size() *
size << endl;
412 std::cout <<
"Element type: " << ( tetra ?
"MBTET" :
"MBHEX" )
413 <<
" order:" << ( quadratic ?
"quadratic" :
"linear" ) << endl;
429 std::cout <<
"merge locally: " << ( clock() - tt ) / (
double)CLOCKS_PER_SEC <<
" seconds" << endl;
440 Range cells, edges, faces;
462 std::cout <<
"parallel mesh merge: " << ( clock() - tt ) / (
double)CLOCKS_PER_SEC <<
" seconds" << endl;
469 rval = pc->resolve_shared_ents(
cset, -1, -1, &new_id_tag );
MB_CHK_SET_ERR( rval,
"Can't resolve shared ents" );
473 std::cout <<
"resolve shared entities: " << ( clock() - tt ) / (
double)CLOCKS_PER_SEC <<
" seconds"
487 rval = pc->delete_entities( toDelete );
MB_CHK_SET_ERR( rval,
"Can't delete entities" );
492 std::cout <<
"delete edges and faces \n";
493 toDelete.
print( std::cout );
495 std::cout << ( clock() - tt ) / (
double)CLOCKS_PER_SEC <<
" seconds" << endl;
503 rval = pc->exchange_ghost_cells( 3,
511 std::cout <<
"exchange " << GL <<
" ghost layer(s) :" << ( clock() - tt ) / (
double)CLOCKS_PER_SEC
512 <<
" seconds" << endl;