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: 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: 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: 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: 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;
204:   moab::ErrorCode          merr;

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:   PetscMalloc1(nnodes, &gids);

213:   merr = iface->get_node_coords(3, nnodes, 0, startv, arrays);
214:   MBERR("Can't get node coords.", merr);

216:   /* will start with the lower corner: */
217:   /* x = ( m * genCtx.A + a) * genCtx.q * genCtx.blockSizeElementXYZ[0]; */
218:   /* y = ( n * genCtx.B + b) * genCtx.q * genCtx.blockSizeElementXYZ[1]; */
219:   /* z = ( k * genCtx.C + c) * genCtx.q * genCtx.blockSizeElementXYZ[2]; */

221:   x = (m * genCtx.A + a) * genCtx.q;
222:   y = (n * genCtx.B + b) * genCtx.q;
223:   z = (k * genCtx.C + c) * genCtx.q;
224:   PetscInfo(NULL, "Starting offset for coordinates := %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", x, y, z);
225:   ix = 0;
226:   moab::Range verts(startv, startv + nnodes - 1);
227:   for (kk = 0; kk < (genCtx.dim > 2 ? genCtx.blockSizeVertexXYZ[2] : 1); kk++) {
228:     for (jj = 0; jj < (genCtx.dim > 1 ? genCtx.blockSizeVertexXYZ[1] : 1); jj++) {
229:       for (ii = 0; ii < genCtx.blockSizeVertexXYZ[0]; ii++, ix++) {
230:         /* set coordinates for the vertices */
231:         arrays[0][ix] = (x + ii) * genCtx.dx + genCtx.xyzbounds[0];
232:         arrays[1][ix] = (y + jj) * genCtx.dy + genCtx.xyzbounds[2];
233:         arrays[2][ix] = (z + kk) * genCtx.dz + genCtx.xyzbounds[4];
234:         PetscInfo(NULL, "Creating vertex with coordinates := %f, %f, %f\n", arrays[0][ix], arrays[1][ix], arrays[2][ix]);

236:         /* If we want to set some tags on the vertices -> use the following entity handle definition:
237:            moab::EntityHandle v = startv + ix;
238:         */
239:         /* compute the global ID for vertex */
240:         gids[ix] = 1 + (x + ii) + (y + jj) * genCtx.NX + (z + kk) * (genCtx.NX * genCtx.NY);
241:       }
242:     }
243:   }
244:   /* set global ID data on vertices */
245:   mbImpl->tag_set_data(global_id_tag, verts, &gids[0]);
246:   verts.swap(uverts);
247:   PetscFree(gids);
248:   return 0;
249: }

251: 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)
252: {
253:   moab::ErrorCode     merr;
254:   PetscInt            ix, ie, xe, ye, ze;
255:   PetscInt            ii, jj, kk, nvperelem;
256:   PetscInt            simplices_per_tensor[4] = {0, 1, 2, 6};
257:   PetscInt            ntensorelems            = genCtx.blockSizeElementXYZ[0] * (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] * (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1) : 1); /*pow(genCtx.blockSizeElement,genCtx.dim);*/
258:   PetscInt            nelems                  = ntensorelems;
259:   moab::EntityHandle  starte; /* connectivity */
260:   moab::EntityHandle *conn;

262:   switch (genCtx.dim) {
263:   case 1:
264:     nvperelem = 2;
265:     merr      = iface->get_element_connect(nelems, 2, moab::MBEDGE, 0, starte, conn);
266:     MBERR("Can't get EDGE2 element connectivity.", merr);
267:     break;
268:   case 2:
269:     if (genCtx.simplex) {
270:       nvperelem = 3;
271:       nelems    = ntensorelems * simplices_per_tensor[genCtx.dim];
272:       merr      = iface->get_element_connect(nelems, 3, moab::MBTRI, 0, starte, conn);
273:       MBERR("Can't get TRI3 element connectivity.", merr);
274:     } else {
275:       nvperelem = 4;
276:       merr      = iface->get_element_connect(nelems, 4, moab::MBQUAD, 0, starte, conn);
277:       MBERR("Can't get QUAD4 element connectivity.", merr);
278:     }
279:     break;
280:   case 3:
281:   default:
282:     if (genCtx.simplex) {
283:       nvperelem = 4;
284:       nelems    = ntensorelems * simplices_per_tensor[genCtx.dim];
285:       merr      = iface->get_element_connect(nelems, 4, moab::MBTET, 0, starte, conn);
286:       MBERR("Can't get TET4 element connectivity.", merr);
287:     } else {
288:       nvperelem = 8;
289:       merr      = iface->get_element_connect(nelems, 8, moab::MBHEX, 0, starte, conn);
290:       MBERR("Can't get HEX8 element connectivity.", merr);
291:     }
292:     break;
293:   }

295:   ix = ie = 0; /* index now in the elements, for global ids */

297:   /* create a temporary range to store local element handles */
298:   moab::Range           tmp(starte, starte + nelems - 1);
299:   std::vector<PetscInt> gids(nelems);

301:   /* identify the elements at the lower corner, for their global ids */
302:   xe = m * genCtx.A * genCtx.blockSizeElementXYZ[0] + a * genCtx.blockSizeElementXYZ[0];
303:   ye = (genCtx.dim > 1 ? n * genCtx.B * genCtx.blockSizeElementXYZ[1] + b * genCtx.blockSizeElementXYZ[1] : 0);
304:   ze = (genCtx.dim > 2 ? k * genCtx.C * genCtx.blockSizeElementXYZ[2] + c * genCtx.blockSizeElementXYZ[2] : 0);

306:   /* create owned elements requested by genCtx */
307:   for (kk = 0; kk < (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1); kk++) {
308:     for (jj = 0; jj < (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] : 1); jj++) {
309:       for (ii = 0; ii < genCtx.blockSizeElementXYZ[0]; ii++) {
310:         moab::EntityHandle corner = startv + genCtx.q * ii + genCtx.q * jj * genCtx.ystride + genCtx.q * kk * genCtx.zstride;

312:         std::pair<PetscInt, PetscInt> entoffset = DMMoab_SetElementConnectivity_Private(genCtx, ix, corner, conn);

314:         for (PetscInt j = 0; j < entoffset.second; j++) {
315:           /* The entity handle for the particular element -> if we want to set some tags is
316:              moab::EntityHandle eh = starte + ie + j;
317:           */
318:           gids[ie + j] = 1 + ((xe + ii) + (ye + jj) * genCtx.nex + (ze + kk) * (genCtx.nex * genCtx.ney));
319:           /* gids[ie+j] = ie + j + ((xe + ii) + (ye + jj) * genCtx.nex + (ze + kk) * (genCtx.nex * genCtx.ney)); */
320:           /* gids[ie+j] = 1 + ie; */
321:           /* ie++; */
322:         }

324:         ix += entoffset.first;
325:         ie += entoffset.second;
326:       }
327:     }
328:   }
329:   if (genCtx.adjEnts) { /* we need to update adjacencies now, because some elements are new */
330:     merr = iface->update_adjacencies(starte, nelems, nvperelem, conn);
331:     MBERR("Can't update adjacencies", merr);
332:   }
333:   tmp.swap(cells);
334:   merr = mbImpl->tag_set_data(global_id_tag, cells, &gids[0]);
335:   MBERR("Can't set global ids to elements.", merr);
336:   return 0;
337: }

339: PetscErrorCode DMMBUtil_InitializeOptions(DMMoabMeshGeneratorCtx &genCtx, PetscInt dim, PetscBool simplex, PetscInt rank, PetscInt nprocs, const PetscReal *bounds, PetscInt nelems)
340: {
341:   /* Initialize all genCtx data */
342:   genCtx.dim            = dim;
343:   genCtx.simplex        = simplex;
344:   genCtx.newMergeMethod = genCtx.keep_skins = genCtx.adjEnts = true;
345:   /* determine other global quantities for the mesh used for nodes increments */
346:   genCtx.q        = 1;
347:   genCtx.fraction = genCtx.remainder = genCtx.cumfraction = 0;

349:   if (!genCtx.usrxyzgrid) { /* not overridden by genCtx - assume nele equally and that genCtx wants a uniform cube mesh */

351:     genCtx.fraction    = nelems / nprocs; /* partition only by the largest dimension */
352:     genCtx.remainder   = nelems % nprocs; /* remainder after partition which gets evenly distributed by round-robin */
353:     genCtx.cumfraction = (rank > 0 ? (genCtx.fraction) * (rank) + (rank - 1 < genCtx.remainder ? rank : genCtx.remainder) : 0);
354:     if (rank < genCtx.remainder) /* This process gets "fraction+1" elements */
355:       genCtx.fraction++;

357:     PetscInfo(NULL, "Fraction = %" PetscInt_FMT ", Remainder = %" PetscInt_FMT ", Cumulative fraction = %" PetscInt_FMT "\n", genCtx.fraction, genCtx.remainder, genCtx.cumfraction);
358:     switch (genCtx.dim) {
359:     case 1:
360:       genCtx.blockSizeElementXYZ[0] = genCtx.fraction;
361:       genCtx.blockSizeElementXYZ[1] = 1;
362:       genCtx.blockSizeElementXYZ[2] = 1;
363:       break;
364:     case 2:
365:       genCtx.blockSizeElementXYZ[0] = nelems;
366:       genCtx.blockSizeElementXYZ[1] = genCtx.fraction;
367:       genCtx.blockSizeElementXYZ[2] = 1;
368:       break;
369:     case 3:
370:     default:
371:       genCtx.blockSizeElementXYZ[0] = nelems;
372:       genCtx.blockSizeElementXYZ[1] = nelems;
373:       genCtx.blockSizeElementXYZ[2] = genCtx.fraction;
374:       break;
375:     }
376:   }

378:   /* partition only by the largest dimension */
379:   /* Total number of local elements := genCtx.blockSizeElementXYZ[0]*(genCtx.dim>1? genCtx.blockSizeElementXYZ[1]*(genCtx.dim>2 ? genCtx.blockSizeElementXYZ[2]:1) :1); */
380:   if (bounds) {
381:     for (PetscInt i = 0; i < 6; i++) genCtx.xyzbounds[i] = bounds[i];
382:   } else {
383:     genCtx.xyzbounds[0] = genCtx.xyzbounds[2] = genCtx.xyzbounds[4] = 0.0;
384:     genCtx.xyzbounds[1] = genCtx.xyzbounds[3] = genCtx.xyzbounds[5] = 1.0;
385:   }

387:   if (!genCtx.usrprocgrid) {
388:     switch (genCtx.dim) {
389:     case 1:
390:       genCtx.M = nprocs;
391:       genCtx.N = genCtx.K = 1;
392:       break;
393:     case 2:
394:       genCtx.N = nprocs;
395:       genCtx.M = genCtx.K = 1;
396:       break;
397:     default:
398:       genCtx.K = nprocs;
399:       genCtx.M = genCtx.N = 1;
400:       break;
401:     }
402:   }

404:   if (!genCtx.usrrefgrid) genCtx.A = genCtx.B = genCtx.C = 1;

406:   /* more default values */
407:   genCtx.nex = genCtx.ney = genCtx.nez = 0;
408:   genCtx.xstride = genCtx.ystride = genCtx.zstride = 0;
409:   genCtx.NX = genCtx.NY = genCtx.NZ = 0;
410:   genCtx.nex = genCtx.ney = genCtx.nez = 0;
411:   genCtx.blockSizeVertexXYZ[0] = genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 1;

413:   switch (genCtx.dim) {
414:   case 3:
415:     genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1;
416:     genCtx.blockSizeVertexXYZ[1] = genCtx.q * genCtx.blockSizeElementXYZ[1] + 1;
417:     genCtx.blockSizeVertexXYZ[2] = genCtx.q * genCtx.blockSizeElementXYZ[2] + 1;

419:     genCtx.nex     = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0];               /* number of elements in x direction, used for global id on element */
420:     genCtx.dx      = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (nelems * genCtx.q); /* distance between 2 nodes in x direction */
421:     genCtx.NX      = (genCtx.q * genCtx.nex + 1);
422:     genCtx.xstride = 1;
423:     genCtx.ney     = genCtx.N * genCtx.B * genCtx.blockSizeElementXYZ[1];               /* number of elements in y direction  .... */
424:     genCtx.dy      = (genCtx.xyzbounds[3] - genCtx.xyzbounds[2]) / (nelems * genCtx.q); /* distance between 2 nodes in y direction */
425:     genCtx.NY      = (genCtx.q * genCtx.ney + 1);
426:     genCtx.ystride = genCtx.blockSizeVertexXYZ[0];
427:     genCtx.nez     = genCtx.K * genCtx.C * genCtx.blockSizeElementXYZ[2];               /* number of elements in z direction  .... */
428:     genCtx.dz      = (genCtx.xyzbounds[5] - genCtx.xyzbounds[4]) / (nelems * genCtx.q); /* distance between 2 nodes in z direction */
429:     genCtx.NZ      = (genCtx.q * genCtx.nez + 1);
430:     genCtx.zstride = genCtx.blockSizeVertexXYZ[0] * genCtx.blockSizeVertexXYZ[1];
431:     break;
432:   case 2:
433:     genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1;
434:     genCtx.blockSizeVertexXYZ[1] = genCtx.q * genCtx.blockSizeElementXYZ[1] + 1;
435:     genCtx.blockSizeVertexXYZ[2] = 0;

437:     genCtx.nex     = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0];                   /* number of elements in x direction, used for global id on element */
438:     genCtx.dx      = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (genCtx.nex * genCtx.q); /* distance between 2 nodes in x direction */
439:     genCtx.NX      = (genCtx.q * genCtx.nex + 1);
440:     genCtx.xstride = 1;
441:     genCtx.ney     = genCtx.N * genCtx.B * genCtx.blockSizeElementXYZ[1];               /* number of elements in y direction  .... */
442:     genCtx.dy      = (genCtx.xyzbounds[3] - genCtx.xyzbounds[2]) / (nelems * genCtx.q); /* distance between 2 nodes in y direction */
443:     genCtx.NY      = (genCtx.q * genCtx.ney + 1);
444:     genCtx.ystride = genCtx.blockSizeVertexXYZ[0];
445:     break;
446:   case 1:
447:     genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 0;
448:     genCtx.blockSizeVertexXYZ[0]                                = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1;

450:     genCtx.nex     = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0];               /* number of elements in x direction, used for global id on element */
451:     genCtx.dx      = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (nelems * genCtx.q); /* distance between 2 nodes in x direction */
452:     genCtx.NX      = (genCtx.q * genCtx.nex + 1);
453:     genCtx.xstride = 1;
454:     break;
455:   }

457:   /* Lets check for some valid input */
460:              genCtx.N, genCtx.K, nprocs);
461:   /* validate the bounds data */

466:   PetscInfo(NULL, "Local elements:= %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.blockSizeElementXYZ[0], genCtx.blockSizeElementXYZ[1], genCtx.blockSizeElementXYZ[2]);
467:   PetscInfo(NULL, "Local vertices:= %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.blockSizeVertexXYZ[0], genCtx.blockSizeVertexXYZ[1], genCtx.blockSizeVertexXYZ[2]);
468:   PetscInfo(NULL, "Local blocks/processors := %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.A, genCtx.B, genCtx.C);
469:   PetscInfo(NULL, "Local processors := %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.M, genCtx.N, genCtx.K);
470:   PetscInfo(NULL, "Local nexyz:= %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.nex, genCtx.ney, genCtx.nez);
471:   PetscInfo(NULL, "Local delxyz:= %g, %g, %g\n", genCtx.dx, genCtx.dy, genCtx.dz);
472:   PetscInfo(NULL, "Local strides:= %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.xstride, genCtx.ystride, genCtx.zstride);
473:   return 0;
474: }

476: /*@C
477:   DMMoabCreateBoxMesh - Creates a mesh on the tensor product (box) of intervals with genCtx specified bounds.

479:   Collective

481:   Input Parameters:
482: + comm - The communicator for the DM object
483: . dim - The spatial dimension
484: . 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
485: . nele - The number of discrete elements in each direction
486: - user_nghost - The number of ghosted layers needed in the partitioned mesh

488:   Output Parameter:
489: . dm  - The DM object

491:   Level: beginner

493: .seealso: `DMSetType()`, `DMCreate()`, `DMMoabLoadFromFile()`
494: @*/
495: PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool useSimplex, const PetscReal *bounds, PetscInt nele, PetscInt nghost, DM *dm)
496: {
497:   moab::ErrorCode  merr;
498:   PetscInt         a, b, c, n, global_size, global_rank;
499:   DM_Moab         *dmmoab;
500:   moab::Interface *mbImpl;
501: #ifdef MOAB_HAVE_MPI
502:   moab::ParallelComm *pcomm;
503: #endif
504:   moab::ReadUtilIface   *readMeshIface;
505:   moab::Range            verts, cells, edges, faces, adj, dim3, dim2;
506:   DMMoabMeshGeneratorCtx genCtx;
507:   const PetscInt         npts = nele + 1; /* Number of points in every dimension */

509:   moab::Tag          global_id_tag, part_tag, geom_tag, mat_tag, dir_tag, neu_tag;
510:   moab::Range        ownedvtx, ownedelms, localvtxs, localelms;
511:   moab::EntityHandle regionset;
512:   PetscInt           ml = 0, nl = 0, kl = 0;


516:   PetscLogEventRegister("GenerateMesh", DM_CLASSID, &genCtx.generateMesh);
517:   PetscLogEventRegister("AddVertices", DM_CLASSID, &genCtx.generateVertices);
518:   PetscLogEventRegister("AddElements", DM_CLASSID, &genCtx.generateElements);
519:   PetscLogEventRegister("ParResolve", DM_CLASSID, &genCtx.parResolve);
520:   PetscLogEventBegin(genCtx.generateMesh, 0, 0, 0, 0);
521:   MPI_Comm_size(comm, &global_size);
522:   /* total number of vertices in all dimensions */
523:   n = pow(npts, dim);

525:   /* do some error checking */

530:   /* Create the basic DMMoab object and keep the default parameters created by DM impls */
531:   DMMoabCreateMoab(comm, NULL, NULL, NULL, dm);

533:   /* get all the necessary handles from the private DM object */
534:   dmmoab = (DM_Moab *)(*dm)->data;
535:   mbImpl = dmmoab->mbiface;
536: #ifdef MOAB_HAVE_MPI
537:   pcomm       = dmmoab->pcomm;
538:   global_rank = pcomm->rank();
539: #else
540:   global_rank = 0;
541:   global_size = 1;
542: #endif
543:   global_id_tag       = dmmoab->ltog_tag;
544:   dmmoab->dim         = dim;
545:   dmmoab->nghostrings = nghost;
546:   dmmoab->refct       = 1;

548:   /* create a file set to associate all entities in current mesh */
549:   merr = mbImpl->create_meshset(moab::MESHSET_SET, dmmoab->fileset);
550:   MBERR("Creating file set failed", merr);

552:   /* No errors yet; proceed with building the mesh */
553:   merr = mbImpl->query_interface(readMeshIface);
554:   MBERRNM(merr);

556:   genCtx.M = genCtx.N = genCtx.K = 1;
557:   genCtx.A = genCtx.B = genCtx.C = 1;
558:   genCtx.blockSizeElementXYZ[0]  = 0;
559:   genCtx.blockSizeElementXYZ[1]  = 0;
560:   genCtx.blockSizeElementXYZ[2]  = 0;

562:   PetscOptionsBegin(comm, "", "DMMoab Creation Options", "DMMOAB");
563:   /* Handle DMMoab spatial resolution */
564:   PetscOptionsInt("-dmb_grid_x", "Number of grid points in x direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[0], &genCtx.blockSizeElementXYZ[0], &genCtx.usrxyzgrid);
565:   if (dim > 1) PetscOptionsInt("-dmb_grid_y", "Number of grid points in y direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[1], &genCtx.blockSizeElementXYZ[1], &genCtx.usrxyzgrid);
566:   if (dim > 2) PetscOptionsInt("-dmb_grid_z", "Number of grid points in z direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[2], &genCtx.blockSizeElementXYZ[2], &genCtx.usrxyzgrid);

568:   /* Handle DMMoab parallel distribution */
569:   PetscOptionsInt("-dmb_processors_x", "Number of processors in x direction", "DMMoabSetNumProcs", genCtx.M, &genCtx.M, &genCtx.usrprocgrid);
570:   if (dim > 1) PetscOptionsInt("-dmb_processors_y", "Number of processors in y direction", "DMMoabSetNumProcs", genCtx.N, &genCtx.N, &genCtx.usrprocgrid);
571:   if (dim > 2) PetscOptionsInt("-dmb_processors_z", "Number of processors in z direction", "DMMoabSetNumProcs", genCtx.K, &genCtx.K, &genCtx.usrprocgrid);

573:   /* Handle DMMoab block refinement */
574:   PetscOptionsInt("-dmb_refine_x", "Number of refinement blocks in x direction", "DMMoabSetRefinement", genCtx.A, &genCtx.A, &genCtx.usrrefgrid);
575:   if (dim > 1) PetscOptionsInt("-dmb_refine_y", "Number of refinement blocks in y direction", "DMMoabSetRefinement", genCtx.B, &genCtx.B, &genCtx.usrrefgrid);
576:   if (dim > 2) PetscOptionsInt("-dmb_refine_z", "Number of refinement blocks in z direction", "DMMoabSetRefinement", genCtx.C, &genCtx.C, &genCtx.usrrefgrid);
577:   PetscOptionsEnd();

579:   DMMBUtil_InitializeOptions(genCtx, dim, useSimplex, global_rank, global_size, bounds, nele);


583:   if (genCtx.adjEnts) genCtx.keep_skins = true; /* do not delete anything - consumes more memory */

585:   /* determine m, n, k for processor rank */
586:   ml = nl = kl = 0;
587:   switch (genCtx.dim) {
588:   case 1:
589:     ml = (genCtx.cumfraction);
590:     break;
591:   case 2:
592:     nl = (genCtx.cumfraction);
593:     break;
594:   default:
595:     kl = (genCtx.cumfraction) / genCtx.q / genCtx.blockSizeElementXYZ[2] / genCtx.C; //genCtx.K
596:     break;
597:   }

599:   /*
600:    * so there are a total of M * A * blockSizeElement elements in x direction (so M * A * blockSizeElement + 1 verts in x direction)
601:    * so there are a total of N * B * blockSizeElement elements in y direction (so N * B * blockSizeElement + 1 verts in y direction)
602:    * so there are a total of K * C * blockSizeElement elements in z direction (so K * C * blockSizeElement + 1 verts in z direction)

604:    * there are ( M * A blockSizeElement)      *  ( N * B * blockSizeElement)      * (K * C * blockSizeElement)    hexas
605:    * there are ( M * A * blockSizeElement + 1) *  ( N * B * blockSizeElement + 1) * (K * C * blockSizeElement + 1) vertices
606:    * x is the first dimension that varies
607:    */

609:   /* generate the block at (a, b, c); it will represent a partition , it will get a partition tag */
610:   PetscInt dum_id = -1;
611:   merr            = mbImpl->tag_get_handle("GLOBAL_ID", 1, moab::MB_TYPE_INTEGER, global_id_tag);
612:   MBERR("Getting Global_ID Tag handle failed", merr);

614:   merr = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, mat_tag);
615:   MBERR("Getting Material set Tag handle failed", merr);
616:   merr = mbImpl->tag_get_handle(DIRICHLET_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, dir_tag);
617:   MBERR("Getting Dirichlet set Tag handle failed", merr);
618:   merr = mbImpl->tag_get_handle(NEUMANN_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, neu_tag);
619:   MBERR("Getting Neumann set Tag handle failed", merr);

621:   merr = mbImpl->tag_get_handle("PARALLEL_PARTITION", 1, moab::MB_TYPE_INTEGER, part_tag, moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &dum_id);
622:   MBERR("Getting Partition Tag handle failed", merr);

624:   /* lets create some sets */
625:   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);
626:   MBERRNM(merr);
627:   merr = mbImpl->create_meshset(moab::MESHSET_SET, regionset);
628:   MBERRNM(merr);
629:   PetscLogEventEnd(genCtx.generateMesh, 0, 0, 0, 0);

631:   for (a = 0; a < (genCtx.dim > 0 ? genCtx.A : genCtx.A); a++) {
632:     for (b = 0; b < (genCtx.dim > 1 ? genCtx.B : 1); b++) {
633:       for (c = 0; c < (genCtx.dim > 2 ? genCtx.C : 1); c++) {
634:         moab::EntityHandle startv;

636:         PetscLogEventBegin(genCtx.generateVertices, 0, 0, 0, 0);
637:         DMMoab_GenerateVertices_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, verts);
638:         PetscLogEventEnd(genCtx.generateVertices, 0, 0, 0, 0);

640:         PetscLogEventBegin(genCtx.generateElements, 0, 0, 0, 0);
641:         DMMoab_GenerateElements_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, cells);
642:         PetscLogEventEnd(genCtx.generateElements, 0, 0, 0, 0);

644:         PetscInt part_num = 0;
645:         switch (genCtx.dim) {
646:         case 3:
647:           part_num += (c + kl * genCtx.C) * (genCtx.M * genCtx.A * genCtx.N * genCtx.B);
648:         case 2:
649:           part_num += (b + nl * genCtx.B) * (genCtx.M * genCtx.A);
650:         case 1:
651:           part_num += (a + ml * genCtx.A);
652:           break;
653:         }

655:         moab::EntityHandle part_set;
656:         merr = mbImpl->create_meshset(moab::MESHSET_SET, part_set);
657:         MBERR("Can't create mesh set.", merr);

659:         merr = mbImpl->add_entities(part_set, verts);
660:         MBERR("Can't add vertices to set.", merr);
661:         merr = mbImpl->add_entities(part_set, cells);
662:         MBERR("Can't add entities to set.", merr);
663:         merr = mbImpl->add_entities(regionset, cells);
664:         MBERR("Can't add entities to set.", merr);

666:         /* if needed, add all edges and faces */
667:         if (genCtx.adjEnts) {
668:           if (genCtx.dim > 1) {
669:             merr = mbImpl->get_adjacencies(cells, 1, true, edges, moab::Interface::UNION);
670:             MBERR("Can't get edges", merr);
671:             merr = mbImpl->add_entities(part_set, edges);
672:             MBERR("Can't add edges to partition set.", merr);
673:           }
674:           if (genCtx.dim > 2) {
675:             merr = mbImpl->get_adjacencies(cells, 2, true, faces, moab::Interface::UNION);
676:             MBERR("Can't get faces", merr);
677:             merr = mbImpl->add_entities(part_set, faces);
678:             MBERR("Can't add faces to partition set.", merr);
679:           }
680:           edges.clear();
681:           faces.clear();
682:         }
683:         verts.clear();
684:         cells.clear();

686:         merr = mbImpl->tag_set_data(part_tag, &part_set, 1, &part_num);
687:         MBERR("Can't set part tag on set", merr);
688:         if (dmmoab->fileset) {
689:           merr = mbImpl->add_parent_child(dmmoab->fileset, part_set);
690:           MBERR("Can't add part set to file set.", merr);
691:           merr = mbImpl->unite_meshset(dmmoab->fileset, part_set);
692:           MBERRNM(merr);
693:         }
694:         merr = mbImpl->add_entities(dmmoab->fileset, &part_set, 1);
695:         MBERRNM(merr);
696:       }
697:     }
698:   }

700:   merr = mbImpl->add_parent_child(dmmoab->fileset, regionset);
701:   MBERRNM(merr);

703:   /* Only in parallel: resolve shared entities between processors and exchange ghost layers */
704:   if (global_size > 1) {
705:     PetscLogEventBegin(genCtx.parResolve, 0, 0, 0, 0);

707:     merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, genCtx.dim, cells);
708:     MBERR("Can't get all d-dimensional elements.", merr);
709:     merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 0, verts);
710:     MBERR("Can't get all vertices.", merr);

712:     if (genCtx.A * genCtx.B * genCtx.C != 1) { //  merge needed
713:       moab::MergeMesh mm(mbImpl);
714:       if (genCtx.newMergeMethod) {
715:         merr = mm.merge_using_integer_tag(verts, global_id_tag);
716:         MBERR("Can't merge with GLOBAL_ID tag", merr);
717:       } else {
718:         merr = mm.merge_entities(cells, 0.0001);
719:         MBERR("Can't merge with coordinates", merr);
720:       }
721:     }

723: #ifdef MOAB_HAVE_MPI
724:     /* check the handles */
725:     merr = pcomm->check_all_shared_handles();
726:     MBERRV(mbImpl, merr);

728:     /* resolve the shared entities by exchanging information to adjacent processors */
729:     merr = pcomm->resolve_shared_ents(dmmoab->fileset, cells, dim, dim - 1, NULL, &global_id_tag);
730:     MBERRV(mbImpl, merr);
731:     if (dmmoab->fileset) {
732:       merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false, &dmmoab->fileset);
733:       MBERRV(mbImpl, merr);
734:     } else {
735:       merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false);
736:       MBERRV(mbImpl, merr);
737:     }

739:     /* Reassign global IDs on all entities. */
740:     merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, false, true, false);
741:     MBERRNM(merr);
742: #endif

744:     PetscLogEventEnd(genCtx.parResolve, 0, 0, 0, 0);
745:   }

747:   if (!genCtx.keep_skins) { // default is to delete the 1- and 2-dimensional entities
748:     // delete all quads and edges
749:     moab::Range toDelete;
750:     if (genCtx.dim > 1) {
751:       merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 1, toDelete);
752:       MBERR("Can't get edges", merr);
753:     }

755:     if (genCtx.dim > 2) {
756:       merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 2, toDelete);
757:       MBERR("Can't get faces", merr);
758:     }

760: #ifdef MOAB_HAVE_MPI
761:     merr = dmmoab->pcomm->delete_entities(toDelete);
762:     MBERR("Can't delete entities", merr);
763: #endif
764:   }

766:   /* set geometric dimension tag for regions */
767:   merr = mbImpl->tag_set_data(geom_tag, &regionset, 1, &dmmoab->dim);
768:   MBERRNM(merr);
769:   /* set default material ID for regions */
770:   int default_material = 1;
771:   merr                 = mbImpl->tag_set_data(mat_tag, &regionset, 1, &default_material);
772:   MBERRNM(merr);
773:   /*
774:     int default_dbc = 0;
775:     merr = mbImpl->tag_set_data(dir_tag, &vertexset, 1, &default_dbc);MBERRNM(merr);
776:   */
777:   return 0;
778: }

780: 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)
781: {
782:   char *ropts;
783:   char  ropts_par[PETSC_MAX_PATH_LEN], ropts_pargh[PETSC_MAX_PATH_LEN];
784:   char  ropts_dbg[PETSC_MAX_PATH_LEN];

786:   PetscMalloc1(PETSC_MAX_PATH_LEN, &ropts);
787:   PetscMemzero(&ropts_par, PETSC_MAX_PATH_LEN);
788:   PetscMemzero(&ropts_pargh, PETSC_MAX_PATH_LEN);
789:   PetscMemzero(&ropts_dbg, PETSC_MAX_PATH_LEN);

791:   /* do parallel read unless using only one processor */
792:   if (numproc > 1) {
793:     // 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":""));
794:     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;" : ""));
795:     if (nghost) PetscSNPrintf(ropts_pargh, PETSC_MAX_PATH_LEN, "PARALLEL_GHOSTS=%" PetscInt_FMT ".0.%" PetscInt_FMT ";", dim, nghost);
796:   }

798:   if (dbglevel) {
799:     if (numproc > 1) {
800:       PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%" PetscInt_FMT ";DEBUG_PIO=%" PetscInt_FMT ";", dbglevel, dbglevel);
801:     } else PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%" PetscInt_FMT ";", dbglevel);
802:   }

804:   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 : ""));
805:   *read_opts = ropts;
806:   return 0;
807: }

809: /*@C
810:   DMMoabLoadFromFile - Creates a DM object by loading the mesh from a user specified file.

812:   Collective

814:   Input Parameters:
815: + comm - The communicator for the DM object
816: . dim - The spatial dimension
817: . filename - The name of the mesh file to be loaded
818: - usrreadopts - The options string to read a MOAB mesh.

820:   Reference (Parallel Mesh Initialization: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html#fivetwo)

822:   Output Parameter:
823: . dm  - The DM object

825:   Level: beginner

827: .seealso: `DMSetType()`, `DMCreate()`, `DMMoabCreateBoxMesh()`
828: @*/
829: PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm, PetscInt dim, PetscInt nghost, const char *filename, const char *usrreadopts, DM *dm)
830: {
831:   moab::ErrorCode  merr;
832:   PetscInt         nprocs;
833:   DM_Moab         *dmmoab;
834:   moab::Interface *mbiface;
835: #ifdef MOAB_HAVE_MPI
836:   moab::ParallelComm *pcomm;
837: #endif
838:   moab::Range verts, elems;
839:   const char *readopts;


843:   /* Create the basic DMMoab object and keep the default parameters created by DM impls */
844:   DMMoabCreateMoab(comm, NULL, NULL, NULL, dm);

846:   /* get all the necessary handles from the private DM object */
847:   dmmoab  = (DM_Moab *)(*dm)->data;
848:   mbiface = dmmoab->mbiface;
849: #ifdef MOAB_HAVE_MPI
850:   pcomm  = dmmoab->pcomm;
851:   nprocs = pcomm->size();
852: #else
853:   nprocs      = 1;
854: #endif
855:   /* TODO: Decipher dimension based on the loaded mesh instead of getting from user */
856:   dmmoab->dim         = dim;
857:   dmmoab->nghostrings = nghost;
858:   dmmoab->refct       = 1;

860:   /* create a file set to associate all entities in current mesh */
861:   merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, dmmoab->fileset);
862:   MBERR("Creating file set failed", merr);

864:   /* add mesh loading options specific to the DM */
865:   DMMoab_GetReadOptions_Private(dmmoab->partition_by_rank, nprocs, dim, nghost, dmmoab->read_mode, dmmoab->rw_dbglevel, dmmoab->extra_read_options, usrreadopts, &readopts);

867:   PetscInfo(*dm, "Reading file %s with options: %s\n", filename, readopts);

869:   /* Load the mesh from a file. */
870:   if (dmmoab->fileset) {
871:     merr = mbiface->load_file(filename, &dmmoab->fileset, readopts);
872:     MBERRVM(mbiface, "Reading MOAB file failed.", merr);
873:   } else {
874:     merr = mbiface->load_file(filename, 0, readopts);
875:     MBERRVM(mbiface, "Reading MOAB file failed.", merr);
876:   }

878: #ifdef MOAB_HAVE_MPI
879:   /* Reassign global IDs on all entities. */
880:   /* merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, true, true, true);MBERRNM(merr); */
881: #endif

883:   /* load the local vertices */
884:   merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);
885:   MBERRNM(merr);
886:   /* load the local elements */
887:   merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);
888:   MBERRNM(merr);

890: #ifdef MOAB_HAVE_MPI
891:   /* Everything is set up, now just do a tag exchange to update tags
892:      on all of the ghost vertexes */
893:   merr = pcomm->exchange_tags(dmmoab->ltog_tag, verts);
894:   MBERRV(mbiface, merr);
895:   merr = pcomm->exchange_tags(dmmoab->ltog_tag, elems);
896:   MBERRV(mbiface, merr);
897:   merr = pcomm->collective_sync_partition();
898:   MBERR("Collective sync failed", merr);
899: #endif

901:   PetscInfo(*dm, "MOAB file '%s' was successfully loaded. Found %zu vertices and %zu elements.\n", filename, verts.size(), elems.size());
902:   PetscFree(readopts);
903:   return 0;
904: }

906: /*@C
907:   DMMoabRenumberMeshEntities - Order and number all entities (vertices->elements) to be contiguously ordered
908:   in parallel

910:   Collective

912:   Input Parameters:
913: . dm  - The DM object

915:   Level: advanced

917: .seealso: `DMSetUp()`, `DMCreate()`
918: @*/
919: PetscErrorCode DMMoabRenumberMeshEntities(DM dm)
920: {
921:   moab::Range verts;


925: #ifdef MOAB_HAVE_MPI
926:   /* Insert new points */
927:   moab::ErrorCode merr;
928:   merr = ((DM_Moab *)dm->data)->pcomm->assign_global_ids(((DM_Moab *)dm->data)->fileset, 3, 0, false, true, false);
929:   MBERRNM(merr);
930: #endif
931:   return 0;
932: }