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:       }
 86:       else {        /* 2 3 0 of a QUAD */
 87:         connectivity[offset + subent_conn[0]] = D;
 88:         connectivity[offset + subent_conn[1]] = A;
 89:         connectivity[offset + subent_conn[2]] = C;
 90:       }
 91:     }
 92:     else if (trigen_opts == 2) {
 93:       if (subelem) { /* 0 1 2 of a QUAD */
 94:         connectivity[offset + subent_conn[0]] = A;
 95:         connectivity[offset + subent_conn[1]] = B;
 96:         connectivity[offset + subent_conn[2]] = D;
 97:       }
 98:       else {        /* 2 3 0 of a QUAD */
 99:         connectivity[offset + subent_conn[0]] = C;
100:         connectivity[offset + subent_conn[1]] = D;
101:         connectivity[offset + subent_conn[2]] = B;
102:       }
103:     }
104:     else {
105:       switch (subelem) { /* 0 1 2 of a QUAD */
106:       case 0:
107:         connectivity[offset + subent_conn[0]] = A;
108:         connectivity[offset + subent_conn[1]] = B;
109:         connectivity[offset + subent_conn[2]] = M;
110:         break;
111:       case 1:
112:         connectivity[offset + subent_conn[0]] = B;
113:         connectivity[offset + subent_conn[1]] = C;
114:         connectivity[offset + subent_conn[2]] = M;
115:         break;
116:       case 2:
117:         connectivity[offset + subent_conn[0]] = C;
118:         connectivity[offset + subent_conn[1]] = D;
119:         connectivity[offset + subent_conn[2]] = M;
120:         break;
121:       case 3:
122:         connectivity[offset + subent_conn[0]] = D;
123:         connectivity[offset + subent_conn[1]] = A;
124:         connectivity[offset + subent_conn[2]] = M;
125:         break;
126:       }
127:     }
128:     break;
129:   case 3:
130:   default:
131:     C = corner + 1 + genCtx.ystride;
132:     D = corner +     genCtx.ystride;
133:     E = corner +                      genCtx.zstride;
134:     F = corner + 1 +                  genCtx.zstride;
135:     G = corner + 1 + genCtx.ystride + genCtx.zstride;
136:     H = corner +     genCtx.ystride + genCtx.zstride;
137:     subent_conn.resize(4);  /* only linear TET supported */
138:     moab::CN::SubEntityVertexIndices(moab::MBTET, 3, 0, subent_conn.data());
139:     switch (subelem) {
140:     case 0: /* 4 3 7 6 of a HEX */
141:       connectivity[offset + subent_conn[0]] = E;
142:       connectivity[offset + subent_conn[1]] = D;
143:       connectivity[offset + subent_conn[2]] = H;
144:       connectivity[offset + subent_conn[3]] = G;
145:       break;
146:     case 1: /* 0 1 2 5 of a HEX */
147:       connectivity[offset + subent_conn[0]] = A;
148:       connectivity[offset + subent_conn[1]] = B;
149:       connectivity[offset + subent_conn[2]] = C;
150:       connectivity[offset + subent_conn[3]] = F;
151:       break;
152:     case 2: /* 0 3 4 5 of a HEX */
153:       connectivity[offset + subent_conn[0]] = A;
154:       connectivity[offset + subent_conn[1]] = D;
155:       connectivity[offset + subent_conn[2]] = E;
156:       connectivity[offset + subent_conn[3]] = F;
157:       break;
158:     case 3: /* 2 6 3 5 of a HEX */
159:       connectivity[offset + subent_conn[0]] = C;
160:       connectivity[offset + subent_conn[1]] = G;
161:       connectivity[offset + subent_conn[2]] = D;
162:       connectivity[offset + subent_conn[3]] = F;
163:       break;
164:     case 4: /* 0 2 3 5 of a HEX */
165:       connectivity[offset + subent_conn[0]] = A;
166:       connectivity[offset + subent_conn[1]] = C;
167:       connectivity[offset + subent_conn[2]] = D;
168:       connectivity[offset + subent_conn[3]] = F;
169:       break;
170:     case 5: /* 3 6 4 5 of a HEX */
171:       connectivity[offset + subent_conn[0]] = D;
172:       connectivity[offset + subent_conn[1]] = G;
173:       connectivity[offset + subent_conn[2]] = E;
174:       connectivity[offset + subent_conn[3]] = F;
175:       break;
176:     }
177:     break;
178:   }
179:   return subent_conn.size();
180: }

182: std::pair<PetscInt, PetscInt> DMMoab_SetElementConnectivity_Private(DMMoabMeshGeneratorCtx& genCtx, PetscInt offset, PetscInt corner, moab::EntityHandle *connectivity)
183: {
184:   PetscInt vcount = 0;
185:   PetscInt simplices_per_tensor[4] = {0, 1, 2, 6};
186:   std::vector<PetscInt> subent_conn;  /* only linear edge, tri, tet supported now */
187:   subent_conn.reserve(27);
188:   PetscInt m, subelem;
189:   if (genCtx.simplex) {
190:     subelem = simplices_per_tensor[genCtx.dim];
191:     for (m = 0; m < subelem; m++) {
192:       vcount = DMMoab_SetSimplexElementConnectivity_Private(genCtx, m, offset, corner, subent_conn, connectivity);
193:       offset += vcount;
194:     }
195:   }
196:   else {
197:     subelem = 1;
198:     vcount = DMMoab_SetTensorElementConnectivity_Private(genCtx, offset, corner, subent_conn, connectivity);
199:   }
200:   return std::pair<PetscInt, PetscInt>(vcount * subelem, subelem);
201: }

203: PetscErrorCode DMMoab_GenerateVertices_Private(moab::Interface *mbImpl, moab::ReadUtilIface *iface, DMMoabMeshGeneratorCtx& genCtx, PetscInt m, PetscInt n, PetscInt k,
204:     PetscInt a, PetscInt b, PetscInt c, moab::Tag& global_id_tag, moab::EntityHandle& startv, moab::Range& uverts)
205: {
206:   PetscInt x, y, z, ix, nnodes;
207:   PetscInt ii, jj, kk;
208:   std::vector<PetscReal*> arrays;
209:   PetscInt* gids;
210:   moab::ErrorCode merr;

212:   /* we will generate (q*block+1)^3 vertices, and block^3 hexas; q is 1 for linear, 2 for quadratic
213:    * the global id of the vertices will come from m, n, k, a, b, c
214:    * x will vary from  m*A*q*block + a*q*block to m*A*q*block+(a+1)*q*block etc.
215:    */
216:   nnodes = genCtx.blockSizeVertexXYZ[0] * (genCtx.dim > 1 ? genCtx.blockSizeVertexXYZ[1] * (genCtx.dim > 2 ? genCtx.blockSizeVertexXYZ[2] : 1) : 1);
217:   PetscMalloc1(nnodes, &gids);

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

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

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

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

256: PetscErrorCode DMMoab_GenerateElements_Private(moab::Interface* mbImpl, moab::ReadUtilIface* iface, DMMoabMeshGeneratorCtx& genCtx, PetscInt m, PetscInt n, PetscInt k,
257:     PetscInt a, PetscInt b, PetscInt c, moab::Tag& global_id_tag, moab::EntityHandle startv, moab::Range& cells)
258: {
259:   moab::ErrorCode merr;
260:   PetscInt ix, ie, xe, ye, ze;
261:   PetscInt ii, jj, kk, nvperelem;
262:   PetscInt simplices_per_tensor[4] = {0, 1, 2, 6};
263:   PetscInt ntensorelems = genCtx.blockSizeElementXYZ[0] * (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] * (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1) : 1); /*pow(genCtx.blockSizeElement,genCtx.dim);*/
264:   PetscInt nelems = ntensorelems;
265:   moab::EntityHandle starte; /* connectivity */
266:   moab::EntityHandle* conn;

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

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

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

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

309:   /* create owned elements requested by genCtx */
310:   for (kk = 0; kk < (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1); kk++) {
311:     for (jj = 0; jj < (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] : 1); jj++) {
312:       for (ii = 0; ii < genCtx.blockSizeElementXYZ[0]; ii++) {

314:         moab::EntityHandle corner = startv + genCtx.q * ii + genCtx.q * jj * genCtx.ystride + genCtx.q * kk * genCtx.zstride;

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

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

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

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

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

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

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

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

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

407:   if (!genCtx.usrrefgrid) {
408:     genCtx.A = genCtx.B = genCtx.C = 1;
409:   }

411:   /* more default values */
412:   genCtx.nex = genCtx.ney = genCtx.nez = 0;
413:   genCtx.xstride = genCtx.ystride = genCtx.zstride = 0;
414:   genCtx.NX = genCtx.NY = genCtx.NZ = 0;
415:   genCtx.nex = genCtx.ney = genCtx.nez = 0;
416:   genCtx.blockSizeVertexXYZ[0] = genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 1;

418:   switch (genCtx.dim) {
419:   case 3:
420:     genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1;
421:     genCtx.blockSizeVertexXYZ[1] = genCtx.q * genCtx.blockSizeElementXYZ[1] + 1;
422:     genCtx.blockSizeVertexXYZ[2] = genCtx.q * genCtx.blockSizeElementXYZ[2] + 1;

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

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

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

462:   /* Lets check for some valid input */
465:   /* validate the bounds data */

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

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

483:   Collective

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

492:   Output Parameter:
493: . dm  - The DM object

495:   Level: beginner

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

514:   moab::Tag              global_id_tag, part_tag, geom_tag, mat_tag, dir_tag, neu_tag;
515:   moab::Range            ownedvtx, ownedelms, localvtxs, localelms;
516:   moab::EntityHandle     regionset;
517:   PetscInt               ml = 0, nl = 0, kl = 0;


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

530:   /* do some error checking */

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

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

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

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

559:   genCtx.M = genCtx.N = genCtx.K = 1;
560:   genCtx.A = genCtx.B = genCtx.C = 1;
561:   genCtx.blockSizeElementXYZ[0] = 0;
562:   genCtx.blockSizeElementXYZ[1] = 0;
563:   genCtx.blockSizeElementXYZ[2] = 0;

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

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

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

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


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

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

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

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

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

616:   merr = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, mat_tag);MBERR("Getting Material set Tag handle failed", merr);
617:   merr = mbImpl->tag_get_handle(DIRICHLET_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, dir_tag);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);MBERR("Getting Neumann set Tag handle failed", merr);

620:   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);

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

627:   for (a = 0; a < (genCtx.dim > 0 ? genCtx.A : genCtx.A); a++) {
628:     for (b = 0; b < (genCtx.dim > 1 ? genCtx.B : 1); b++) {
629:       for (c = 0; c < (genCtx.dim > 2 ? genCtx.C : 1); c++) {

631:         moab::EntityHandle startv;

633:         PetscLogEventBegin(genCtx.generateVertices, 0, 0, 0, 0);
634:         DMMoab_GenerateVertices_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, verts);
635:         PetscLogEventEnd(genCtx.generateVertices, 0, 0, 0, 0);

637:         PetscLogEventBegin(genCtx.generateElements, 0, 0, 0, 0);
638:         DMMoab_GenerateElements_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, cells);
639:         PetscLogEventEnd(genCtx.generateElements, 0, 0, 0, 0);

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

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

655:         merr = mbImpl->add_entities(part_set, verts);MBERR("Can't add vertices to set.", merr);
656:         merr = mbImpl->add_entities(part_set, cells);MBERR("Can't add entities to set.", merr);
657:         merr = mbImpl->add_entities(regionset, cells);MBERR("Can't add entities to set.", merr);

659:         /* if needed, add all edges and faces */
660:         if (genCtx.adjEnts)
661:         {
662:           if (genCtx.dim > 1) {
663:             merr = mbImpl->get_adjacencies(cells, 1, true, edges, moab::Interface::UNION);MBERR("Can't get edges", merr);
664:             merr = mbImpl->add_entities(part_set, edges);MBERR("Can't add edges to partition set.", merr);
665:           }
666:           if (genCtx.dim > 2) {
667:             merr = mbImpl->get_adjacencies(cells, 2, true, faces, moab::Interface::UNION);MBERR("Can't get faces", merr);
668:             merr = mbImpl->add_entities(part_set, faces);MBERR("Can't add faces to partition set.", merr);
669:           }
670:           edges.clear();
671:           faces.clear();
672:         }
673:         verts.clear(); cells.clear();

675:         merr = mbImpl->tag_set_data(part_tag, &part_set, 1, &part_num);MBERR("Can't set part tag on set", merr);
676:         if (dmmoab->fileset) {
677:           merr = mbImpl->add_parent_child(dmmoab->fileset, part_set);MBERR("Can't add part set to file set.", merr);
678:           merr = mbImpl->unite_meshset(dmmoab->fileset, part_set);MBERRNM(merr);
679:         }
680:         merr = mbImpl->add_entities(dmmoab->fileset, &part_set, 1);MBERRNM(merr);
681:       }
682:     }
683:   }

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

687:   /* Only in parallel: resolve shared entities between processors and exchange ghost layers */
688:   if (global_size > 1) {

690:     PetscLogEventBegin(genCtx.parResolve, 0, 0, 0, 0);

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

695:     if (genCtx.A * genCtx.B * genCtx.C != 1) { //  merge needed
696:       moab::MergeMesh mm(mbImpl);
697:       if (genCtx.newMergeMethod) {
698:         merr = mm.merge_using_integer_tag(verts, global_id_tag);MBERR("Can't merge with GLOBAL_ID tag", merr);
699:       }
700:       else {
701:         merr = mm.merge_entities(cells, 0.0001);MBERR("Can't merge with coordinates", merr);
702:       }
703:     }

705: #ifdef MOAB_HAVE_MPI
706:     /* check the handles */
707:     merr = pcomm->check_all_shared_handles();MBERRV(mbImpl, merr);

709:     /* resolve the shared entities by exchanging information to adjacent processors */
710:     merr = pcomm->resolve_shared_ents(dmmoab->fileset, cells, dim, dim - 1, NULL, &global_id_tag);MBERRV(mbImpl, merr);
711:     if (dmmoab->fileset) {
712:       merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false, &dmmoab->fileset);MBERRV(mbImpl, merr);
713:     }
714:     else {
715:       merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false);MBERRV(mbImpl, merr);
716:     }

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

722:     PetscLogEventEnd(genCtx.parResolve, 0, 0, 0, 0);
723:   }

725:   if (!genCtx.keep_skins) { // default is to delete the 1- and 2-dimensional entities
726:     // delete all quads and edges
727:     moab::Range toDelete;
728:     if (genCtx.dim > 1) {
729:       merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 1, toDelete);MBERR("Can't get edges", merr);
730:     }

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

736: #ifdef MOAB_HAVE_MPI
737:     merr = dmmoab->pcomm->delete_entities(toDelete) ;MBERR("Can't delete entities", merr);
738: #endif
739:   }

741:   /* set geometric dimension tag for regions */
742:   merr = mbImpl->tag_set_data(geom_tag, &regionset, 1, &dmmoab->dim);MBERRNM(merr);
743:   /* set default material ID for regions */
744:   int default_material = 1;
745:   merr = mbImpl->tag_set_data(mat_tag, &regionset, 1, &default_material);MBERRNM(merr);
746:   /*
747:     int default_dbc = 0;
748:     merr = mbImpl->tag_set_data(dir_tag, &vertexset, 1, &default_dbc);MBERRNM(merr);
749:   */
750:   return 0;
751: }

753: 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)
754: {
755:   char *ropts;
756:   char  ropts_par[PETSC_MAX_PATH_LEN], ropts_pargh[PETSC_MAX_PATH_LEN];
757:   char  ropts_dbg[PETSC_MAX_PATH_LEN];

759:   PetscMalloc1(PETSC_MAX_PATH_LEN, &ropts);
760:   PetscMemzero(&ropts_par, PETSC_MAX_PATH_LEN);
761:   PetscMemzero(&ropts_pargh, PETSC_MAX_PATH_LEN);
762:   PetscMemzero(&ropts_dbg, PETSC_MAX_PATH_LEN);

764:   /* do parallel read unless using only one processor */
765:   if (numproc > 1) {
766:     // 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":""));
767:     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;" : ""));
768:     if (nghost) {
769:       PetscSNPrintf(ropts_pargh, PETSC_MAX_PATH_LEN, "PARALLEL_GHOSTS=%" PetscInt_FMT ".0.%" PetscInt_FMT ";", dim, nghost);
770:     }
771:   }

773:   if (dbglevel) {
774:     if (numproc > 1) {
775:       PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%" PetscInt_FMT ";DEBUG_PIO=%" PetscInt_FMT ";", dbglevel, dbglevel);
776:     }
777:     else PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%" PetscInt_FMT ";", dbglevel);
778:   }

780:   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 : ""));
781:   *read_opts = ropts;
782:   return 0;
783: }

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

788:   Collective

790:   Input Parameters:
791: + comm - The communicator for the DM object
792: . dim - The spatial dimension
793: . filename - The name of the mesh file to be loaded
794: - usrreadopts - The options string to read a MOAB mesh.

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

798:   Output Parameter:
799: . dm  - The DM object

801:   Level: beginner

803: .seealso: DMSetType(), DMCreate(), DMMoabCreateBoxMesh()
804: @*/
805: PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm, PetscInt dim, PetscInt nghost, const char* filename, const char* usrreadopts, DM *dm)
806: {
807:   moab::ErrorCode     merr;
808:   PetscInt            nprocs;
809:   DM_Moab            *dmmoab;
810:   moab::Interface    *mbiface;
811: #ifdef MOAB_HAVE_MPI
812:   moab::ParallelComm *pcomm;
813: #endif
814:   moab::Range         verts, elems;
815:   const char         *readopts;


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

822:   /* get all the necessary handles from the private DM object */
823:   dmmoab = (DM_Moab*)(*dm)->data;
824:   mbiface = dmmoab->mbiface;
825: #ifdef MOAB_HAVE_MPI
826:   pcomm = dmmoab->pcomm;
827:   nprocs = pcomm->size();
828: #else
829:   nprocs = 1;
830: #endif
831:   /* TODO: Decipher dimension based on the loaded mesh instead of getting from user */
832:   dmmoab->dim = dim;
833:   dmmoab->nghostrings = nghost;
834:   dmmoab->refct = 1;

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

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

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

845:   /* Load the mesh from a file. */
846:   if (dmmoab->fileset) {
847:     merr = mbiface->load_file(filename, &dmmoab->fileset, readopts);MBERRVM(mbiface, "Reading MOAB file failed.", merr);
848:   }
849:   else {
850:     merr = mbiface->load_file(filename, 0, readopts);MBERRVM(mbiface, "Reading MOAB file failed.", merr);
851:   }

853: #ifdef MOAB_HAVE_MPI
854:   /* Reassign global IDs on all entities. */
855:   /* merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, true, true, true);MBERRNM(merr); */
856: #endif

858:   /* load the local vertices */
859:   merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);MBERRNM(merr);
860:   /* load the local elements */
861:   merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);MBERRNM(merr);

863: #ifdef MOAB_HAVE_MPI
864:   /* Everything is set up, now just do a tag exchange to update tags
865:      on all of the ghost vertexes */
866:   merr = pcomm->exchange_tags(dmmoab->ltog_tag, verts);MBERRV(mbiface, merr);
867:   merr = pcomm->exchange_tags(dmmoab->ltog_tag, elems);MBERRV(mbiface, merr);
868:   merr = pcomm->collective_sync_partition();MBERR("Collective sync failed", merr);
869: #endif

871:   PetscInfo(*dm, "MOAB file '%s' was successfully loaded. Found %D vertices and %D elements.\n", filename, verts.size(), elems.size());
872:   PetscFree(readopts);
873:   return 0;
874: }

876: /*@C
877:   DMMoabRenumberMeshEntities - Order and number all entities (vertices->elements) to be contiguously ordered
878:   in parallel

880:   Collective

882:   Input Parameters:
883: . dm  - The DM object

885:   Level: advanced

887: .seealso: DMSetUp(), DMCreate()
888: @*/
889: PetscErrorCode DMMoabRenumberMeshEntities(DM dm)
890: {
891:   moab::Range         verts;


895: #ifdef MOAB_HAVE_MPI
896:   /* Insert new points */
897:   moab::ErrorCode     merr;
898:   merr = ((DM_Moab*) dm->data)->pcomm->assign_global_ids(((DM_Moab*) dm->data)->fileset, 3, 0, false, true, false);MBERRNM(merr);
899: #endif
900:   return 0;
901: }