54 int A = opts.A, B = opts.B, C = opts.C, M = opts.M, N = opts.N, K = opts.K;
55 int blockSize = opts.blockSize;
56 double xsize = opts.xsize, ysize = opts.ysize, zsize = opts.zsize;
57 bool newMergeMethod = opts.newMergeMethod;
58 bool quadratic = opts.quadratic;
59 bool tetra = opts.tetra;
60 bool adjEnts = opts.adjEnts;
61 bool parmerge = opts.parmerge;
65 bool keep_skins = opts.keep_skins;
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;
117 ReadUtilIface*
iface;
129 "Can't get handle id tag" );
135 "Can't get parallel partition tag" );
140 for(
int a = 0; a < A; a++ )
142 for(
int b = 0; b < B; b++ )
144 for(
int c = 0; c < C; c++ )
149 int num_nodes = blockSize1 * blockSize1 * blockSize1;
151 vector< double* > arrays;
153 MB_CHK_SET_ERR(
iface->get_node_coords( 3, num_nodes, 0, startv, arrays ),
"Can't get node coords" );
156 int x = m * A * q * blockSize + a * q * blockSize;
157 int y = n * B * q * blockSize + b * q * blockSize;
158 int z = k * C * q * blockSize + c * q * blockSize;
160 vector< int > gids( num_nodes );
161 vector< long > lgids( num_nodes );
162 Range verts( startv, startv + num_nodes - 1 );
163 for(
int kk = 0; kk < blockSize1; kk++ )
165 for(
int jj = 0; jj < blockSize1; jj++ )
167 for(
int ii = 0; ii < blockSize1; ii++ )
169 arrays[0][ix] = ( x + ii ) * dx;
170 arrays[1][ix] = ( y + jj ) * dy;
171 arrays[2][ix] = ( z + kk ) * dz;
172 gids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + ( z + kk ) * ( NX * NY );
174 lgids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + (long)( z + kk ) * ( NX * NY );
183 "Can't set global ids to vertices" );
187 "Can't set the new handle id tags" );
189 localVerts.merge( verts );
190 int num_hexas = blockSize * blockSize * blockSize;
191 int num_el = num_hexas * factor;
195 int num_v_per_elem = 8;
200 "Can't get element connectivity" );
206 "Can't get element connectivity" );
211 "Can't get element connectivity" );
214 Range cells( starte, starte + num_el - 1 );
218 int xe = m * A * blockSize + a * blockSize;
219 int ye = n * B * blockSize + b * blockSize;
220 int ze = k * C * blockSize + c * blockSize;
221 gids.resize( num_el );
222 lgids.resize( num_el );
224 for(
int kk = 0; kk < blockSize; kk++ )
226 for(
int jj = 0; jj < blockSize; jj++ )
228 for(
int ii = 0; ii < blockSize; ii++ )
230 EntityHandle corner = startv + q * ii + q * jj * ystride + q * kk * zstride;
232 gids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + ( ze + kk ) * ( nex * ney ) ) *
234 lgids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + (long)( ze + kk ) * ( nex * ney ) ) *
258 conn[ix + 1] = corner + 2;
259 conn[ix + 2] = corner + 2 + 2 * ystride;
260 conn[ix + 3] = corner + 2 * ystride;
261 conn[ix + 4] = corner + 2 * zstride;
262 conn[ix + 5] = corner + 2 + 2 * zstride;
263 conn[ix + 6] = corner + 2 + 2 * ystride + 2 * zstride;
264 conn[ix + 7] = corner + 2 * ystride + 2 * zstride;
265 conn[ix + 8] = corner + 1;
266 conn[ix + 9] = corner + 2 + ystride;
267 conn[ix + 10] = corner + 1 + 2 * ystride;
268 conn[ix + 11] = corner + ystride;
269 conn[ix + 12] = corner + zstride;
270 conn[ix + 13] = corner + 2 + zstride;
271 conn[ix + 14] = corner + 2 + 2 * ystride + zstride;
272 conn[ix + 15] = corner + 2 * ystride + zstride;
273 conn[ix + 16] = corner + 1 + 2 * zstride;
274 conn[ix + 17] = corner + 2 + ystride + 2 * zstride;
275 conn[ix + 18] = corner + 1 + 2 * ystride + 2 * zstride;
276 conn[ix + 19] = corner + ystride + 2 * zstride;
277 conn[ix + 20] = corner + 1 + zstride;
278 conn[ix + 21] = corner + 2 + ystride + zstride;
279 conn[ix + 22] = corner + 1 + 2 * ystride + zstride;
280 conn[ix + 23] = corner + ystride + zstride;
281 conn[ix + 24] = corner + 1 + ystride;
282 conn[ix + 25] = corner + 1 + ystride + 2 * zstride;
283 conn[ix + 26] = corner + 1 + ystride + zstride;
338 for(
int ff = 0; ff < factor - 1; ff++ )
340 gids[ie] = gids[ie - 1] + 1;
350 conn[ix + 1] = corner + 1;
351 conn[ix + 2] = corner + 1 + ystride;
352 conn[ix + 3] = corner + ystride;
353 conn[ix + 4] = corner + zstride;
354 conn[ix + 5] = corner + 1 + zstride;
355 conn[ix + 6] = corner + 1 + ystride + zstride;
356 conn[ix + 7] = corner + ystride + zstride;
366 all3dcells.merge( cells );
369 "Can't update adjacencies" );
383 "Can't set global ids to elements" );
388 int part_num = a + m * A + ( b + n * B ) * ( M * A ) + ( c + k * C ) * ( M * A * N * B );
390 wsets.insert( part_set );
398 pc->partition_sets() = wsets;
415 std::cout <<
"generate local mesh: " << ( clock() - tt ) / (
double)CLOCKS_PER_SEC <<
" seconds" << endl;
419 std::cout <<
"number of elements on rank 0: " << all3dcells.size() << endl;
420 std::cout <<
"Total number of elements " << all3dcells.size() * size << endl;
421 std::cout <<
"Element type: " << ( tetra ?
"MBTET" :
"MBHEX" )
422 <<
" order:" << ( quadratic ?
"quadratic" :
"linear" ) << endl;
429 MB_CHK_SET_ERR( mm.merge_using_integer_tag( localVerts, global_id_tag ),
"Can't merge" );
433 MB_CHK_SET_ERR( mm.merge_entities( all3dcells, 0.0001 ),
"Can't merge" );
438 std::cout <<
"merge locally: " << ( clock() - tt ) / (
double)CLOCKS_PER_SEC <<
" seconds" << endl;
446 for(
Range::iterator wsit = wsets.begin(); wsit != wsets.end(); ++wsit )
449 Range cells, edges, faces;
466 ParallelMergeMesh pm( pc, 0.00001 );
471 std::cout <<
"parallel mesh merge: " << ( clock() - tt ) / (
double)CLOCKS_PER_SEC <<
" seconds" << endl;
478 MB_CHK_SET_ERR( pc->resolve_shared_ents(
cset, -1, -1, &new_id_tag ),
"Can't resolve shared ents" );
482 std::cout <<
"resolve shared entities: " << ( clock() - tt ) / (
double)CLOCKS_PER_SEC <<
" seconds"
496 MB_CHK_SET_ERR( pc->delete_entities( toDelete ),
"Can't delete entities" );
501 std::cout <<
"delete edges and faces \n";
502 toDelete.print( std::cout );
504 std::cout << ( clock() - tt ) / (
double)CLOCKS_PER_SEC <<
" seconds" << endl;
520 std::cout <<
"exchange " << GL <<
" ghost layer(s) :" << ( clock() - tt ) / (
double)CLOCKS_PER_SEC
521 <<
" seconds" << endl;