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