Branch data Line data Source code
1 : : /*
2 : : * MGen.cpp
3 : : *
4 : : */
5 : :
6 : : #include "moab/MeshGeneration.hpp"
7 : : #include "moab/MergeMesh.hpp"
8 : : #include <iostream>
9 : : #include <ctime>
10 : : #include <vector>
11 : : #ifdef WIN32 /* windows */
12 : : #include <time.h>
13 : : #endif
14 : :
15 : : #ifdef MOAB_HAVE_MPI
16 : : #include "moab_mpi.h"
17 : : #include "moab/ParallelComm.hpp"
18 : : #include "MBParallelConventions.h"
19 : : #include "moab/ParallelMergeMesh.hpp"
20 : : #endif
21 : : #include "moab/ReadUtilIface.hpp"
22 : :
23 : : using std::endl;
24 : : using std::string;
25 : : using std::vector;
26 : :
27 : : namespace moab
28 : : {
29 : :
30 : 1 : MeshGeneration::MeshGeneration( Interface* impl, ParallelComm* comm, EntityHandle rset )
31 : 1 : : mb( impl ), pc( comm ), cset( rset )
32 : : {
33 : : // ErrorCode error;
34 : :
35 : : #ifdef MOAB_HAVE_MPI
36 : : // Get the Parallel Comm instance to prepare all new sets to work in parallel
37 : : // in case the user did not provide any arguments
38 [ - + ]: 1 : if( !comm ) pc = moab::ParallelComm::get_pcomm( mb, 0 );
39 : : #endif
40 : 1 : }
41 : :
42 [ # # ]: 0 : MeshGeneration::~MeshGeneration() {}
43 : :
44 : 1 : ErrorCode MeshGeneration::BrickInstance( MeshGeneration::BrickOpts& opts )
45 : : {
46 : 1 : int A = opts.A, B = opts.B, C = opts.C, M = opts.M, N = opts.N, K = opts.K;
47 : 1 : int blockSize = opts.blockSize;
48 : 1 : double xsize = opts.xsize, ysize = opts.ysize, zsize = opts.zsize; // The size of the region
49 : 1 : int GL = opts.GL; // number of ghost layers
50 : :
51 : 1 : bool newMergeMethod = opts.newMergeMethod;
52 : 1 : bool quadratic = opts.quadratic;
53 : 1 : bool keep_skins = opts.keep_skins;
54 : 1 : bool tetra = opts.tetra;
55 : 1 : bool adjEnts = opts.adjEnts;
56 : 1 : bool parmerge = opts.parmerge;
57 : :
58 : 1 : int rank = 0, size = 1;
59 : 1 : clock_t tt = clock();
60 : :
61 : : #ifdef MOAB_HAVE_MPI
62 [ + - ]: 1 : rank = pc->rank();
63 [ + - ]: 1 : size = pc->size();
64 : : #endif
65 : :
66 [ - + ]: 1 : if( M * N * K != size )
67 : : {
68 [ # # ][ # # ]: 0 : if( 0 == rank ) std::cout << "M*N*K = " << M * N * K << " != size = " << size << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
69 : :
70 : 0 : return MB_FAILURE;
71 : : }
72 : : // Determine m, n, k for processor rank
73 : : int m, n, k;
74 : 1 : k = rank / ( M * N );
75 : 1 : int leftover = rank % ( M * N );
76 : 1 : n = leftover / M;
77 : 1 : m = leftover % M;
78 : :
79 : : // Used for nodes increments
80 [ - + ]: 1 : int q = ( quadratic ) ? 2 : 1;
81 : : // Used for element increments
82 [ - + ]: 1 : int factor = ( tetra ) ? 6 : 1;
83 : :
84 : 1 : double dx = xsize / ( A * M * blockSize * q ); // Distance between 2 nodes in x direction
85 : 1 : double dy = ysize / ( B * N * blockSize * q ); // Distance between 2 nodes in y direction
86 : 1 : double dz = zsize / ( C * K * blockSize * q ); // Distance between 2 nodes in z direction
87 : :
88 : 1 : int NX = ( q * M * A * blockSize + 1 );
89 : 1 : int NY = ( q * N * B * blockSize + 1 );
90 : 1 : int nex = M * A * blockSize; // Number of elements in x direction, used for global id on element
91 : 1 : int ney = N * B * blockSize; // Number of elements in y direction ...
92 : : // int NZ = (K * C * blockSize + 1); // Not used
93 : 1 : int blockSize1 = q * blockSize + 1; // Used for vertices
94 : 1 : long num_total_verts = (long)NX * NY * ( K * C * blockSize + 1 );
95 [ + - ]: 1 : if( 0 == rank )
96 : : {
97 [ + - ][ + - ]: 1 : std::cout << "Generate mesh on " << size << " processors \n";
[ + - ]
98 [ + - ][ + - ]: 1 : std::cout << "Total number of vertices: " << num_total_verts << "\n";
[ + - ]
99 : : }
100 : : // int xstride = 1;
101 : 1 : int ystride = blockSize1;
102 : :
103 : 1 : int zstride = blockSize1 * blockSize1;
104 : : // Generate the block at (a, b, c); it will represent a partition, it will get a partition tag
105 : :
106 : : ReadUtilIface* iface;
107 [ + - ][ - + ]: 1 : ErrorCode rval = mb->query_interface( iface );MB_CHK_SET_ERR( rval, "Can't get reader interface" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
108 : :
109 : : Tag global_id_tag;
110 [ + - ][ - + ]: 1 : 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" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
111 : :
112 : : // set global ids
113 : : Tag new_id_tag;
114 [ + - ]: 1 : if( !parmerge )
115 : : {
116 : : rval =
117 [ + - ][ - + ]: 1 : 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" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
118 : : }
119 : : Tag part_tag;
120 : 1 : int dum_id = -1;
121 : : rval =
122 [ + - ][ - + ]: 1 : 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" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
123 : :
124 [ + - ]: 1 : Range wsets; // write only part sets
125 [ + - ]: 2 : Range localVerts;
126 [ + - ]: 2 : Range all3dcells;
127 [ + + ]: 3 : for( int a = 0; a < A; a++ )
128 : : {
129 [ + + ]: 6 : for( int b = 0; b < B; b++ )
130 : : {
131 [ + + ]: 12 : for( int c = 0; c < C; c++ )
132 : : {
133 : : // We will generate (q*block + 1)^3 vertices, and block^3 hexas; q is 1 for linear,
134 : : // 2 for quadratic the global id of the vertices will come from m, n, k, a, b, c x
135 : : // will vary from m*A*q*block + a*q*block to m*A*q*block + (a+1)*q*block etc;
136 : 8 : int num_nodes = blockSize1 * blockSize1 * blockSize1;
137 : :
138 [ + - ]: 8 : vector< double* > arrays;
139 : : EntityHandle startv;
140 [ + - ][ - + ]: 8 : rval = iface->get_node_coords( 3, num_nodes, 0, startv, arrays );MB_CHK_SET_ERR( rval, "Can't get node coords" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
141 : :
142 : : // Will start with the lower corner:
143 : 8 : int x = m * A * q * blockSize + a * q * blockSize;
144 : 8 : int y = n * B * q * blockSize + b * q * blockSize;
145 : 8 : int z = k * C * q * blockSize + c * q * blockSize;
146 : 8 : int ix = 0;
147 [ + - ][ + - ]: 16 : vector< int > gids( num_nodes );
148 [ + - ][ + - ]: 16 : vector< long > lgids( num_nodes );
149 [ + - ][ + - ]: 16 : Range verts( startv, startv + num_nodes - 1 );
150 [ + + ]: 48 : for( int kk = 0; kk < blockSize1; kk++ )
151 : : {
152 [ + + ]: 240 : for( int jj = 0; jj < blockSize1; jj++ )
153 : : {
154 [ + + ]: 1200 : for( int ii = 0; ii < blockSize1; ii++ )
155 : : {
156 [ + - ]: 1000 : arrays[0][ix] = ( x + ii ) * dx;
157 [ + - ]: 1000 : arrays[1][ix] = ( y + jj ) * dy;
158 [ + - ]: 1000 : arrays[2][ix] = ( z + kk ) * dz;
159 [ + - ]: 1000 : gids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + ( z + kk ) * ( NX * NY );
160 [ + - ]: 1000 : if( !parmerge )
161 [ + - ]: 1000 : lgids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + (long)( z + kk ) * ( NX * NY );
162 : : // Set int tags, some nice values?
163 : :
164 : 1000 : ix++;
165 : : }
166 : : }
167 : : }
168 : :
169 [ + - ][ + - ]: 8 : rval = mb->tag_set_data( global_id_tag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids to vertices" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
170 [ + - ]: 8 : if( !parmerge )
171 : : {
172 [ + - ][ + - ]: 8 : rval = mb->tag_set_data( new_id_tag, verts, &lgids[0] );MB_CHK_SET_ERR( rval, "Can't set the new handle id tags" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
173 : : }
174 [ + - ]: 8 : localVerts.merge( verts );
175 : 8 : int num_hexas = blockSize * blockSize * blockSize;
176 : 8 : int num_el = num_hexas * factor;
177 : :
178 : : EntityHandle starte; // Connectivity
179 : : EntityHandle* conn;
180 : 8 : int num_v_per_elem = 8;
181 [ - + ]: 8 : if( quadratic )
182 : : {
183 : 0 : num_v_per_elem = 27;
184 [ # # ][ # # ]: 0 : rval = iface->get_element_connect( num_el, 27, MBHEX, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
185 : : }
186 [ - + ]: 8 : else if( tetra )
187 : : {
188 : 0 : num_v_per_elem = 4;
189 [ # # ][ # # ]: 0 : rval = iface->get_element_connect( num_el, 4, MBTET, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
190 : : }
191 : : else
192 : : {
193 [ + - ][ - + ]: 8 : rval = iface->get_element_connect( num_el, 8, MBHEX, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
194 : : }
195 : :
196 [ + - ][ + - ]: 16 : Range cells( starte, starte + num_el - 1 ); // Should be elements
197 : : // Fill cells
198 : 8 : ix = 0;
199 : : // Identify the elements at the lower corner, for their global ids
200 : 8 : int xe = m * A * blockSize + a * blockSize;
201 : 8 : int ye = n * B * blockSize + b * blockSize;
202 : 8 : int ze = k * C * blockSize + c * blockSize;
203 [ + - ]: 8 : gids.resize( num_el );
204 [ + - ]: 8 : lgids.resize( num_el );
205 : 8 : int ie = 0; // Index now in the elements, for global ids
206 [ + + ]: 40 : for( int kk = 0; kk < blockSize; kk++ )
207 : : {
208 [ + + ]: 160 : for( int jj = 0; jj < blockSize; jj++ )
209 : : {
210 [ + + ]: 640 : for( int ii = 0; ii < blockSize; ii++ )
211 : : {
212 : 512 : EntityHandle corner = startv + q * ii + q * jj * ystride + q * kk * zstride;
213 : : // These could overflow for large numbers
214 [ + - ]: 512 : gids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + ( ze + kk ) * ( nex * ney ) ) *
215 : 512 : factor; // 6 more for tetra
216 [ + - ]: 512 : lgids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + (long)( ze + kk ) * ( nex * ney ) ) *
217 : 512 : factor; // 6 more for tetra
218 : : // EntityHandle eh = starte + ie;
219 : :
220 : 512 : ie++;
221 [ - + ]: 512 : if( quadratic )
222 : : {
223 : : // 4 ----- 19 ----- 7
224 : : // . | . |
225 : : // 16 25 18 |
226 : : // . | . |
227 : : // 5 ----- 17 ----- 6 |
228 : : // | 12 | 23 15
229 : : // | | |
230 : : // | 20 | 26 | 22 |
231 : : // | | |
232 : : // 13 21 | 14 |
233 : : // | 0 ----- 11 ----- 3
234 : : // | . | .
235 : : // | 8 24 | 10
236 : : // | . | .
237 : : // 1 ----- 9 ----- 2
238 : : //
239 : 0 : conn[ix] = corner;
240 : 0 : conn[ix + 1] = corner + 2;
241 : 0 : conn[ix + 2] = corner + 2 + 2 * ystride;
242 : 0 : conn[ix + 3] = corner + 2 * ystride;
243 : 0 : conn[ix + 4] = corner + 2 * zstride;
244 : 0 : conn[ix + 5] = corner + 2 + 2 * zstride;
245 : 0 : conn[ix + 6] = corner + 2 + 2 * ystride + 2 * zstride;
246 : 0 : conn[ix + 7] = corner + 2 * ystride + 2 * zstride;
247 : 0 : conn[ix + 8] = corner + 1; // 0-1
248 : 0 : conn[ix + 9] = corner + 2 + ystride; // 1-2
249 : 0 : conn[ix + 10] = corner + 1 + 2 * ystride; // 2-3
250 : 0 : conn[ix + 11] = corner + ystride; // 3-0
251 : 0 : conn[ix + 12] = corner + zstride; // 0-4
252 : 0 : conn[ix + 13] = corner + 2 + zstride; // 1-5
253 : 0 : conn[ix + 14] = corner + 2 + 2 * ystride + zstride; // 2-6
254 : 0 : conn[ix + 15] = corner + 2 * ystride + zstride; // 3-7
255 : 0 : conn[ix + 16] = corner + 1 + 2 * zstride; // 4-5
256 : 0 : conn[ix + 17] = corner + 2 + ystride + 2 * zstride; // 5-6
257 : 0 : conn[ix + 18] = corner + 1 + 2 * ystride + 2 * zstride; // 6-7
258 : 0 : conn[ix + 19] = corner + ystride + 2 * zstride; // 4-7
259 : 0 : conn[ix + 20] = corner + 1 + zstride; // 0154
260 : 0 : conn[ix + 21] = corner + 2 + ystride + zstride; // 1265
261 : 0 : conn[ix + 22] = corner + 1 + 2 * ystride + zstride; // 2376
262 : 0 : conn[ix + 23] = corner + ystride + zstride; // 0374
263 : 0 : conn[ix + 24] = corner + 1 + ystride; // 0123
264 : 0 : conn[ix + 25] = corner + 1 + ystride + 2 * zstride; // 4567
265 : 0 : conn[ix + 26] = corner + 1 + ystride + zstride; // center
266 : 0 : ix += 27;
267 : : }
268 [ - + ]: 512 : else if( tetra )
269 : : {
270 : : // E H
271 : : // F G
272 : : //
273 : : // A D
274 : : // B C
275 : 0 : EntityHandle AA = corner;
276 : 0 : EntityHandle BB = corner + 1;
277 : 0 : EntityHandle CC = corner + 1 + ystride;
278 : 0 : EntityHandle D = corner + ystride;
279 : 0 : EntityHandle E = corner + zstride;
280 : 0 : EntityHandle F = corner + 1 + zstride;
281 : 0 : EntityHandle G = corner + 1 + ystride + zstride;
282 : 0 : EntityHandle H = corner + ystride + zstride;
283 : :
284 : : // tet EDHG
285 : 0 : conn[ix] = E;
286 : 0 : conn[ix + 1] = D;
287 : 0 : conn[ix + 2] = H;
288 : 0 : conn[ix + 3] = G;
289 : :
290 : : // tet ABCF
291 : 0 : conn[ix + 4] = AA;
292 : 0 : conn[ix + 5] = BB;
293 : 0 : conn[ix + 6] = CC;
294 : 0 : conn[ix + 7] = F;
295 : :
296 : : // tet ADEF
297 : 0 : conn[ix + 8] = AA;
298 : 0 : conn[ix + 9] = D;
299 : 0 : conn[ix + 10] = E;
300 : 0 : conn[ix + 11] = F;
301 : :
302 : : // tet CGDF
303 : 0 : conn[ix + 12] = CC;
304 : 0 : conn[ix + 13] = G;
305 : 0 : conn[ix + 14] = D;
306 : 0 : conn[ix + 15] = F;
307 : :
308 : : // tet ACDF
309 : 0 : conn[ix + 16] = AA;
310 : 0 : conn[ix + 17] = CC;
311 : 0 : conn[ix + 18] = D;
312 : 0 : conn[ix + 19] = F;
313 : :
314 : : // tet DGEF
315 : 0 : conn[ix + 20] = D;
316 : 0 : conn[ix + 21] = G;
317 : 0 : conn[ix + 22] = E;
318 : 0 : conn[ix + 23] = F;
319 : 0 : ix += 24;
320 [ # # ]: 0 : for( int ff = 0; ff < factor - 1; ff++ )
321 : : {
322 [ # # ][ # # ]: 0 : gids[ie] = gids[ie - 1] + 1; // 6 more for tetra
323 : :
324 : : // eh = starte + ie;
325 : :
326 : 0 : ie++;
327 : : }
328 : : }
329 : : else
330 : : { // Linear hex
331 : 512 : conn[ix] = corner;
332 : 512 : conn[ix + 1] = corner + 1;
333 : 512 : conn[ix + 2] = corner + 1 + ystride;
334 : 512 : conn[ix + 3] = corner + ystride;
335 : 512 : conn[ix + 4] = corner + zstride;
336 : 512 : conn[ix + 5] = corner + 1 + zstride;
337 : 512 : conn[ix + 6] = corner + 1 + ystride + zstride;
338 : 512 : conn[ix + 7] = corner + ystride + zstride;
339 : 512 : ix += 8;
340 : : }
341 : : }
342 : : }
343 : : }
344 : :
345 : : EntityHandle part_set;
346 [ + - ][ - + ]: 8 : rval = mb->create_meshset( MESHSET_SET, part_set );MB_CHK_SET_ERR( rval, "Can't create mesh set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
347 [ + - ][ - + ]: 8 : rval = mb->add_entities( part_set, cells );MB_CHK_SET_ERR( rval, "Can't add entities to set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
348 [ + - ]: 8 : all3dcells.merge( cells );
349 : : // update adjacencies now, because some elements are new;
350 [ + - ][ - + ]: 8 : rval = iface->update_adjacencies( starte, num_el, num_v_per_elem, conn );MB_CHK_SET_ERR( rval, "Can't update adjacencies" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
351 : : // If needed, add all edges and faces
352 [ - + ]: 8 : if( adjEnts )
353 : : {
354 : : // Generate all adj entities dimension 1 and 2 (edges and faces/ tri or qua)
355 [ # # ][ # # ]: 0 : Range edges, faces;
[ # # ]
356 [ # # ][ # # ]: 0 : rval = mb->get_adjacencies( cells, 1, true, edges, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get edges" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
357 [ # # ][ # # ]: 0 : rval = mb->get_adjacencies( cells, 2, true, faces, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get faces" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
358 : : // rval = mb->add_entities(part_set, edges);MB_CHK_SET_ERR(rval, "Can't add
359 : : // edges to partition set"); rval = mb->add_entities(part_set,
360 : : // faces);MB_CHK_SET_ERR(rval, "Can't add faces to partition set");
361 : : }
362 : :
363 [ + - ][ + - ]: 8 : rval = mb->tag_set_data( global_id_tag, cells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids to elements" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
364 [ + - ]: 8 : if( !parmerge )
365 : : {
366 [ + - ][ + - ]: 8 : rval = mb->tag_set_data( new_id_tag, cells, &lgids[0] );MB_CHK_SET_ERR( rval, "Can't set new ids to elements" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
367 : : }
368 : 8 : int part_num = a + m * A + ( b + n * B ) * ( M * A ) + ( c + k * C ) * ( M * A * N * B );
369 [ + - ][ - + ]: 8 : rval = mb->tag_set_data( part_tag, &part_set, 1, &part_num );MB_CHK_SET_ERR( rval, "Can't set part tag on set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
370 [ + - ][ + - ]: 8 : wsets.insert( part_set );
371 : 8 : }
372 : : }
373 : : }
374 : :
375 [ + - ]: 1 : mb->add_entities( cset, all3dcells );
376 [ + - ][ - + ]: 1 : rval = mb->add_entities( cset, wsets );MB_CHK_SET_ERR( rval, "Can't add entity sets" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
377 : : #ifdef MOAB_HAVE_MPI
378 [ + - ][ + - ]: 1 : pc->partition_sets() = wsets;
379 : : #endif
380 : :
381 : : /*
382 : : // Before merge locally
383 : : rval = mb->write_file("test0.h5m", 0, ";;PARALLEL=WRITE_PART");MB_CHK_SET_ERR(rval, "Can't write
384 : : in parallel, before merging");
385 : : */
386 : : // After the mesh is generated on each proc, merge the vertices
387 [ + - ]: 2 : MergeMesh mm( mb );
388 : :
389 : : // rval = mb->get_entities_by_dimension(0, 3, all3dcells);MB_CHK_SET_ERR(rval, "Can't get all 3d
390 : : // cells elements");
391 : :
392 [ + - ]: 1 : if( 0 == rank )
393 : : {
394 [ + - ][ + - ]: 1 : std::cout << "generate local mesh: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
[ + - ][ + - ]
395 : 1 : tt = clock();
396 [ + - ][ + - ]: 1 : std::cout << "number of elements on rank 0: " << all3dcells.size() << endl;
[ + - ][ + - ]
397 [ + - ][ + - ]: 1 : std::cout << "Total number of elements " << all3dcells.size() * size << endl;
[ + - ][ + - ]
398 [ - + ][ + - ]: 1 : std::cout << "Element type: " << ( tetra ? "MBTET" : "MBHEX" )
[ + - ]
399 [ - + ][ + - ]: 2 : << " order:" << ( quadratic ? "quadratic" : "linear" ) << endl;
[ + - ][ + - ]
400 : : }
401 : :
402 [ + - ]: 1 : if( A * B * C != 1 )
403 : : { // Merge needed
404 [ - + ]: 1 : if( newMergeMethod )
405 : : {
406 [ # # ][ # # ]: 0 : rval = mm.merge_using_integer_tag( localVerts, global_id_tag );MB_CHK_SET_ERR( rval, "Can't merge" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
407 : : }
408 : : else
409 : : {
410 [ + - ][ - + ]: 1 : rval = mm.merge_entities( all3dcells, 0.0001 );MB_CHK_SET_ERR( rval, "Can't merge" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
411 : : }
412 : :
413 [ + - ]: 1 : if( 0 == rank )
414 : : {
415 [ + - ][ + - ]: 1 : std::cout << "merge locally: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
[ + - ][ + - ]
416 : 1 : tt = clock();
417 : : }
418 : : }
419 : : // if adjEnts, add now to each set
420 [ - + ]: 1 : if( adjEnts )
421 : : {
422 [ # # ][ # # ]: 0 : for( Range::iterator wsit = wsets.begin(); wsit != wsets.end(); ++wsit )
[ # # ][ # # ]
[ # # ]
423 : : {
424 [ # # ]: 0 : EntityHandle ws = *wsit; // write set
425 [ # # ][ # # ]: 0 : Range cells, edges, faces;
[ # # ][ # # ]
[ # # ]
426 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( ws, 3, cells );MB_CHK_SET_ERR( rval, "Can't get cells" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
427 [ # # ][ # # ]: 0 : rval = mb->get_adjacencies( cells, 1, false, edges, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get edges" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
428 [ # # ][ # # ]: 0 : rval = mb->get_adjacencies( cells, 2, false, faces, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get faces" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
429 [ # # ][ # # ]: 0 : rval = mb->add_entities( ws, edges );MB_CHK_SET_ERR( rval, "Can't add edges to partition set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
430 [ # # ][ # # ]: 0 : rval = mb->add_entities( ws, faces );MB_CHK_SET_ERR( rval, "Can't add faces to partition set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
431 : 0 : }
432 : : }
433 : : #ifdef MOAB_HAVE_MPI
434 [ - + ]: 1 : if( size > 1 )
435 : : {
436 : :
437 : : // rval = mb->create_meshset(MESHSET_SET, mesh_set);MB_CHK_SET_ERR(rval, "Can't create new
438 : : // set");
439 : :
440 [ # # ]: 0 : if( parmerge )
441 : : {
442 [ # # ]: 0 : ParallelMergeMesh pm( pc, 0.00001 );
443 [ # # ][ # # ]: 0 : rval = pm.merge();MB_CHK_SET_ERR( rval, "Can't resolve shared ents" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
444 [ # # ]: 0 : if( 0 == rank )
445 : : {
446 [ # # ][ # # ]: 0 : std::cout << "parallel mesh merge: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
[ # # ][ # # ]
447 [ # # ]: 0 : tt = clock();
448 : 0 : }
449 : : }
450 : : else
451 : : {
452 [ # # ][ # # ]: 0 : rval = pc->resolve_shared_ents( cset, -1, -1, &new_id_tag );MB_CHK_SET_ERR( rval, "Can't resolve shared ents" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
453 : :
454 [ # # ]: 0 : if( 0 == rank )
455 : : {
456 [ # # ][ # # ]: 0 : std::cout << "resolve shared entities: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds"
[ # # ]
457 [ # # ]: 0 : << endl;
458 : 0 : tt = clock();
459 : : }
460 : : }
461 [ # # ]: 0 : if( !keep_skins )
462 : : { // Default is to delete the 1- and 2-dimensional entities
463 : : // Delete all quads and edges
464 [ # # ]: 0 : Range toDelete;
465 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( cset, 1, toDelete );MB_CHK_SET_ERR( rval, "Can't get edges" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
466 : :
467 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( cset, 2, toDelete );MB_CHK_SET_ERR( rval, "Can't get faces" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
468 : :
469 [ # # ][ # # ]: 0 : rval = pc->delete_entities( toDelete );MB_CHK_SET_ERR( rval, "Can't delete entities" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
470 [ # # ][ # # ]: 0 : rval = mb->remove_entities( cset, toDelete );MB_CHK_SET_ERR( rval, "Can't remove entities from base set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
471 [ # # ]: 0 : if( 0 == rank )
472 : : {
473 : :
474 [ # # ]: 0 : std::cout << "delete edges and faces \n";
475 [ # # ]: 0 : toDelete.print( std::cout );
476 : :
477 [ # # ][ # # ]: 0 : std::cout << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
[ # # ]
478 [ # # ]: 0 : tt = clock();
479 : 0 : }
480 : : }
481 : : // do some ghosting if required
482 [ # # ]: 0 : if( GL > 0 )
483 : : {
484 : : rval = pc->exchange_ghost_cells( 3, // int ghost_dim
485 : : 0, // int bridge_dim
486 : : GL, // int num_layers
487 : : 0, // int addl_ents
488 [ # # ][ # # ]: 0 : true );MB_CHK_ERR( rval ); // bool store_remote_handles
[ # # ][ # # ]
489 [ # # ]: 0 : if( 0 == rank )
490 : : {
491 [ # # ][ # # ]: 0 : std::cout << "exchange " << GL << " ghost layer(s) :" << ( clock() - tt ) / (double)CLOCKS_PER_SEC
[ # # ][ # # ]
492 [ # # ][ # # ]: 0 : << " seconds" << endl;
493 : 0 : tt = clock();
494 : : }
495 : : }
496 : : }
497 : : #endif
498 : :
499 [ + - ]: 1 : if( !parmerge )
500 : : {
501 [ + - ][ - + ]: 1 : rval = mb->tag_delete( new_id_tag );MB_CHK_SET_ERR( rval, "Can't delete new ID tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
502 : : }
503 : :
504 : 2 : return MB_SUCCESS;
505 : : }
506 : :
507 [ + - ][ + - ]: 4 : } /* namespace moab */
|
