1 #include <petsc/private/dmmbimpl.h> /*I "petscdmmoab.h" I*/
2 #include <petsc/private/vecimpl.h>
3
4 #include <petscdmmoab.h>
5 #include <MBTagConventions.hpp>
6 #include <moab/ReadUtilIface.hpp>
7 #include <moab/MergeMesh.hpp>
8 #include <moab/CN.hpp>
9
10 typedef struct {
11 // options
12 PetscInt A, B, C, M, N, K, dim;
13 PetscInt blockSizeVertexXYZ[3]; // Number of element blocks per partition
14 PetscInt blockSizeElementXYZ[3];
15 PetscReal xyzbounds[6]; // the physical size of the domain
16 bool newMergeMethod, keep_skins, simplex, adjEnts;
17
18 // compute params
19 PetscReal dx, dy, dz;
20 PetscInt NX, NY, NZ, nex, ney, nez;
21 PetscInt q, xstride, ystride, zstride;
22 PetscBool usrxyzgrid, usrprocgrid, usrrefgrid;
23 PetscInt fraction, remainder, cumfraction;
24 PetscLogEvent generateMesh, generateElements, generateVertices, parResolve;
25
26 } DMMoabMeshGeneratorCtx;
27
28
DMMoab_SetTensorElementConnectivity_Private(DMMoabMeshGeneratorCtx & genCtx,PetscInt offset,PetscInt corner,std::vector<PetscInt> & subent_conn,moab::EntityHandle * connectivity)29 PetscInt DMMoab_SetTensorElementConnectivity_Private(DMMoabMeshGeneratorCtx& genCtx, PetscInt offset, PetscInt corner, std::vector<PetscInt>& subent_conn, moab::EntityHandle *connectivity)
30 {
31 switch (genCtx.dim) {
32 case 1:
33 subent_conn.resize(2);
34 moab::CN::SubEntityVertexIndices(moab::MBEDGE, 1, 0, subent_conn.data());
35 connectivity[offset + subent_conn[0]] = corner;
36 connectivity[offset + subent_conn[1]] = corner + 1;
37 break;
38 case 2:
39 subent_conn.resize(4);
40 moab::CN::SubEntityVertexIndices(moab::MBQUAD, 2, 0, subent_conn.data());
41 connectivity[offset + subent_conn[0]] = corner;
42 connectivity[offset + subent_conn[1]] = corner + 1;
43 connectivity[offset + subent_conn[2]] = corner + 1 + genCtx.ystride;
44 connectivity[offset + subent_conn[3]] = corner + genCtx.ystride;
45 break;
46 case 3:
47 default:
48 subent_conn.resize(8);
49 moab::CN::SubEntityVertexIndices(moab::MBHEX, 3, 0, subent_conn.data());
50 connectivity[offset + subent_conn[0]] = corner;
51 connectivity[offset + subent_conn[1]] = corner + 1;
52 connectivity[offset + subent_conn[2]] = corner + 1 + genCtx.ystride;
53 connectivity[offset + subent_conn[3]] = corner + genCtx.ystride;
54 connectivity[offset + subent_conn[4]] = corner + genCtx.zstride;
55 connectivity[offset + subent_conn[5]] = corner + 1 + genCtx.zstride;
56 connectivity[offset + subent_conn[6]] = corner + 1 + genCtx.ystride + genCtx.zstride;
57 connectivity[offset + subent_conn[7]] = corner + genCtx.ystride + genCtx.zstride;
58 break;
59 }
60 return subent_conn.size();
61 }
62
63
DMMoab_SetSimplexElementConnectivity_Private(DMMoabMeshGeneratorCtx & genCtx,PetscInt subelem,PetscInt offset,PetscInt corner,std::vector<PetscInt> & subent_conn,moab::EntityHandle * connectivity)64 PetscInt DMMoab_SetSimplexElementConnectivity_Private(DMMoabMeshGeneratorCtx& genCtx, PetscInt subelem, PetscInt offset, PetscInt corner, std::vector<PetscInt>& subent_conn, moab::EntityHandle *connectivity)
65 {
66 PetscInt A, B, C, D, E, F, G, H, M;
67 const PetscInt trigen_opts = 1; /* 1 - Aligned diagonally to right, 2 - Aligned diagonally to left, 3 - 4 elements per quad */
68 A = corner;
69 B = corner + 1;
70 switch (genCtx.dim) {
71 case 1:
72 subent_conn.resize(2); /* only linear EDGE supported now */
73 moab::CN::SubEntityVertexIndices(moab::MBEDGE, 1, 0, subent_conn.data());
74 connectivity[offset + subent_conn[0]] = A;
75 connectivity[offset + subent_conn[1]] = B;
76 break;
77 case 2:
78 C = corner + 1 + genCtx.ystride;
79 D = corner + genCtx.ystride;
80 M = corner + 0.5; /* technically -- need to modify vertex generation */
81 subent_conn.resize(3); /* only linear TRI supported */
82 moab::CN::SubEntityVertexIndices(moab::MBTRI, 2, 0, subent_conn.data());
83 if (trigen_opts == 1) {
84 if (subelem) { /* 0 1 2 of a QUAD */
85 connectivity[offset + subent_conn[0]] = B;
86 connectivity[offset + subent_conn[1]] = C;
87 connectivity[offset + subent_conn[2]] = A;
88 }
89 else { /* 2 3 0 of a QUAD */
90 connectivity[offset + subent_conn[0]] = D;
91 connectivity[offset + subent_conn[1]] = A;
92 connectivity[offset + subent_conn[2]] = C;
93 }
94 }
95 else if (trigen_opts == 2) {
96 if (subelem) { /* 0 1 2 of a QUAD */
97 connectivity[offset + subent_conn[0]] = A;
98 connectivity[offset + subent_conn[1]] = B;
99 connectivity[offset + subent_conn[2]] = D;
100 }
101 else { /* 2 3 0 of a QUAD */
102 connectivity[offset + subent_conn[0]] = C;
103 connectivity[offset + subent_conn[1]] = D;
104 connectivity[offset + subent_conn[2]] = B;
105 }
106 }
107 else {
108 switch (subelem) { /* 0 1 2 of a QUAD */
109 case 0:
110 connectivity[offset + subent_conn[0]] = A;
111 connectivity[offset + subent_conn[1]] = B;
112 connectivity[offset + subent_conn[2]] = M;
113 break;
114 case 1:
115 connectivity[offset + subent_conn[0]] = B;
116 connectivity[offset + subent_conn[1]] = C;
117 connectivity[offset + subent_conn[2]] = M;
118 break;
119 case 2:
120 connectivity[offset + subent_conn[0]] = C;
121 connectivity[offset + subent_conn[1]] = D;
122 connectivity[offset + subent_conn[2]] = M;
123 break;
124 case 3:
125 connectivity[offset + subent_conn[0]] = D;
126 connectivity[offset + subent_conn[1]] = A;
127 connectivity[offset + subent_conn[2]] = M;
128 break;
129 }
130 }
131 break;
132 case 3:
133 default:
134 C = corner + 1 + genCtx.ystride;
135 D = corner + genCtx.ystride;
136 E = corner + genCtx.zstride;
137 F = corner + 1 + genCtx.zstride;
138 G = corner + 1 + genCtx.ystride + genCtx.zstride;
139 H = corner + genCtx.ystride + genCtx.zstride;
140 subent_conn.resize(4); /* only linear TET supported */
141 moab::CN::SubEntityVertexIndices(moab::MBTET, 3, 0, subent_conn.data());
142 switch (subelem) {
143 case 0: /* 4 3 7 6 of a HEX */
144 connectivity[offset + subent_conn[0]] = E;
145 connectivity[offset + subent_conn[1]] = D;
146 connectivity[offset + subent_conn[2]] = H;
147 connectivity[offset + subent_conn[3]] = G;
148 break;
149 case 1: /* 0 1 2 5 of a HEX */
150 connectivity[offset + subent_conn[0]] = A;
151 connectivity[offset + subent_conn[1]] = B;
152 connectivity[offset + subent_conn[2]] = C;
153 connectivity[offset + subent_conn[3]] = F;
154 break;
155 case 2: /* 0 3 4 5 of a HEX */
156 connectivity[offset + subent_conn[0]] = A;
157 connectivity[offset + subent_conn[1]] = D;
158 connectivity[offset + subent_conn[2]] = E;
159 connectivity[offset + subent_conn[3]] = F;
160 break;
161 case 3: /* 2 6 3 5 of a HEX */
162 connectivity[offset + subent_conn[0]] = C;
163 connectivity[offset + subent_conn[1]] = G;
164 connectivity[offset + subent_conn[2]] = D;
165 connectivity[offset + subent_conn[3]] = F;
166 break;
167 case 4: /* 0 2 3 5 of a HEX */
168 connectivity[offset + subent_conn[0]] = A;
169 connectivity[offset + subent_conn[1]] = C;
170 connectivity[offset + subent_conn[2]] = D;
171 connectivity[offset + subent_conn[3]] = F;
172 break;
173 case 5: /* 3 6 4 5 of a HEX */
174 connectivity[offset + subent_conn[0]] = D;
175 connectivity[offset + subent_conn[1]] = G;
176 connectivity[offset + subent_conn[2]] = E;
177 connectivity[offset + subent_conn[3]] = F;
178 break;
179 }
180 break;
181 }
182 return subent_conn.size();
183 }
184
185
DMMoab_SetElementConnectivity_Private(DMMoabMeshGeneratorCtx & genCtx,PetscInt offset,PetscInt corner,moab::EntityHandle * connectivity)186 std::pair<PetscInt, PetscInt> DMMoab_SetElementConnectivity_Private(DMMoabMeshGeneratorCtx& genCtx, PetscInt offset, PetscInt corner, moab::EntityHandle *connectivity)
187 {
188 PetscInt vcount = 0;
189 PetscInt simplices_per_tensor[4] = {0, 1, 2, 6};
190 std::vector<PetscInt> subent_conn; /* only linear edge, tri, tet supported now */
191 subent_conn.reserve(27);
192 PetscInt m, subelem;
193 if (genCtx.simplex) {
194 subelem = simplices_per_tensor[genCtx.dim];
195 for (m = 0; m < subelem; m++) {
196 vcount = DMMoab_SetSimplexElementConnectivity_Private(genCtx, m, offset, corner, subent_conn, connectivity);
197 offset += vcount;
198 }
199 }
200 else {
201 subelem = 1;
202 vcount = DMMoab_SetTensorElementConnectivity_Private(genCtx, offset, corner, subent_conn, connectivity);
203 }
204 return std::pair<PetscInt, PetscInt>(vcount * subelem, subelem);
205 }
206
207
DMMoab_GenerateVertices_Private(moab::Interface * mbImpl,moab::ReadUtilIface * iface,DMMoabMeshGeneratorCtx & genCtx,PetscInt m,PetscInt n,PetscInt k,PetscInt a,PetscInt b,PetscInt c,moab::Tag & global_id_tag,moab::EntityHandle & startv,moab::Range & uverts)208 PetscErrorCode DMMoab_GenerateVertices_Private(moab::Interface *mbImpl, moab::ReadUtilIface *iface, DMMoabMeshGeneratorCtx& genCtx, PetscInt m, PetscInt n, PetscInt k,
209 PetscInt a, PetscInt b, PetscInt c, moab::Tag& global_id_tag, moab::EntityHandle& startv, moab::Range& uverts)
210 {
211 PetscInt x, y, z, ix, nnodes;
212 PetscErrorCode ierr;
213 PetscInt ii, jj, kk;
214 std::vector<PetscReal*> arrays;
215 PetscInt* gids;
216 moab::ErrorCode merr;
217
218 PetscFunctionBegin;
219 /* we will generate (q*block+1)^3 vertices, and block^3 hexas; q is 1 for linear, 2 for quadratic
220 * the global id of the vertices will come from m, n, k, a, b, c
221 * x will vary from m*A*q*block + a*q*block to m*A*q*block+(a+1)*q*block etc.
222 */
223 nnodes = genCtx.blockSizeVertexXYZ[0] * (genCtx.dim > 1 ? genCtx.blockSizeVertexXYZ[1] * (genCtx.dim > 2 ? genCtx.blockSizeVertexXYZ[2] : 1) : 1);
224 ierr = PetscMalloc1(nnodes, &gids);CHKERRQ(ierr);
225
226 merr = iface->get_node_coords(3, nnodes, 0, startv, arrays);MBERR("Can't get node coords.", merr);
227
228 /* will start with the lower corner: */
229 /* x = ( m * genCtx.A + a) * genCtx.q * genCtx.blockSizeElementXYZ[0]; */
230 /* y = ( n * genCtx.B + b) * genCtx.q * genCtx.blockSizeElementXYZ[1]; */
231 /* z = ( k * genCtx.C + c) * genCtx.q * genCtx.blockSizeElementXYZ[2]; */
232
233 x = ( m * genCtx.A + a) * genCtx.q;
234 y = ( n * genCtx.B + b) * genCtx.q;
235 z = ( k * genCtx.C + c) * genCtx.q;
236 PetscInfo3(NULL, "Starting offset for coordinates := %d, %d, %d\n", x, y, z);
237 ix = 0;
238 moab::Range verts(startv, startv + nnodes - 1);
239 for (kk = 0; kk < (genCtx.dim > 2 ? genCtx.blockSizeVertexXYZ[2] : 1); kk++) {
240 for (jj = 0; jj < (genCtx.dim > 1 ? genCtx.blockSizeVertexXYZ[1] : 1); jj++) {
241 for (ii = 0; ii < genCtx.blockSizeVertexXYZ[0]; ii++, ix++) {
242 /* set coordinates for the vertices */
243 arrays[0][ix] = (x + ii) * genCtx.dx + genCtx.xyzbounds[0];
244 arrays[1][ix] = (y + jj) * genCtx.dy + genCtx.xyzbounds[2];
245 arrays[2][ix] = (z + kk) * genCtx.dz + genCtx.xyzbounds[4];
246 PetscInfo3(NULL, "Creating vertex with coordinates := %f, %f, %f\n", arrays[0][ix], arrays[1][ix], arrays[2][ix]);
247
248 /* If we want to set some tags on the vertices -> use the following entity handle definition:
249 moab::EntityHandle v = startv + ix;
250 */
251 /* compute the global ID for vertex */
252 gids[ix] = 1 + (x + ii) + (y + jj) * genCtx.NX + (z + kk) * (genCtx.NX * genCtx.NY);
253 }
254 }
255 }
256 /* set global ID data on vertices */
257 mbImpl->tag_set_data(global_id_tag, verts, &gids[0]);
258 verts.swap(uverts);
259 ierr = PetscFree(gids);CHKERRQ(ierr);
260 PetscFunctionReturn(0);
261 }
262
DMMoab_GenerateElements_Private(moab::Interface * mbImpl,moab::ReadUtilIface * iface,DMMoabMeshGeneratorCtx & genCtx,PetscInt m,PetscInt n,PetscInt k,PetscInt a,PetscInt b,PetscInt c,moab::Tag & global_id_tag,moab::EntityHandle startv,moab::Range & cells)263 PetscErrorCode DMMoab_GenerateElements_Private(moab::Interface* mbImpl, moab::ReadUtilIface* iface, DMMoabMeshGeneratorCtx& genCtx, PetscInt m, PetscInt n, PetscInt k,
264 PetscInt a, PetscInt b, PetscInt c, moab::Tag& global_id_tag, moab::EntityHandle startv, moab::Range& cells)
265 {
266 moab::ErrorCode merr;
267 PetscInt ix, ie, xe, ye, ze;
268 PetscInt ii, jj, kk, nvperelem;
269 PetscInt simplices_per_tensor[4] = {0, 1, 2, 6};
270 PetscInt ntensorelems = genCtx.blockSizeElementXYZ[0] * (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] * (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1) : 1); /*pow(genCtx.blockSizeElement,genCtx.dim);*/
271 PetscInt nelems = ntensorelems;
272 moab::EntityHandle starte; /* connectivity */
273 moab::EntityHandle* conn;
274
275 PetscFunctionBegin;
276 switch (genCtx.dim) {
277 case 1:
278 nvperelem = 2;
279 merr = iface->get_element_connect(nelems, 2, moab::MBEDGE, 0, starte, conn);MBERR("Can't get EDGE2 element connectivity.", merr);
280 break;
281 case 2:
282 if (genCtx.simplex) {
283 nvperelem = 3;
284 nelems = ntensorelems * simplices_per_tensor[genCtx.dim];
285 merr = iface->get_element_connect(nelems, 3, moab::MBTRI, 0, starte, conn);MBERR("Can't get TRI3 element connectivity.", merr);
286 }
287 else {
288 nvperelem = 4;
289 merr = iface->get_element_connect(nelems, 4, moab::MBQUAD, 0, starte, conn);MBERR("Can't get QUAD4 element connectivity.", merr);
290 }
291 break;
292 case 3:
293 default:
294 if (genCtx.simplex) {
295 nvperelem = 4;
296 nelems = ntensorelems * simplices_per_tensor[genCtx.dim];
297 merr = iface->get_element_connect(nelems, 4, moab::MBTET, 0, starte, conn);MBERR("Can't get TET4 element connectivity.", merr);
298 }
299 else {
300 nvperelem = 8;
301 merr = iface->get_element_connect(nelems, 8, moab::MBHEX, 0, starte, conn);MBERR("Can't get HEX8 element connectivity.", merr);
302 }
303 break;
304 }
305
306 ix = ie = 0; /* index now in the elements, for global ids */
307
308 /* create a temporary range to store local element handles */
309 moab::Range tmp(starte, starte + nelems - 1);
310 std::vector<PetscInt> gids(nelems);
311
312 /* identify the elements at the lower corner, for their global ids */
313 xe = m * genCtx.A * genCtx.blockSizeElementXYZ[0] + a * genCtx.blockSizeElementXYZ[0];
314 ye = (genCtx.dim > 1 ? n * genCtx.B * genCtx.blockSizeElementXYZ[1] + b * genCtx.blockSizeElementXYZ[1] : 0);
315 ze = (genCtx.dim > 2 ? k * genCtx.C * genCtx.blockSizeElementXYZ[2] + c * genCtx.blockSizeElementXYZ[2] : 0);
316
317 /* create owned elements requested by genCtx */
318 for (kk = 0; kk < (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1); kk++) {
319 for (jj = 0; jj < (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] : 1); jj++) {
320 for (ii = 0; ii < genCtx.blockSizeElementXYZ[0]; ii++) {
321
322 moab::EntityHandle corner = startv + genCtx.q * ii + genCtx.q * jj * genCtx.ystride + genCtx.q * kk * genCtx.zstride;
323
324 std::pair<PetscInt, PetscInt> entoffset = DMMoab_SetElementConnectivity_Private(genCtx, ix, corner, conn);
325
326 for (PetscInt j = 0; j < entoffset.second; j++) {
327 /* The entity handle for the particular element -> if we want to set some tags is
328 moab::EntityHandle eh = starte + ie + j;
329 */
330 gids[ie + j] = 1 + ((xe + ii) + (ye + jj) * genCtx.nex + (ze + kk) * (genCtx.nex * genCtx.ney));
331 /* gids[ie+j] = ie + j + ((xe + ii) + (ye + jj) * genCtx.nex + (ze + kk) * (genCtx.nex * genCtx.ney)); */
332 /* gids[ie+j] = 1 + ie; */
333 /* ie++; */
334 }
335
336 ix += entoffset.first;
337 ie += entoffset.second;
338 }
339 }
340 }
341 if (genCtx.adjEnts) { /* we need to update adjacencies now, because some elements are new */
342 merr = iface->update_adjacencies(starte, nelems, nvperelem, conn);MBERR("Can't update adjacencies", merr);
343 }
344 tmp.swap(cells);
345 merr = mbImpl->tag_set_data(global_id_tag, cells, &gids[0]);MBERR("Can't set global ids to elements.", merr);
346 PetscFunctionReturn(0);
347 }
348
DMMBUtil_InitializeOptions(DMMoabMeshGeneratorCtx & genCtx,PetscInt dim,PetscBool simplex,PetscInt rank,PetscInt nprocs,const PetscReal * bounds,PetscInt nelems)349 PetscErrorCode DMMBUtil_InitializeOptions(DMMoabMeshGeneratorCtx& genCtx, PetscInt dim, PetscBool simplex, PetscInt rank, PetscInt nprocs, const PetscReal* bounds, PetscInt nelems)
350 {
351 PetscFunctionBegin;
352 /* Initialize all genCtx data */
353 genCtx.dim = dim;
354 genCtx.simplex = simplex;
355 genCtx.newMergeMethod = genCtx.keep_skins = genCtx.adjEnts = true;
356 /* determine other global quantities for the mesh used for nodes increments */
357 genCtx.q = 1;
358 genCtx.fraction = genCtx.remainder = genCtx.cumfraction = 0;
359
360 if (!genCtx.usrxyzgrid) { /* not overridden by genCtx - assume nele equally and that genCtx wants a uniform cube mesh */
361
362 genCtx.fraction = nelems / nprocs; /* partition only by the largest dimension */
363 genCtx.remainder = nelems % nprocs; /* remainder after partition which gets evenly distributed by round-robin */
364 genCtx.cumfraction = (rank > 0 ? (genCtx.fraction) * (rank) + (rank - 1 < genCtx.remainder ? rank : genCtx.remainder) : 0);
365 if (rank < genCtx.remainder) /* This process gets "fraction+1" elements */
366 genCtx.fraction++;
367
368 PetscInfo3(NULL, "Fraction = %D, Remainder = %D, Cumulative fraction = %D\n", genCtx.fraction, genCtx.remainder, genCtx.cumfraction);
369 switch (genCtx.dim) {
370 case 1:
371 genCtx.blockSizeElementXYZ[0] = genCtx.fraction;
372 genCtx.blockSizeElementXYZ[1] = 1;
373 genCtx.blockSizeElementXYZ[2] = 1;
374 break;
375 case 2:
376 genCtx.blockSizeElementXYZ[0] = nelems;
377 genCtx.blockSizeElementXYZ[1] = genCtx.fraction;
378 genCtx.blockSizeElementXYZ[2] = 1;
379 break;
380 case 3:
381 default:
382 genCtx.blockSizeElementXYZ[0] = nelems;
383 genCtx.blockSizeElementXYZ[1] = nelems;
384 genCtx.blockSizeElementXYZ[2] = genCtx.fraction;
385 break;
386 }
387 }
388
389 /* partition only by the largest dimension */
390 /* Total number of local elements := genCtx.blockSizeElementXYZ[0]*(genCtx.dim>1? genCtx.blockSizeElementXYZ[1]*(genCtx.dim>2 ? genCtx.blockSizeElementXYZ[2]:1) :1); */
391 if (bounds) {
392 for (PetscInt i = 0; i < 6; i++)
393 genCtx.xyzbounds[i] = bounds[i];
394 }
395 else {
396 genCtx.xyzbounds[0] = genCtx.xyzbounds[2] = genCtx.xyzbounds[4] = 0.0;
397 genCtx.xyzbounds[1] = genCtx.xyzbounds[3] = genCtx.xyzbounds[5] = 1.0;
398 }
399
400 if (!genCtx.usrprocgrid) {
401 switch (genCtx.dim) {
402 case 1:
403 genCtx.M = nprocs;
404 genCtx.N = genCtx.K = 1;
405 break;
406 case 2:
407 genCtx.N = nprocs;
408 genCtx.M = genCtx.K = 1;
409 break;
410 default:
411 genCtx.K = nprocs;
412 genCtx.M = genCtx.N = 1;
413 break;
414 }
415 }
416
417 if (!genCtx.usrrefgrid) {
418 genCtx.A = genCtx.B = genCtx.C = 1;
419 }
420
421 /* more default values */
422 genCtx.nex = genCtx.ney = genCtx.nez = 0;
423 genCtx.xstride = genCtx.ystride = genCtx.zstride = 0;
424 genCtx.NX = genCtx.NY = genCtx.NZ = 0;
425 genCtx.nex = genCtx.ney = genCtx.nez = 0;
426 genCtx.blockSizeVertexXYZ[0] = genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 1;
427
428 switch (genCtx.dim) {
429 case 3:
430 genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1;
431 genCtx.blockSizeVertexXYZ[1] = genCtx.q * genCtx.blockSizeElementXYZ[1] + 1;
432 genCtx.blockSizeVertexXYZ[2] = genCtx.q * genCtx.blockSizeElementXYZ[2] + 1;
433
434 genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0]; /* number of elements in x direction, used for global id on element */
435 genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (nelems * genCtx.q); /* distance between 2 nodes in x direction */
436 genCtx.NX = (genCtx.q * genCtx.nex + 1);
437 genCtx.xstride = 1;
438 genCtx.ney = genCtx.N * genCtx.B * genCtx.blockSizeElementXYZ[1]; /* number of elements in y direction .... */
439 genCtx.dy = (genCtx.xyzbounds[3] - genCtx.xyzbounds[2]) / (nelems * genCtx.q); /* distance between 2 nodes in y direction */
440 genCtx.NY = (genCtx.q * genCtx.ney + 1);
441 genCtx.ystride = genCtx.blockSizeVertexXYZ[0];
442 genCtx.nez = genCtx.K * genCtx.C * genCtx.blockSizeElementXYZ[2]; /* number of elements in z direction .... */
443 genCtx.dz = (genCtx.xyzbounds[5] - genCtx.xyzbounds[4]) / (nelems * genCtx.q); /* distance between 2 nodes in z direction */
444 genCtx.NZ = (genCtx.q * genCtx.nez + 1);
445 genCtx.zstride = genCtx.blockSizeVertexXYZ[0] * genCtx.blockSizeVertexXYZ[1];
446 break;
447 case 2:
448 genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1;
449 genCtx.blockSizeVertexXYZ[1] = genCtx.q * genCtx.blockSizeElementXYZ[1] + 1;
450 genCtx.blockSizeVertexXYZ[2] = 0;
451
452 genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0]; /* number of elements in x direction, used for global id on element */
453 genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (genCtx.nex * genCtx.q); /* distance between 2 nodes in x direction */
454 genCtx.NX = (genCtx.q * genCtx.nex + 1);
455 genCtx.xstride = 1;
456 genCtx.ney = genCtx.N * genCtx.B * genCtx.blockSizeElementXYZ[1]; /* number of elements in y direction .... */
457 genCtx.dy = (genCtx.xyzbounds[3] - genCtx.xyzbounds[2]) / (nelems * genCtx.q); /* distance between 2 nodes in y direction */
458 genCtx.NY = (genCtx.q * genCtx.ney + 1);
459 genCtx.ystride = genCtx.blockSizeVertexXYZ[0];
460 break;
461 case 1:
462 genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 0;
463 genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1;
464
465 genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0]; /* number of elements in x direction, used for global id on element */
466 genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (nelems * genCtx.q); /* distance between 2 nodes in x direction */
467 genCtx.NX = (genCtx.q * genCtx.nex + 1);
468 genCtx.xstride = 1;
469 break;
470 }
471
472 /* Lets check for some valid input */
473 if (genCtx.dim < 1 || genCtx.dim > 3) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid topological dimension specified: %d.\n", genCtx.dim);
474 if (genCtx.M * genCtx.N * genCtx.K != nprocs) SETERRQ4(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid [m, n, k] data: %d, %d, %d. Product must be equal to global size = %d.\n", genCtx.M, genCtx.N, genCtx.K, nprocs);
475 /* validate the bounds data */
476 if (genCtx.xyzbounds[0] >= genCtx.xyzbounds[1]) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "X-dim: Left boundary cannot be greater than right. [%G >= %G]\n", genCtx.xyzbounds[0], genCtx.xyzbounds[1]);
477 if (genCtx.dim > 1 && (genCtx.xyzbounds[2] >= genCtx.xyzbounds[3])) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Y-dim: Left boundary cannot be greater than right. [%G >= %G]\n", genCtx.xyzbounds[2], genCtx.xyzbounds[3]);
478 if (genCtx.dim > 2 && (genCtx.xyzbounds[4] >= genCtx.xyzbounds[5])) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Z-dim: Left boundary cannot be greater than right. [%G >= %G]\n", genCtx.xyzbounds[4], genCtx.xyzbounds[5]);
479
480 PetscInfo3(NULL, "Local elements:= %d, %d, %d\n", genCtx.blockSizeElementXYZ[0], genCtx.blockSizeElementXYZ[1], genCtx.blockSizeElementXYZ[2]);
481 PetscInfo3(NULL, "Local vertices:= %d, %d, %d\n", genCtx.blockSizeVertexXYZ[0], genCtx.blockSizeVertexXYZ[1], genCtx.blockSizeVertexXYZ[2]);
482 PetscInfo3(NULL, "Local blocks/processors := %d, %d, %d\n", genCtx.A, genCtx.B, genCtx.C);
483 PetscInfo3(NULL, "Local processors := %d, %d, %d\n", genCtx.M, genCtx.N, genCtx.K);
484 PetscInfo3(NULL, "Local nexyz:= %d, %d, %d\n", genCtx.nex, genCtx.ney, genCtx.nez);
485 PetscInfo3(NULL, "Local delxyz:= %g, %g, %g\n", genCtx.dx, genCtx.dy, genCtx.dz);
486 PetscInfo3(NULL, "Local strides:= %d, %d, %d\n", genCtx.xstride, genCtx.ystride, genCtx.zstride);
487 PetscFunctionReturn(0);
488 }
489
490 /*@C
491 DMMoabCreateBoxMesh - Creates a mesh on the tensor product (box) of intervals with genCtx specified bounds.
492
493 Collective
494
495 Input Parameters:
496 + comm - The communicator for the DM object
497 . dim - The spatial dimension
498 . bounds - The bounds of the box specified with [x-left, x-right, y-bottom, y-top, z-bottom, z-top] depending on the spatial dimension
499 . nele - The number of discrete elements in each direction
500 - user_nghost - The number of ghosted layers needed in the partitioned mesh
501
502 Output Parameter:
503 . dm - The DM object
504
505 Level: beginner
506
507 .seealso: DMSetType(), DMCreate(), DMMoabLoadFromFile()
508 @*/
DMMoabCreateBoxMesh(MPI_Comm comm,PetscInt dim,PetscBool useSimplex,const PetscReal * bounds,PetscInt nele,PetscInt nghost,DM * dm)509 PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool useSimplex, const PetscReal* bounds, PetscInt nele, PetscInt nghost, DM *dm)
510 {
511 PetscErrorCode ierr;
512 moab::ErrorCode merr;
513 PetscInt a, b, c, n, global_size, global_rank;
514 DM_Moab *dmmoab;
515 moab::Interface *mbImpl;
516 #ifdef MOAB_HAVE_MPI
517 moab::ParallelComm *pcomm;
518 #endif
519 moab::ReadUtilIface *readMeshIface;
520 moab::Range verts, cells, edges, faces, adj, dim3, dim2;
521 DMMoabMeshGeneratorCtx genCtx;
522 const PetscInt npts = nele + 1; /* Number of points in every dimension */
523
524 moab::Tag global_id_tag, part_tag, geom_tag, mat_tag, dir_tag, neu_tag;
525 moab::Range ownedvtx, ownedelms, localvtxs, localelms;
526 moab::EntityHandle regionset;
527 PetscInt ml = 0, nl = 0, kl = 0;
528
529 PetscFunctionBegin;
530 if (dim < 1 || dim > 3) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension argument for mesh: dim=[1,3].\n");
531
532 ierr = PetscLogEventRegister("GenerateMesh", DM_CLASSID, &genCtx.generateMesh);CHKERRQ(ierr);
533 ierr = PetscLogEventRegister("AddVertices", DM_CLASSID, &genCtx.generateVertices);CHKERRQ(ierr);
534 ierr = PetscLogEventRegister("AddElements", DM_CLASSID, &genCtx.generateElements);CHKERRQ(ierr);
535 ierr = PetscLogEventRegister("ParResolve", DM_CLASSID, &genCtx.parResolve);CHKERRQ(ierr);
536 ierr = PetscLogEventBegin(genCtx.generateMesh, 0, 0, 0, 0);CHKERRQ(ierr);
537 ierr = MPI_Comm_size(comm, &global_size);CHKERRQ(ierr);
538 /* total number of vertices in all dimensions */
539 n = pow(npts, dim);
540
541 /* do some error checking */
542 if (n < 2) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be >= 2.\n");
543 if (global_size > n) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of processors must be less than or equal to number of elements.\n");
544 if (nghost < 0) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of ghost layers cannot be negative.\n");
545
546 /* Create the basic DMMoab object and keep the default parameters created by DM impls */
547 ierr = DMMoabCreateMoab(comm, NULL, NULL, NULL, dm);CHKERRQ(ierr);
548
549 /* get all the necessary handles from the private DM object */
550 dmmoab = (DM_Moab*)(*dm)->data;
551 mbImpl = dmmoab->mbiface;
552 #ifdef MOAB_HAVE_MPI
553 pcomm = dmmoab->pcomm;
554 global_rank = pcomm->rank();
555 #else
556 global_rank = 0;
557 global_size = 1;
558 #endif
559 global_id_tag = dmmoab->ltog_tag;
560 dmmoab->dim = dim;
561 dmmoab->nghostrings = nghost;
562 dmmoab->refct = 1;
563
564 /* create a file set to associate all entities in current mesh */
565 merr = mbImpl->create_meshset(moab::MESHSET_SET, dmmoab->fileset);MBERR("Creating file set failed", merr);
566
567 /* No errors yet; proceed with building the mesh */
568 merr = mbImpl->query_interface(readMeshIface);MBERRNM(merr);
569
570 genCtx.M = genCtx.N = genCtx.K = 1;
571 genCtx.A = genCtx.B = genCtx.C = 1;
572 genCtx.blockSizeElementXYZ[0] = 0;
573 genCtx.blockSizeElementXYZ[1] = 0;
574 genCtx.blockSizeElementXYZ[2] = 0;
575
576 ierr = PetscOptionsBegin(comm, "", "DMMoab Creation Options", "DMMOAB");CHKERRQ(ierr);
577 /* Handle DMMoab spatial resolution */
578 ierr = PetscOptionsInt("-dmb_grid_x", "Number of grid points in x direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[0], &genCtx.blockSizeElementXYZ[0], &genCtx.usrxyzgrid);CHKERRQ(ierr);
579 if (dim > 1) {ierr = PetscOptionsInt("-dmb_grid_y", "Number of grid points in y direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[1], &genCtx.blockSizeElementXYZ[1], &genCtx.usrxyzgrid);CHKERRQ(ierr);}
580 if (dim > 2) {ierr = PetscOptionsInt("-dmb_grid_z", "Number of grid points in z direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[2], &genCtx.blockSizeElementXYZ[2], &genCtx.usrxyzgrid);CHKERRQ(ierr);}
581
582 /* Handle DMMoab parallel distibution */
583 ierr = PetscOptionsInt("-dmb_processors_x", "Number of processors in x direction", "DMMoabSetNumProcs", genCtx.M, &genCtx.M, &genCtx.usrprocgrid);CHKERRQ(ierr);
584 if (dim > 1) {ierr = PetscOptionsInt("-dmb_processors_y", "Number of processors in y direction", "DMMoabSetNumProcs", genCtx.N, &genCtx.N, &genCtx.usrprocgrid);CHKERRQ(ierr);}
585 if (dim > 2) {ierr = PetscOptionsInt("-dmb_processors_z", "Number of processors in z direction", "DMMoabSetNumProcs", genCtx.K, &genCtx.K, &genCtx.usrprocgrid);CHKERRQ(ierr);}
586
587 /* Handle DMMoab block refinement */
588 ierr = PetscOptionsInt("-dmb_refine_x", "Number of refinement blocks in x direction", "DMMoabSetRefinement", genCtx.A, &genCtx.A, &genCtx.usrrefgrid);
589 if (dim > 1) {ierr = PetscOptionsInt("-dmb_refine_y", "Number of refinement blocks in y direction", "DMMoabSetRefinement", genCtx.B, &genCtx.B, &genCtx.usrrefgrid);CHKERRQ(ierr);}
590 if (dim > 2) {ierr = PetscOptionsInt("-dmb_refine_z", "Number of refinement blocks in z direction", "DMMoabSetRefinement", genCtx.C, &genCtx.C, &genCtx.usrrefgrid);CHKERRQ(ierr);}
591 PetscOptionsEnd();
592
593 ierr = DMMBUtil_InitializeOptions(genCtx, dim, useSimplex, global_rank, global_size, bounds, nele);CHKERRQ(ierr);
594
595 //if (nele<nprocs) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The dimensional discretization size should be greater or equal to number of processors: %D < %D",nele,nprocs);
596
597 if (genCtx.adjEnts) genCtx.keep_skins = true; /* do not delete anything - consumes more memory */
598
599 /* determine m, n, k for processor rank */
600 ml = nl = kl = 0;
601 switch (genCtx.dim) {
602 case 1:
603 ml = (genCtx.cumfraction);
604 break;
605 case 2:
606 nl = (genCtx.cumfraction);
607 break;
608 default:
609 kl = (genCtx.cumfraction) / genCtx.q / genCtx.blockSizeElementXYZ[2] / genCtx.C; //genCtx.K
610 break;
611 }
612
613 /*
614 * so there are a total of M * A * blockSizeElement elements in x direction (so M * A * blockSizeElement + 1 verts in x direction)
615 * so there are a total of N * B * blockSizeElement elements in y direction (so N * B * blockSizeElement + 1 verts in y direction)
616 * so there are a total of K * C * blockSizeElement elements in z direction (so K * C * blockSizeElement + 1 verts in z direction)
617
618 * there are ( M * A blockSizeElement) * ( N * B * blockSizeElement) * (K * C * blockSizeElement) hexas
619 * there are ( M * A * blockSizeElement + 1) * ( N * B * blockSizeElement + 1) * (K * C * blockSizeElement + 1) vertices
620 * x is the first dimension that varies
621 */
622
623 /* generate the block at (a, b, c); it will represent a partition , it will get a partition tag */
624 PetscInt dum_id = -1;
625 merr = mbImpl->tag_get_handle("GLOBAL_ID", 1, moab::MB_TYPE_INTEGER, global_id_tag);MBERR("Getting Global_ID Tag handle failed", merr);
626
627 merr = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, mat_tag);MBERR("Getting Material set Tag handle failed", merr);
628 merr = mbImpl->tag_get_handle(DIRICHLET_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, dir_tag);MBERR("Getting Dirichlet set Tag handle failed", merr);
629 merr = mbImpl->tag_get_handle(NEUMANN_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, neu_tag);MBERR("Getting Neumann set Tag handle failed", merr);
630
631 merr = mbImpl->tag_get_handle("PARALLEL_PARTITION", 1, moab::MB_TYPE_INTEGER, part_tag, moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &dum_id);MBERR("Getting Partition Tag handle failed", merr);
632
633 /* lets create some sets */
634 merr = mbImpl->tag_get_handle(GEOM_DIMENSION_TAG_NAME, 1, moab::MB_TYPE_INTEGER, geom_tag, moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &dum_id);MBERRNM(merr);
635 merr = mbImpl->create_meshset(moab::MESHSET_SET, regionset);MBERRNM(merr);
636 ierr = PetscLogEventEnd(genCtx.generateMesh, 0, 0, 0, 0);CHKERRQ(ierr);
637
638 for (a = 0; a < (genCtx.dim > 0 ? genCtx.A : genCtx.A); a++) {
639 for (b = 0; b < (genCtx.dim > 1 ? genCtx.B : 1); b++) {
640 for (c = 0; c < (genCtx.dim > 2 ? genCtx.C : 1); c++) {
641
642 moab::EntityHandle startv;
643
644 ierr = PetscLogEventBegin(genCtx.generateVertices, 0, 0, 0, 0);CHKERRQ(ierr);
645 ierr = DMMoab_GenerateVertices_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, verts);CHKERRQ(ierr);
646 ierr = PetscLogEventEnd(genCtx.generateVertices, 0, 0, 0, 0);CHKERRQ(ierr);
647
648 ierr = PetscLogEventBegin(genCtx.generateElements, 0, 0, 0, 0);CHKERRQ(ierr);
649 ierr = DMMoab_GenerateElements_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, cells);CHKERRQ(ierr);
650 ierr = PetscLogEventEnd(genCtx.generateElements, 0, 0, 0, 0);CHKERRQ(ierr);
651
652 PetscInt part_num = 0;
653 switch (genCtx.dim) {
654 case 3:
655 part_num += (c + kl * genCtx.C) * (genCtx.M * genCtx.A * genCtx.N * genCtx.B);
656 case 2:
657 part_num += (b + nl * genCtx.B) * (genCtx.M * genCtx.A);
658 case 1:
659 part_num += (a + ml * genCtx.A);
660 break;
661 }
662
663 moab::EntityHandle part_set;
664 merr = mbImpl->create_meshset(moab::MESHSET_SET, part_set);MBERR("Can't create mesh set.", merr);
665
666 merr = mbImpl->add_entities(part_set, verts);MBERR("Can't add vertices to set.", merr);
667 merr = mbImpl->add_entities(part_set, cells);MBERR("Can't add entities to set.", merr);
668 merr = mbImpl->add_entities(regionset, cells);MBERR("Can't add entities to set.", merr);
669
670 /* if needed, add all edges and faces */
671 if (genCtx.adjEnts)
672 {
673 if (genCtx.dim > 1) {
674 merr = mbImpl->get_adjacencies(cells, 1, true, edges, moab::Interface::UNION);MBERR("Can't get edges", merr);
675 merr = mbImpl->add_entities(part_set, edges);MBERR("Can't add edges to partition set.", merr);
676 }
677 if (genCtx.dim > 2) {
678 merr = mbImpl->get_adjacencies(cells, 2, true, faces, moab::Interface::UNION);MBERR("Can't get faces", merr);
679 merr = mbImpl->add_entities(part_set, faces);MBERR("Can't add faces to partition set.", merr);
680 }
681 edges.clear();
682 faces.clear();
683 }
684 verts.clear(); cells.clear();
685
686 merr = mbImpl->tag_set_data(part_tag, &part_set, 1, &part_num);MBERR("Can't set part tag on set", merr);
687 if (dmmoab->fileset) {
688 merr = mbImpl->add_parent_child(dmmoab->fileset, part_set);MBERR("Can't add part set to file set.", merr);
689 merr = mbImpl->unite_meshset(dmmoab->fileset, part_set);MBERRNM(merr);
690 }
691 merr = mbImpl->add_entities(dmmoab->fileset, &part_set, 1);MBERRNM(merr);
692 }
693 }
694 }
695
696 merr = mbImpl->add_parent_child(dmmoab->fileset, regionset);MBERRNM(merr);
697
698 /* Only in parallel: resolve shared entities between processors and exchange ghost layers */
699 if (global_size > 1) {
700
701 ierr = PetscLogEventBegin(genCtx.parResolve, 0, 0, 0, 0);CHKERRQ(ierr);
702
703 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, genCtx.dim, cells);MBERR("Can't get all d-dimensional elements.", merr);
704 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 0, verts);MBERR("Can't get all vertices.", merr);
705
706 if (genCtx.A * genCtx.B * genCtx.C != 1) { // merge needed
707 moab::MergeMesh mm(mbImpl);
708 if (genCtx.newMergeMethod) {
709 merr = mm.merge_using_integer_tag(verts, global_id_tag);MBERR("Can't merge with GLOBAL_ID tag", merr);
710 }
711 else {
712 merr = mm.merge_entities(cells, 0.0001);MBERR("Can't merge with coordinates", merr);
713 }
714 }
715
716 #ifdef MOAB_HAVE_MPI
717 /* check the handles */
718 merr = pcomm->check_all_shared_handles();MBERRV(mbImpl, merr);
719
720 /* resolve the shared entities by exchanging information to adjacent processors */
721 merr = pcomm->resolve_shared_ents(dmmoab->fileset, cells, dim, dim - 1, NULL, &global_id_tag);MBERRV(mbImpl, merr);
722 if (dmmoab->fileset) {
723 merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false, &dmmoab->fileset);MBERRV(mbImpl, merr);
724 }
725 else {
726 merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false);MBERRV(mbImpl, merr);
727 }
728
729 /* Reassign global IDs on all entities. */
730 merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, false, true, false);MBERRNM(merr);
731 #endif
732
733 ierr = PetscLogEventEnd(genCtx.parResolve, 0, 0, 0, 0);CHKERRQ(ierr);
734 }
735
736 if (!genCtx.keep_skins) { // default is to delete the 1- and 2-dimensional entities
737 // delete all quads and edges
738 moab::Range toDelete;
739 if (genCtx.dim > 1) {
740 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 1, toDelete);MBERR("Can't get edges", merr);
741 }
742
743 if (genCtx.dim > 2) {
744 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 2, toDelete);MBERR("Can't get faces", merr);
745 }
746
747 #ifdef MOAB_HAVE_MPI
748 merr = dmmoab->pcomm->delete_entities(toDelete) ;MBERR("Can't delete entities", merr);
749 #endif
750 }
751
752 /* set geometric dimension tag for regions */
753 merr = mbImpl->tag_set_data(geom_tag, ®ionset, 1, &dmmoab->dim);MBERRNM(merr);
754 /* set default material ID for regions */
755 int default_material = 1;
756 merr = mbImpl->tag_set_data(mat_tag, ®ionset, 1, &default_material);MBERRNM(merr);
757 /*
758 int default_dbc = 0;
759 merr = mbImpl->tag_set_data(dir_tag, &vertexset, 1, &default_dbc);MBERRNM(merr);
760 */
761 PetscFunctionReturn(0);
762 }
763
764
DMMoab_GetReadOptions_Private(PetscBool by_rank,PetscInt numproc,PetscInt dim,PetscInt nghost,MoabReadMode mode,PetscInt dbglevel,const char * dm_opts,const char * extra_opts,const char ** read_opts)765 PetscErrorCode DMMoab_GetReadOptions_Private(PetscBool by_rank, PetscInt numproc, PetscInt dim, PetscInt nghost, MoabReadMode mode, PetscInt dbglevel, const char* dm_opts, const char* extra_opts, const char** read_opts)
766 {
767 char *ropts;
768 char ropts_par[PETSC_MAX_PATH_LEN], ropts_pargh[PETSC_MAX_PATH_LEN];
769 char ropts_dbg[PETSC_MAX_PATH_LEN];
770 PetscErrorCode ierr;
771
772 PetscFunctionBegin;
773 ierr = PetscMalloc1(PETSC_MAX_PATH_LEN, &ropts);CHKERRQ(ierr);
774 ierr = PetscMemzero(&ropts_par, PETSC_MAX_PATH_LEN);CHKERRQ(ierr);
775 ierr = PetscMemzero(&ropts_pargh, PETSC_MAX_PATH_LEN);CHKERRQ(ierr);
776 ierr = PetscMemzero(&ropts_dbg, PETSC_MAX_PATH_LEN);CHKERRQ(ierr);
777
778 /* do parallel read unless using only one processor */
779 if (numproc > 1) {
780 // ierr = PetscSNPrintf(ropts_par, PETSC_MAX_PATH_LEN, "PARALLEL=%s;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=%d.0.1%s;",MoabReadModes[mode],dim,(by_rank ? ";PARTITION_BY_RANK":""));CHKERRQ(ierr);
781 ierr = PetscSNPrintf(ropts_par, PETSC_MAX_PATH_LEN, "PARALLEL=%s;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS;%s", MoabReadModes[mode], (by_rank ? "PARTITION_BY_RANK;" : ""));CHKERRQ(ierr);
782 if (nghost) {
783 ierr = PetscSNPrintf(ropts_pargh, PETSC_MAX_PATH_LEN, "PARALLEL_GHOSTS=%d.0.%d;", dim, nghost);CHKERRQ(ierr);
784 }
785 }
786
787 if (dbglevel) {
788 if (numproc > 1) {
789 ierr = PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%d;DEBUG_PIO=%d;", dbglevel, dbglevel);CHKERRQ(ierr);
790 }
791 else {
792 ierr = PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%d;", dbglevel);CHKERRQ(ierr);
793 }
794 }
795
796 ierr = PetscSNPrintf(ropts, PETSC_MAX_PATH_LEN, "%s%s%s%s%s", ropts_par, (nghost ? ropts_pargh : ""), ropts_dbg, (extra_opts ? extra_opts : ""), (dm_opts ? dm_opts : ""));CHKERRQ(ierr);
797 *read_opts = ropts;
798 PetscFunctionReturn(0);
799 }
800
801
802 /*@C
803 DMMoabLoadFromFile - Creates a DM object by loading the mesh from a user specified file.
804
805 Collective
806
807 Input Parameters:
808 + comm - The communicator for the DM object
809 . dim - The spatial dimension
810 . filename - The name of the mesh file to be loaded
811 - usrreadopts - The options string to read a MOAB mesh.
812
813 Reference (Parallel Mesh Initialization: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html#fivetwo)
814
815 Output Parameter:
816 . dm - The DM object
817
818 Level: beginner
819
820 .seealso: DMSetType(), DMCreate(), DMMoabCreateBoxMesh()
821 @*/
DMMoabLoadFromFile(MPI_Comm comm,PetscInt dim,PetscInt nghost,const char * filename,const char * usrreadopts,DM * dm)822 PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm, PetscInt dim, PetscInt nghost, const char* filename, const char* usrreadopts, DM *dm)
823 {
824 moab::ErrorCode merr;
825 PetscInt nprocs;
826 DM_Moab *dmmoab;
827 moab::Interface *mbiface;
828 #ifdef MOAB_HAVE_MPI
829 moab::ParallelComm *pcomm;
830 #endif
831 moab::Range verts, elems;
832 const char *readopts;
833 PetscErrorCode ierr;
834
835 PetscFunctionBegin;
836 PetscValidPointer(dm, 6);
837
838 /* Create the basic DMMoab object and keep the default parameters created by DM impls */
839 ierr = DMMoabCreateMoab(comm, NULL, NULL, NULL, dm);CHKERRQ(ierr);
840
841 /* get all the necessary handles from the private DM object */
842 dmmoab = (DM_Moab*)(*dm)->data;
843 mbiface = dmmoab->mbiface;
844 #ifdef MOAB_HAVE_MPI
845 pcomm = dmmoab->pcomm;
846 nprocs = pcomm->size();
847 #else
848 nprocs = 1;
849 #endif
850 /* TODO: Decipher dimension based on the loaded mesh instead of getting from user */
851 dmmoab->dim = dim;
852 dmmoab->nghostrings = nghost;
853 dmmoab->refct = 1;
854
855 /* create a file set to associate all entities in current mesh */
856 merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, dmmoab->fileset);MBERR("Creating file set failed", merr);
857
858 /* add mesh loading options specific to the DM */
859 ierr = DMMoab_GetReadOptions_Private(dmmoab->partition_by_rank, nprocs, dim, nghost, dmmoab->read_mode,
860 dmmoab->rw_dbglevel, dmmoab->extra_read_options, usrreadopts, &readopts);CHKERRQ(ierr);
861
862 PetscInfo2(*dm, "Reading file %s with options: %s\n", filename, readopts);
863
864 /* Load the mesh from a file. */
865 if (dmmoab->fileset) {
866 merr = mbiface->load_file(filename, &dmmoab->fileset, readopts);MBERRVM(mbiface, "Reading MOAB file failed.", merr);
867 }
868 else {
869 merr = mbiface->load_file(filename, 0, readopts);MBERRVM(mbiface, "Reading MOAB file failed.", merr);
870 }
871
872 #ifdef MOAB_HAVE_MPI
873 /* Reassign global IDs on all entities. */
874 /* merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, true, true, true);MBERRNM(merr); */
875 #endif
876
877 /* load the local vertices */
878 merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);MBERRNM(merr);
879 /* load the local elements */
880 merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);MBERRNM(merr);
881
882 #ifdef MOAB_HAVE_MPI
883 /* Everything is set up, now just do a tag exchange to update tags
884 on all of the ghost vertexes */
885 merr = pcomm->exchange_tags(dmmoab->ltog_tag, verts);MBERRV(mbiface, merr);
886 merr = pcomm->exchange_tags(dmmoab->ltog_tag, elems);MBERRV(mbiface, merr);
887 merr = pcomm->collective_sync_partition();MBERR("Collective sync failed", merr);
888 #endif
889
890 PetscInfo3(*dm, "MOAB file '%s' was successfully loaded. Found %D vertices and %D elements.\n", filename, verts.size(), elems.size());
891 ierr = PetscFree(readopts);CHKERRQ(ierr);
892 PetscFunctionReturn(0);
893 }
894
895
896 /*@C
897 DMMoabRenumberMeshEntities - Order and number all entities (vertices->elements) to be contiguously ordered
898 in parallel
899
900 Collective
901
902 Input Parameters:
903 . dm - The DM object
904
905 Level: advanced
906
907 .seealso: DMSetUp(), DMCreate()
908 @*/
DMMoabRenumberMeshEntities(DM dm)909 PetscErrorCode DMMoabRenumberMeshEntities(DM dm)
910 {
911 moab::Range verts;
912
913 PetscFunctionBegin;
914 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
915
916 #ifdef MOAB_HAVE_MPI
917 /* Insert new points */
918 moab::ErrorCode merr;
919 merr = ((DM_Moab*) dm->data)->pcomm->assign_global_ids(((DM_Moab*) dm->data)->fileset, 3, 0, false, true, false);MBERRNM(merr);
920 #endif
921 PetscFunctionReturn(0);
922 }
923
924