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