Actual source code: dmmbutil.cxx

petsc-3.10.5 2019-03-28
Report Typos and Errors
  1:  #include <petsc/private/dmmbimpl.h>
  2:  #include <petsc/private/vecimpl.h>

  4:  #include <petscdmmoab.h>
  5: #include <MBTagConventions.hpp>
  6: #include <moab/ReadUtilIface.hpp>
  7: #include <moab/MergeMesh.hpp>
  8: #include <moab/CN.hpp>

 10: typedef struct {
 11:   // options
 12:   PetscInt  A, B, C, M, N, K, dim;
 13:   PetscInt  blockSizeVertexXYZ[3];              // Number of element blocks per partition
 14:   PetscInt  blockSizeElementXYZ[3];
 15:   PetscReal xyzbounds[6]; // the physical size of the domain
 16:   bool      newMergeMethod, keep_skins, simplex, adjEnts;

 18:   // compute params
 19:   PetscReal dx, dy, dz;
 20:   PetscInt  NX, NY, NZ, nex, ney, nez;
 21:   PetscInt  q, xstride, ystride, zstride;
 22:   PetscBool usrxyzgrid, usrprocgrid, usrrefgrid;
 23:   PetscInt  fraction, remainder, cumfraction;
 24:   PetscLogEvent generateMesh, generateElements, generateVertices, parResolve;

 26: } DMMoabMeshGeneratorCtx;


 29: PetscInt DMMoab_SetTensorElementConnectivity_Private(DMMoabMeshGeneratorCtx& genCtx, PetscInt offset, PetscInt corner, std::vector<PetscInt>& subent_conn, moab::EntityHandle *connectivity)
 30: {
 31:   switch (genCtx.dim) {
 32:   case 1:
 33:     subent_conn.resize(2);
 34:     moab::CN::SubEntityVertexIndices(moab::MBEDGE, 1, 0, subent_conn.data());
 35:     connectivity[offset + subent_conn[0]] = corner;
 36:     connectivity[offset + subent_conn[1]] = corner + 1;
 37:     break;
 38:   case 2:
 39:     subent_conn.resize(4);
 40:     moab::CN::SubEntityVertexIndices(moab::MBQUAD, 2, 0, subent_conn.data());
 41:     connectivity[offset + subent_conn[0]] = corner;
 42:     connectivity[offset + subent_conn[1]] = corner + 1;
 43:     connectivity[offset + subent_conn[2]] = corner + 1 + genCtx.ystride;
 44:     connectivity[offset + subent_conn[3]] = corner + genCtx.ystride;
 45:     break;
 46:   case 3:
 47:   default:
 48:     subent_conn.resize(8);
 49:     moab::CN::SubEntityVertexIndices(moab::MBHEX, 3, 0, subent_conn.data());
 50:     connectivity[offset + subent_conn[0]] = corner;
 51:     connectivity[offset + subent_conn[1]] = corner + 1;
 52:     connectivity[offset + subent_conn[2]] = corner + 1 + genCtx.ystride;
 53:     connectivity[offset + subent_conn[3]] = corner + genCtx.ystride;
 54:     connectivity[offset + subent_conn[4]] = corner + genCtx.zstride;
 55:     connectivity[offset + subent_conn[5]] = corner + 1 + genCtx.zstride;
 56:     connectivity[offset + subent_conn[6]] = corner + 1 + genCtx.ystride + genCtx.zstride;
 57:     connectivity[offset + subent_conn[7]] = corner + genCtx.ystride + genCtx.zstride;
 58:     break;
 59:   }
 60:   return subent_conn.size();
 61: }


 64: PetscInt DMMoab_SetSimplexElementConnectivity_Private(DMMoabMeshGeneratorCtx& genCtx, PetscInt subelem, PetscInt offset, PetscInt corner, std::vector<PetscInt>& subent_conn, moab::EntityHandle *connectivity)
 65: {
 66:   PetscInt A, B, C, D, E, F, G, H, M;
 67:   const PetscInt trigen_opts = 1; /* 1 - Aligned diagonally to right, 2 - Aligned diagonally to left, 3 - 4 elements per quad */
 68:   A = corner;
 69:   B = corner + 1;
 70:   switch (genCtx.dim) {
 71:   case 1:
 72:     subent_conn.resize(2);  /* only linear EDGE supported now */
 73:     moab::CN::SubEntityVertexIndices(moab::MBEDGE, 1, 0, subent_conn.data());
 74:     connectivity[offset + subent_conn[0]] = A;
 75:     connectivity[offset + subent_conn[1]] = B;
 76:     break;
 77:   case 2:
 78:     C = corner + 1 + genCtx.ystride;
 79:     D = corner +     genCtx.ystride;
 80:     M = corner + 0.5; /* technically -- need to modify vertex generation */
 81:     subent_conn.resize(3);  /* only linear TRI supported */
 82:     moab::CN::SubEntityVertexIndices(moab::MBTRI, 2, 0, subent_conn.data());
 83:     if (trigen_opts == 1) {
 84:       if (subelem) { /* 0 1 2 of a QUAD */
 85:         connectivity[offset + subent_conn[0]] = B;
 86:         connectivity[offset + subent_conn[1]] = C;
 87:         connectivity[offset + subent_conn[2]] = A;
 88:       }
 89:       else {        /* 2 3 0 of a QUAD */
 90:         connectivity[offset + subent_conn[0]] = D;
 91:         connectivity[offset + subent_conn[1]] = A;
 92:         connectivity[offset + subent_conn[2]] = C;
 93:       }
 94:     }
 95:     else if (trigen_opts == 2) {
 96:       if (subelem) { /* 0 1 2 of a QUAD */
 97:         connectivity[offset + subent_conn[0]] = A;
 98:         connectivity[offset + subent_conn[1]] = B;
 99:         connectivity[offset + subent_conn[2]] = D;
100:       }
101:       else {        /* 2 3 0 of a QUAD */
102:         connectivity[offset + subent_conn[0]] = C;
103:         connectivity[offset + subent_conn[1]] = D;
104:         connectivity[offset + subent_conn[2]] = B;
105:       }
106:     }
107:     else {
108:       switch (subelem) { /* 0 1 2 of a QUAD */
109:       case 0:
110:         connectivity[offset + subent_conn[0]] = A;
111:         connectivity[offset + subent_conn[1]] = B;
112:         connectivity[offset + subent_conn[2]] = M;
113:         break;
114:       case 1:
115:         connectivity[offset + subent_conn[0]] = B;
116:         connectivity[offset + subent_conn[1]] = C;
117:         connectivity[offset + subent_conn[2]] = M;
118:         break;
119:       case 2:
120:         connectivity[offset + subent_conn[0]] = C;
121:         connectivity[offset + subent_conn[1]] = D;
122:         connectivity[offset + subent_conn[2]] = M;
123:         break;
124:       case 3:
125:         connectivity[offset + subent_conn[0]] = D;
126:         connectivity[offset + subent_conn[1]] = A;
127:         connectivity[offset + subent_conn[2]] = M;
128:         break;
129:       }
130:     }
131:     break;
132:   case 3:
133:   default:
134:     C = corner + 1 + genCtx.ystride;
135:     D = corner +     genCtx.ystride;
136:     E = corner +                      genCtx.zstride;
137:     F = corner + 1 +                  genCtx.zstride;
138:     G = corner + 1 + genCtx.ystride + genCtx.zstride;
139:     H = corner +     genCtx.ystride + genCtx.zstride;
140:     subent_conn.resize(4);  /* only linear TET supported */
141:     moab::CN::SubEntityVertexIndices(moab::MBTET, 3, 0, subent_conn.data());
142:     switch (subelem) {
143:     case 0: /* 4 3 7 6 of a HEX */
144:       connectivity[offset + subent_conn[0]] = E;
145:       connectivity[offset + subent_conn[1]] = D;
146:       connectivity[offset + subent_conn[2]] = H;
147:       connectivity[offset + subent_conn[3]] = G;
148:       break;
149:     case 1: /* 0 1 2 5 of a HEX */
150:       connectivity[offset + subent_conn[0]] = A;
151:       connectivity[offset + subent_conn[1]] = B;
152:       connectivity[offset + subent_conn[2]] = C;
153:       connectivity[offset + subent_conn[3]] = F;
154:       break;
155:     case 2: /* 0 3 4 5 of a HEX */
156:       connectivity[offset + subent_conn[0]] = A;
157:       connectivity[offset + subent_conn[1]] = D;
158:       connectivity[offset + subent_conn[2]] = E;
159:       connectivity[offset + subent_conn[3]] = F;
160:       break;
161:     case 3: /* 2 6 3 5 of a HEX */
162:       connectivity[offset + subent_conn[0]] = C;
163:       connectivity[offset + subent_conn[1]] = G;
164:       connectivity[offset + subent_conn[2]] = D;
165:       connectivity[offset + subent_conn[3]] = F;
166:       break;
167:     case 4: /* 0 2 3 5 of a HEX */
168:       connectivity[offset + subent_conn[0]] = A;
169:       connectivity[offset + subent_conn[1]] = C;
170:       connectivity[offset + subent_conn[2]] = D;
171:       connectivity[offset + subent_conn[3]] = F;
172:       break;
173:     case 5: /* 3 6 4 5 of a HEX */
174:       connectivity[offset + subent_conn[0]] = D;
175:       connectivity[offset + subent_conn[1]] = G;
176:       connectivity[offset + subent_conn[2]] = E;
177:       connectivity[offset + subent_conn[3]] = F;
178:       break;
179:     }
180:     break;
181:   }
182:   return subent_conn.size();
183: }


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


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

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

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

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

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

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

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

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

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

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

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

317:   /* create owned elements requested by genCtx */
318:   for (kk = 0; kk < (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1); kk++) {
319:     for (jj = 0; jj < (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] : 1); jj++) {
320:       for (ii = 0; ii < genCtx.blockSizeElementXYZ[0]; ii++) {

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

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

326:         for (PetscInt j = 0; j < entoffset.second; j++) {
327:           /* The entity handle for the particular element -> if we want to set some tags is
328:              moab::EntityHandle eh = starte + ie + j;
329:           */
330:           gids[ie + j] = 1 + ((xe + ii) + (ye + jj) * genCtx.nex + (ze + kk) * (genCtx.nex * genCtx.ney));
331:           /* gids[ie+j] = ie + j + ((xe + ii) + (ye + jj) * genCtx.nex + (ze + kk) * (genCtx.nex * genCtx.ney)); */
332:           /* gids[ie+j] = 1 + ie; */
333:           /* ie++; */
334:         }

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

349: PetscErrorCode DMMBUtil_InitializeOptions(DMMoabMeshGeneratorCtx& genCtx, PetscInt dim, PetscBool simplex, PetscInt rank, PetscInt nprocs, const PetscReal* bounds, PetscInt nelems)
350: {
352:   /* Initialize all genCtx data */
353:   genCtx.dim = dim;
354:   genCtx.simplex = simplex;
355:   genCtx.newMergeMethod = genCtx.keep_skins = genCtx.adjEnts = true;
356:   /* determine other global quantities for the mesh used for nodes increments */
357:   genCtx.q = 1;
358:   genCtx.fraction = genCtx.remainder = genCtx.cumfraction = 0;

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

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

368:     PetscInfo3(NULL, "Fraction = %D, Remainder = %D, Cumulative fraction = %D\n", genCtx.fraction, genCtx.remainder, genCtx.cumfraction);
369:     switch (genCtx.dim) {
370:     case 1:
371:       genCtx.blockSizeElementXYZ[0] = genCtx.fraction;
372:       genCtx.blockSizeElementXYZ[1] = 1;
373:       genCtx.blockSizeElementXYZ[2] = 1;
374:       break;
375:     case 2:
376:       genCtx.blockSizeElementXYZ[0] = nelems;
377:       genCtx.blockSizeElementXYZ[1] = genCtx.fraction;
378:       genCtx.blockSizeElementXYZ[2] = 1;
379:       break;
380:     case 3:
381:     default:
382:       genCtx.blockSizeElementXYZ[0] = nelems;
383:       genCtx.blockSizeElementXYZ[1] = nelems;
384:       genCtx.blockSizeElementXYZ[2] = genCtx.fraction;
385:       break;
386:     }
387:   }

389:   /* partition only by the largest dimension */
390:   /* Total number of local elements := genCtx.blockSizeElementXYZ[0]*(genCtx.dim>1? genCtx.blockSizeElementXYZ[1]*(genCtx.dim>2 ? genCtx.blockSizeElementXYZ[2]:1) :1); */
391:   if (bounds) {
392:     for (PetscInt i = 0; i < 6; i++)
393:       genCtx.xyzbounds[i] = bounds[i];
394:   }
395:   else {
396:     genCtx.xyzbounds[0] = genCtx.xyzbounds[2] = genCtx.xyzbounds[4] = 0.0;
397:     genCtx.xyzbounds[1] = genCtx.xyzbounds[3] = genCtx.xyzbounds[5] = 1.0;
398:   }

400:   if (!genCtx.usrprocgrid) {
401:     switch (genCtx.dim) {
402:     case 1:
403:       genCtx.M = nprocs;
404:       genCtx.N = genCtx.K = 1;
405:       break;
406:     case 2:
407:       genCtx.N = nprocs;
408:       genCtx.M = genCtx.K = 1;
409:       break;
410:     default:
411:       genCtx.K = nprocs;
412:       genCtx.M = genCtx.N = 1;
413:       break;
414:     }
415:   }

417:   if (!genCtx.usrrefgrid) {
418:     genCtx.A = genCtx.B = genCtx.C = 1;
419:   }

421:   /* more default values */
422:   genCtx.nex = genCtx.ney = genCtx.nez = 0;
423:   genCtx.xstride = genCtx.ystride = genCtx.zstride = 0;
424:   genCtx.NX = genCtx.NY = genCtx.NZ = 0;
425:   genCtx.nex = genCtx.ney = genCtx.nez = 0;
426:   genCtx.blockSizeVertexXYZ[0] = genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 1;

428:   switch (genCtx.dim) {
429:   case 3:
430:     genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1;
431:     genCtx.blockSizeVertexXYZ[1] = genCtx.q * genCtx.blockSizeElementXYZ[1] + 1;
432:     genCtx.blockSizeVertexXYZ[2] = genCtx.q * genCtx.blockSizeElementXYZ[2] + 1;

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

452:     genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0];   /* number of elements in x direction, used for global id on element */
453:     genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (genCtx.nex * genCtx.q); /* distance between 2 nodes in x direction */
454:     genCtx.NX = (genCtx.q * genCtx.nex + 1);
455:     genCtx.xstride = 1;
456:     genCtx.ney = genCtx.N * genCtx.B * genCtx.blockSizeElementXYZ[1];  /* number of elements in y direction  .... */
457:     genCtx.dy = (genCtx.xyzbounds[3] - genCtx.xyzbounds[2]) / (nelems * genCtx.q); /* distance between 2 nodes in y direction */
458:     genCtx.NY = (genCtx.q * genCtx.ney + 1);
459:     genCtx.ystride = genCtx.blockSizeVertexXYZ[0];
460:     break;
461:   case 1:
462:     genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 0;
463:     genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1;

465:     genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0];   /* number of elements in x direction, used for global id on element */
466:     genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (nelems * genCtx.q); /* distance between 2 nodes in x direction */
467:     genCtx.NX = (genCtx.q * genCtx.nex + 1);
468:     genCtx.xstride = 1;
469:     break;
470:   }

472:   /* Lets check for some valid input */
473:   if (genCtx.dim < 1 || genCtx.dim > 3) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid topological dimension specified: %d.\n", genCtx.dim);
474:   if (genCtx.M * genCtx.N * genCtx.K != nprocs) SETERRQ4(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid [m, n, k] data: %d, %d, %d. Product must be equal to global size = %d.\n", genCtx.M, genCtx.N, genCtx.K, nprocs);
475:   /* validate the bounds data */
476:   if (genCtx.xyzbounds[0] >= genCtx.xyzbounds[1]) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "X-dim: Left boundary cannot be greater than right. [%G >= %G]\n", genCtx.xyzbounds[0], genCtx.xyzbounds[1]);
477:   if (genCtx.dim > 1 && (genCtx.xyzbounds[2] >= genCtx.xyzbounds[3])) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Y-dim: Left boundary cannot be greater than right. [%G >= %G]\n", genCtx.xyzbounds[2], genCtx.xyzbounds[3]);
478:   if (genCtx.dim > 2 && (genCtx.xyzbounds[4] >= genCtx.xyzbounds[5])) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Z-dim: Left boundary cannot be greater than right. [%G >= %G]\n", genCtx.xyzbounds[4], genCtx.xyzbounds[5]);

480:   PetscInfo3(NULL, "Local elements:= %d, %d, %d\n", genCtx.blockSizeElementXYZ[0], genCtx.blockSizeElementXYZ[1], genCtx.blockSizeElementXYZ[2]);
481:   PetscInfo3(NULL, "Local vertices:= %d, %d, %d\n", genCtx.blockSizeVertexXYZ[0], genCtx.blockSizeVertexXYZ[1], genCtx.blockSizeVertexXYZ[2]);
482:   PetscInfo3(NULL, "Local blocks/processors := %d, %d, %d\n", genCtx.A, genCtx.B, genCtx.C);
483:   PetscInfo3(NULL, "Local processors := %d, %d, %d\n", genCtx.M, genCtx.N, genCtx.K);
484:   PetscInfo3(NULL, "Local nexyz:= %d, %d, %d\n", genCtx.nex, genCtx.ney, genCtx.nez);
485:   PetscInfo3(NULL, "Local delxyz:= %g, %g, %g\n", genCtx.dx, genCtx.dy, genCtx.dz);
486:   PetscInfo3(NULL, "Local strides:= %d, %d, %d\n", genCtx.xstride, genCtx.ystride, genCtx.zstride);
487:   return(0);
488: }

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

493:   Collective on MPI_Comm

495:   Input Parameters:
496: + comm - The communicator for the DM object
497: . dim - The spatial dimension
498: . bounds - The bounds of the box specified with [x-left, x-right, y-bottom, y-top, z-bottom, z-top] depending on the spatial dimension
499: . nele - The number of discrete elements in each direction
500: . user_nghost - The number of ghosted layers needed in the partitioned mesh

502:   Output Parameter:
503: . dm  - The DM object

505:   Level: beginner

507: .keywords: DM, create
508: .seealso: DMSetType(), DMCreate(), DMMoabLoadFromFile()
509: @*/
510: PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool useSimplex, const PetscReal* bounds, PetscInt nele, PetscInt nghost, DM *dm)
511: {
512:   PetscErrorCode         ierr;
513:   moab::ErrorCode        merr;
514:   PetscInt               a, b, c, n, global_size, global_rank;
515:   DM_Moab               *dmmoab;
516:   moab::Interface       *mbImpl;
517: #ifdef MOAB_HAVE_MPI
518:   moab::ParallelComm    *pcomm;
519: #endif
520:   moab::ReadUtilIface   *readMeshIface;
521:   moab::Range            verts, cells, edges, faces, adj, dim3, dim2;
522:   DMMoabMeshGeneratorCtx genCtx;
523:   const PetscInt         npts = nele + 1;    /* Number of points in every dimension */

525:   moab::Tag              global_id_tag, part_tag, geom_tag, mat_tag, dir_tag, neu_tag;
526:   moab::Range            ownedvtx, ownedelms, localvtxs, localelms;
527:   moab::EntityHandle     regionset;
528:   PetscInt               ml = 0, nl = 0, kl = 0;

531:   if (dim < 1 || dim > 3) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension argument for mesh: dim=[1,3].\n");

533:   PetscLogEventRegister("GenerateMesh", DM_CLASSID,   &genCtx.generateMesh);
534:   PetscLogEventRegister("AddVertices", DM_CLASSID,   &genCtx.generateVertices);
535:   PetscLogEventRegister("AddElements", DM_CLASSID,   &genCtx.generateElements);
536:   PetscLogEventRegister("ParResolve", DM_CLASSID,   &genCtx.parResolve);
537:   PetscLogEventBegin(genCtx.generateMesh, 0, 0, 0, 0);
538:   MPI_Comm_size(comm, &global_size);
539:   /* total number of vertices in all dimensions */
540:   n = pow(npts, dim);

542:   /* do some error checking */
543:   if (n < 2) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be >= 2.\n");
544:   if (global_size > n) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of processors must be less than or equal to number of elements.\n");
545:   if (nghost < 0) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of ghost layers cannot be negative.\n");

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

550:   /* get all the necessary handles from the private DM object */
551:   dmmoab = (DM_Moab*)(*dm)->data;
552:   mbImpl = dmmoab->mbiface;
553: #ifdef MOAB_HAVE_MPI
554:   pcomm = dmmoab->pcomm;
555:   global_rank = pcomm->rank();
556: #else
557:   global_rank = 0;
558:   global_size = 1;
559: #endif
560:   global_id_tag = dmmoab->ltog_tag;
561:   dmmoab->dim = dim;
562:   dmmoab->nghostrings = nghost;
563:   dmmoab->refct = 1;

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

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

571:   genCtx.M = genCtx.N = genCtx.K = 1;
572:   genCtx.A = genCtx.B = genCtx.C = 1;

574:   PetscOptionsBegin(comm, "", "DMMoab Creation Options", "DMMOAB");
575:   /* Handle DMMoab spatial resolution */
576:   PetscOptionsInt("-dmb_grid_x", "Number of grid points in x direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[0], &genCtx.blockSizeElementXYZ[0], &genCtx.usrxyzgrid);
577:   if (dim > 1) {PetscOptionsInt("-dmb_grid_y", "Number of grid points in y direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[1], &genCtx.blockSizeElementXYZ[1], &genCtx.usrxyzgrid);}
578:   if (dim > 2) {PetscOptionsInt("-dmb_grid_z", "Number of grid points in z direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[2], &genCtx.blockSizeElementXYZ[2], &genCtx.usrxyzgrid);}

580:   /* Handle DMMoab parallel distibution */
581:   PetscOptionsInt("-dmb_processors_x", "Number of processors in x direction", "DMMoabSetNumProcs", genCtx.M, &genCtx.M, &genCtx.usrprocgrid);
582:   if (dim > 1) {PetscOptionsInt("-dmb_processors_y", "Number of processors in y direction", "DMMoabSetNumProcs", genCtx.N, &genCtx.N, &genCtx.usrprocgrid);}
583:   if (dim > 2) {PetscOptionsInt("-dmb_processors_z", "Number of processors in z direction", "DMMoabSetNumProcs", genCtx.K, &genCtx.K, &genCtx.usrprocgrid);}

585:   /* Handle DMMoab block refinement */
586:   PetscOptionsInt("-dmb_refine_x", "Number of refinement blocks in x direction", "DMMoabSetRefinement", genCtx.A, &genCtx.A, &genCtx.usrrefgrid);
587:   if (dim > 1) {PetscOptionsInt("-dmb_refine_y", "Number of refinement blocks in y direction", "DMMoabSetRefinement", genCtx.B, &genCtx.B, &genCtx.usrrefgrid);}
588:   if (dim > 2) {PetscOptionsInt("-dmb_refine_z", "Number of refinement blocks in z direction", "DMMoabSetRefinement", genCtx.C, &genCtx.C, &genCtx.usrrefgrid);}
589:   PetscOptionsEnd();

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

593:   //if(nele<nprocs) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The dimensional discretization size should be greater or equal to number of processors: %D < %D",nele,nprocs);

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

597:   /* determine m, n, k for processor rank */
598:   ml = nl = kl = 0;
599:   switch (genCtx.dim) {
600:   case 1:
601:     ml = (genCtx.cumfraction);
602:     break;
603:   case 2:
604:     nl = (genCtx.cumfraction);
605:     break;
606:   default:
607:     kl = (genCtx.cumfraction) / genCtx.q / genCtx.blockSizeElementXYZ[2] / genCtx.C; //genCtx.K
608:     break;
609:   }

611:   /*
612:    * so there are a total of M * A * blockSizeElement elements in x direction (so M * A * blockSizeElement + 1 verts in x direction)
613:    * so there are a total of N * B * blockSizeElement elements in y direction (so N * B * blockSizeElement + 1 verts in y direction)
614:    * so there are a total of K * C * blockSizeElement elements in z direction (so K * C * blockSizeElement + 1 verts in z direction)

616:    * there are ( M * A blockSizeElement )      *  ( N * B * blockSizeElement)      * (K * C * blockSizeElement )    hexas
617:    * there are ( M * A * blockSizeElement + 1) *  ( N * B * blockSizeElement + 1 ) * (K * C * blockSizeElement + 1) vertices
618:    * x is the first dimension that varies
619:    */

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

625:   merr = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, mat_tag);MBERR("Getting Material set Tag handle failed", merr);
626:   merr = mbImpl->tag_get_handle(DIRICHLET_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, dir_tag);MBERR("Getting Dirichlet set Tag handle failed", merr);
627:   merr = mbImpl->tag_get_handle(NEUMANN_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, neu_tag);MBERR("Getting Neumann set Tag handle failed", merr);

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

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

636:   for (a = 0; a < (genCtx.dim > 0 ? genCtx.A : genCtx.A); a++) {
637:     for (b = 0; b < (genCtx.dim > 1 ? genCtx.B : 1); b++) {
638:       for (c = 0; c < (genCtx.dim > 2 ? genCtx.C : 1); c++) {

640:         moab::EntityHandle startv;

642:         PetscLogEventBegin(genCtx.generateVertices, 0, 0, 0, 0);
643:         DMMoab_GenerateVertices_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, verts);
644:         PetscLogEventEnd(genCtx.generateVertices, 0, 0, 0, 0);

646:         PetscLogEventBegin(genCtx.generateElements, 0, 0, 0, 0);
647:         DMMoab_GenerateElements_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, cells);
648:         PetscLogEventEnd(genCtx.generateElements, 0, 0, 0, 0);

650:         PetscInt part_num = 0;
651:         switch (genCtx.dim) {
652:         case 3:
653:           part_num += (c + kl * genCtx.C) * (genCtx.M * genCtx.A * genCtx.N * genCtx.B);
654:         case 2:
655:           part_num += (b + nl * genCtx.B) * (genCtx.M * genCtx.A);
656:         case 1:
657:           part_num += (a + ml * genCtx.A);
658:           break;
659:         }

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

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

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

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

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

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

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

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

704:     if (genCtx.A * genCtx.B * genCtx.C != 1) { //  merge needed
705:       moab::MergeMesh mm(mbImpl);
706:       if (genCtx.newMergeMethod) {
707:         merr = mm.merge_using_integer_tag(verts, global_id_tag);MBERR("Can't merge with GLOBAL_ID tag", merr);
708:       }
709:       else {
710:         merr = mm.merge_entities(cells, 0.0001);MBERR("Can't merge with coordinates", merr);
711:       }
712:     }

714: #ifdef MOAB_HAVE_MPI
715:     /* check the handles */
716:     merr = pcomm->check_all_shared_handles();MBERRV(mbImpl, merr);

718:     /* resolve the shared entities by exchanging information to adjacent processors */
719:     merr = pcomm->resolve_shared_ents(dmmoab->fileset, cells, dim, dim - 1, NULL, &global_id_tag);MBERRV(mbImpl, merr);
720:     if (dmmoab->fileset) {
721:       merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false, &dmmoab->fileset);MBERRV(mbImpl, merr);
722:     }
723:     else {
724:       merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false);MBERRV(mbImpl, merr);
725:     }

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

731:     PetscLogEventEnd(genCtx.parResolve, 0, 0, 0, 0);
732:   }

734:   if (!genCtx.keep_skins) { // default is to delete the 1- and 2-dimensional entities
735:     // delete all quads and edges
736:     moab::Range toDelete;
737:     if (genCtx.dim > 1) {
738:       merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 1, toDelete);MBERR("Can't get edges", merr);
739:     }

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

745: #ifdef MOAB_HAVE_MPI
746:     merr = dmmoab->pcomm->delete_entities(toDelete) ;MBERR("Can't delete entities", merr);
747: #endif
748:   }

750:   /* set geometric dimension tag for regions */
751:   merr = mbImpl->tag_set_data(geom_tag, &regionset, 1, &dmmoab->dim);MBERRNM(merr);
752:   /* set default material ID for regions */
753:   int default_material = 1;
754:   merr = mbImpl->tag_set_data(mat_tag, &regionset, 1, &default_material);MBERRNM(merr);
755:   /*
756:     int default_dbc = 0;
757:     merr = mbImpl->tag_set_data(dir_tag, &vertexset, 1, &default_dbc);MBERRNM(merr);
758:   */
759:   return(0);
760: }


763: 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)
764: {
765:   char           *ropts;
766:   char           ropts_par[PETSC_MAX_PATH_LEN], ropts_pargh[PETSC_MAX_PATH_LEN];
767:   char           ropts_dbg[PETSC_MAX_PATH_LEN];

771:   PetscMalloc1(PETSC_MAX_PATH_LEN, &ropts);
772:   PetscMemzero(&ropts_par, PETSC_MAX_PATH_LEN);
773:   PetscMemzero(&ropts_pargh, PETSC_MAX_PATH_LEN);
774:   PetscMemzero(&ropts_dbg, PETSC_MAX_PATH_LEN);

776:   /* do parallel read unless using only one processor */
777:   if (numproc > 1) {
778:     // 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":""));
779:     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;" : ""));
780:     if (nghost) {
781:       PetscSNPrintf(ropts_pargh, PETSC_MAX_PATH_LEN, "PARALLEL_GHOSTS=%d.0.%d;", dim, nghost);
782:     }
783:   }

785:   if (dbglevel) {
786:     if (numproc > 1) {
787:       PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%d;DEBUG_PIO=%d;", dbglevel, dbglevel);
788:     }
789:     else {
790:       PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%d;", dbglevel);
791:     }
792:   }

794:   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 : ""));
795:   *read_opts = ropts;
796:   return(0);
797: }


800: /*@
801:   DMMoabLoadFromFile - Creates a DM object by loading the mesh from a user specified file.

803:   Collective on MPI_Comm

805:   Input Parameters:
806: + comm - The communicator for the DM object
807: . dim - The spatial dimension
808: . filename - The name of the mesh file to be loaded
809: . usrreadopts - The options string to read a MOAB mesh.
810:   Reference (Parallel Mesh Initialization: http://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html#fivetwo)

812:   Output Parameter:
813: . dm  - The DM object

815:   Level: beginner

817: .keywords: DM, create

819: .seealso: DMSetType(), DMCreate(), DMMoabCreateBoxMesh()
820: @*/
821: PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm, PetscInt dim, PetscInt nghost, const char* filename, const char* usrreadopts, DM *dm)
822: {
823:   moab::ErrorCode     merr;
824:   PetscInt            nprocs;
825:   DM_Moab            *dmmoab;
826:   moab::Interface    *mbiface;
827: #ifdef MOAB_HAVE_MPI
828:   moab::ParallelComm *pcomm;
829: #endif
830:   moab::Range         verts, elems;
831:   const char         *readopts;
832:   PetscErrorCode      ierr;


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

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

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

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

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

863:   /* Load the mesh from a file. */
864:   if (dmmoab->fileset) {
865:     merr = mbiface->load_file(filename, &dmmoab->fileset, readopts);MBERRVM(mbiface, "Reading MOAB file failed.", merr);
866:   }
867:   else {
868:     merr = mbiface->load_file(filename, 0, readopts);MBERRVM(mbiface, "Reading MOAB file failed.", merr);
869:   }

871: #ifdef MOAB_HAVE_MPI
872:   /* Reassign global IDs on all entities. */
873:   /* merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, true, true, true);MBERRNM(merr); */
874: #endif

876:   /* load the local vertices */
877:   merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);MBERRNM(merr);
878:   /* load the local elements */
879:   merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);MBERRNM(merr);

881: #ifdef MOAB_HAVE_MPI
882:   /* Everything is set up, now just do a tag exchange to update tags
883:      on all of the ghost vertexes */
884:   merr = pcomm->exchange_tags(dmmoab->ltog_tag, verts);MBERRV(mbiface, merr);
885:   merr = pcomm->exchange_tags(dmmoab->ltog_tag, elems);MBERRV(mbiface, merr);
886:   merr = pcomm->collective_sync_partition();MBERR("Collective sync failed", merr);
887: #endif

889:   PetscInfo3(*dm, "MOAB file '%s' was successfully loaded. Found %D vertices and %D elements.\n", filename, verts.size(), elems.size());
890:   PetscFree(readopts);
891:   return(0);
892: }


895: /*@
896:   DMMoabRenumberMeshEntities - Order and number all entities (vertices->elements) to be contiguously ordered
897:   in parallel

899:   Collective on MPI_Comm

901:   Input Parameters:
902: + dm  - The DM object

904:   Level: advanced

906: .keywords: DM, reorder

908: .seealso: DMSetUp(), DMCreate()
909: @*/
910: PetscErrorCode DMMoabRenumberMeshEntities(DM dm)
911: {
912:   moab::Range         verts;


917: #ifdef MOAB_HAVE_MPI
918:   /* Insert new points */
919:   moab::ErrorCode     merr;
920:   merr = ((DM_Moab*) dm->data)->pcomm->assign_global_ids(((DM_Moab*) dm->data)->fileset, 3, 0, false, true, false);MBERRNM(merr);
921: #endif
922:   return(0);
923: }