1 #include <petsc/private/dmmbimpl.h> /*I "petscdmmoab.h" I*/ 2 #include <petsc/private/vecimpl.h> 3 4 #include <petscdmmoab.h> 5 #include <MBTagConventions.hpp> 6 #include <moab/ReadUtilIface.hpp> 7 #include <moab/MergeMesh.hpp> 8 #include <moab/CN.hpp> 9 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; 17 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; 26 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 } 60 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 } 177 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 } 197 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; 205 206 PetscFunctionBegin; 207 /* we will generate (q*block+1)^3 vertices, and block^3 hexas; q is 1 for linear, 2 for quadratic 208 * the global id of the vertices will come from m, n, k, a, b, c 209 * x will vary from m*A*q*block + a*q*block to m*A*q*block+(a+1)*q*block etc. 210 */ 211 nnodes = genCtx.blockSizeVertexXYZ[0] * (genCtx.dim > 1 ? genCtx.blockSizeVertexXYZ[1] * (genCtx.dim > 2 ? genCtx.blockSizeVertexXYZ[2] : 1) : 1); 212 PetscCall(PetscMalloc1(nnodes, &gids)); 213 214 merr = iface->get_node_coords(3, nnodes, 0, startv, arrays); 215 MBERR("Can't get node coords.", merr); 216 217 /* will start with the lower corner: */ 218 /* x = ( m * genCtx.A + a) * genCtx.q * genCtx.blockSizeElementXYZ[0]; */ 219 /* y = ( n * genCtx.B + b) * genCtx.q * genCtx.blockSizeElementXYZ[1]; */ 220 /* z = ( k * genCtx.C + c) * genCtx.q * genCtx.blockSizeElementXYZ[2]; */ 221 222 x = (m * genCtx.A + a) * genCtx.q; 223 y = (n * genCtx.B + b) * genCtx.q; 224 z = (k * genCtx.C + c) * genCtx.q; 225 PetscCall(PetscInfo(NULL, "Starting offset for coordinates := %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", x, y, z)); 226 ix = 0; 227 moab::Range verts(startv, startv + nnodes - 1); 228 for (kk = 0; kk < (genCtx.dim > 2 ? genCtx.blockSizeVertexXYZ[2] : 1); kk++) { 229 for (jj = 0; jj < (genCtx.dim > 1 ? genCtx.blockSizeVertexXYZ[1] : 1); jj++) { 230 for (ii = 0; ii < genCtx.blockSizeVertexXYZ[0]; ii++, ix++) { 231 /* set coordinates for the vertices */ 232 arrays[0][ix] = (x + ii) * genCtx.dx + genCtx.xyzbounds[0]; 233 arrays[1][ix] = (y + jj) * genCtx.dy + genCtx.xyzbounds[2]; 234 arrays[2][ix] = (z + kk) * genCtx.dz + genCtx.xyzbounds[4]; 235 PetscCall(PetscInfo(NULL, "Creating vertex with coordinates := %f, %f, %f\n", arrays[0][ix], arrays[1][ix], arrays[2][ix])); 236 237 /* If we want to set some tags on the vertices -> use the following entity handle definition: 238 moab::EntityHandle v = startv + ix; 239 */ 240 /* compute the global ID for vertex */ 241 gids[ix] = 1 + (x + ii) + (y + jj) * genCtx.NX + (z + kk) * (genCtx.NX * genCtx.NY); 242 } 243 } 244 } 245 /* set global ID data on vertices */ 246 mbImpl->tag_set_data(global_id_tag, verts, &gids[0]); 247 verts.swap(uverts); 248 PetscCall(PetscFree(gids)); 249 PetscFunctionReturn(0); 250 } 251 252 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) 253 { 254 moab::ErrorCode merr; 255 PetscInt ix, ie, xe, ye, ze; 256 PetscInt ii, jj, kk, nvperelem; 257 PetscInt simplices_per_tensor[4] = {0, 1, 2, 6}; 258 PetscInt ntensorelems = genCtx.blockSizeElementXYZ[0] * (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] * (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1) : 1); /*pow(genCtx.blockSizeElement,genCtx.dim);*/ 259 PetscInt nelems = ntensorelems; 260 moab::EntityHandle starte; /* connectivity */ 261 moab::EntityHandle *conn; 262 263 PetscFunctionBegin; 264 switch (genCtx.dim) { 265 case 1: 266 nvperelem = 2; 267 merr = iface->get_element_connect(nelems, 2, moab::MBEDGE, 0, starte, conn); 268 MBERR("Can't get EDGE2 element connectivity.", merr); 269 break; 270 case 2: 271 if (genCtx.simplex) { 272 nvperelem = 3; 273 nelems = ntensorelems * simplices_per_tensor[genCtx.dim]; 274 merr = iface->get_element_connect(nelems, 3, moab::MBTRI, 0, starte, conn); 275 MBERR("Can't get TRI3 element connectivity.", merr); 276 } else { 277 nvperelem = 4; 278 merr = iface->get_element_connect(nelems, 4, moab::MBQUAD, 0, starte, conn); 279 MBERR("Can't get QUAD4 element connectivity.", merr); 280 } 281 break; 282 case 3: 283 default: 284 if (genCtx.simplex) { 285 nvperelem = 4; 286 nelems = ntensorelems * simplices_per_tensor[genCtx.dim]; 287 merr = iface->get_element_connect(nelems, 4, moab::MBTET, 0, starte, conn); 288 MBERR("Can't get TET4 element connectivity.", merr); 289 } else { 290 nvperelem = 8; 291 merr = iface->get_element_connect(nelems, 8, moab::MBHEX, 0, starte, conn); 292 MBERR("Can't get HEX8 element connectivity.", merr); 293 } 294 break; 295 } 296 297 ix = ie = 0; /* index now in the elements, for global ids */ 298 299 /* create a temporary range to store local element handles */ 300 moab::Range tmp(starte, starte + nelems - 1); 301 std::vector<PetscInt> gids(nelems); 302 303 /* identify the elements at the lower corner, for their global ids */ 304 xe = m * genCtx.A * genCtx.blockSizeElementXYZ[0] + a * genCtx.blockSizeElementXYZ[0]; 305 ye = (genCtx.dim > 1 ? n * genCtx.B * genCtx.blockSizeElementXYZ[1] + b * genCtx.blockSizeElementXYZ[1] : 0); 306 ze = (genCtx.dim > 2 ? k * genCtx.C * genCtx.blockSizeElementXYZ[2] + c * genCtx.blockSizeElementXYZ[2] : 0); 307 308 /* create owned elements requested by genCtx */ 309 for (kk = 0; kk < (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1); kk++) { 310 for (jj = 0; jj < (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] : 1); jj++) { 311 for (ii = 0; ii < genCtx.blockSizeElementXYZ[0]; ii++) { 312 moab::EntityHandle corner = startv + genCtx.q * ii + genCtx.q * jj * genCtx.ystride + genCtx.q * kk * genCtx.zstride; 313 314 std::pair<PetscInt, PetscInt> entoffset = DMMoab_SetElementConnectivity_Private(genCtx, ix, corner, conn); 315 316 for (PetscInt j = 0; j < entoffset.second; j++) { 317 /* The entity handle for the particular element -> if we want to set some tags is 318 moab::EntityHandle eh = starte + ie + j; 319 */ 320 gids[ie + j] = 1 + ((xe + ii) + (ye + jj) * genCtx.nex + (ze + kk) * (genCtx.nex * genCtx.ney)); 321 /* gids[ie+j] = ie + j + ((xe + ii) + (ye + jj) * genCtx.nex + (ze + kk) * (genCtx.nex * genCtx.ney)); */ 322 /* gids[ie+j] = 1 + ie; */ 323 /* ie++; */ 324 } 325 326 ix += entoffset.first; 327 ie += entoffset.second; 328 } 329 } 330 } 331 if (genCtx.adjEnts) { /* we need to update adjacencies now, because some elements are new */ 332 merr = iface->update_adjacencies(starte, nelems, nvperelem, conn); 333 MBERR("Can't update adjacencies", merr); 334 } 335 tmp.swap(cells); 336 merr = mbImpl->tag_set_data(global_id_tag, cells, &gids[0]); 337 MBERR("Can't set global ids to elements.", merr); 338 PetscFunctionReturn(0); 339 } 340 341 PetscErrorCode DMMBUtil_InitializeOptions(DMMoabMeshGeneratorCtx &genCtx, PetscInt dim, PetscBool simplex, PetscInt rank, PetscInt nprocs, const PetscReal *bounds, PetscInt nelems) 342 { 343 PetscFunctionBegin; 344 /* Initialize all genCtx data */ 345 genCtx.dim = dim; 346 genCtx.simplex = simplex; 347 genCtx.newMergeMethod = genCtx.keep_skins = genCtx.adjEnts = true; 348 /* determine other global quantities for the mesh used for nodes increments */ 349 genCtx.q = 1; 350 genCtx.fraction = genCtx.remainder = genCtx.cumfraction = 0; 351 352 if (!genCtx.usrxyzgrid) { /* not overridden by genCtx - assume nele equally and that genCtx wants a uniform cube mesh */ 353 354 genCtx.fraction = nelems / nprocs; /* partition only by the largest dimension */ 355 genCtx.remainder = nelems % nprocs; /* remainder after partition which gets evenly distributed by round-robin */ 356 genCtx.cumfraction = (rank > 0 ? (genCtx.fraction) * (rank) + (rank - 1 < genCtx.remainder ? rank : genCtx.remainder) : 0); 357 if (rank < genCtx.remainder) /* This process gets "fraction+1" elements */ 358 genCtx.fraction++; 359 360 PetscCall(PetscInfo(NULL, "Fraction = %" PetscInt_FMT ", Remainder = %" PetscInt_FMT ", Cumulative fraction = %" PetscInt_FMT "\n", genCtx.fraction, genCtx.remainder, genCtx.cumfraction)); 361 switch (genCtx.dim) { 362 case 1: 363 genCtx.blockSizeElementXYZ[0] = genCtx.fraction; 364 genCtx.blockSizeElementXYZ[1] = 1; 365 genCtx.blockSizeElementXYZ[2] = 1; 366 break; 367 case 2: 368 genCtx.blockSizeElementXYZ[0] = nelems; 369 genCtx.blockSizeElementXYZ[1] = genCtx.fraction; 370 genCtx.blockSizeElementXYZ[2] = 1; 371 break; 372 case 3: 373 default: 374 genCtx.blockSizeElementXYZ[0] = nelems; 375 genCtx.blockSizeElementXYZ[1] = nelems; 376 genCtx.blockSizeElementXYZ[2] = genCtx.fraction; 377 break; 378 } 379 } 380 381 /* partition only by the largest dimension */ 382 /* Total number of local elements := genCtx.blockSizeElementXYZ[0]*(genCtx.dim>1? genCtx.blockSizeElementXYZ[1]*(genCtx.dim>2 ? genCtx.blockSizeElementXYZ[2]:1) :1); */ 383 if (bounds) { 384 for (PetscInt i = 0; i < 6; i++) genCtx.xyzbounds[i] = bounds[i]; 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 } 389 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 } 406 407 if (!genCtx.usrrefgrid) genCtx.A = genCtx.B = genCtx.C = 1; 408 409 /* more default values */ 410 genCtx.nex = genCtx.ney = genCtx.nez = 0; 411 genCtx.xstride = genCtx.ystride = genCtx.zstride = 0; 412 genCtx.NX = genCtx.NY = genCtx.NZ = 0; 413 genCtx.nex = genCtx.ney = genCtx.nez = 0; 414 genCtx.blockSizeVertexXYZ[0] = genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 1; 415 416 switch (genCtx.dim) { 417 case 3: 418 genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1; 419 genCtx.blockSizeVertexXYZ[1] = genCtx.q * genCtx.blockSizeElementXYZ[1] + 1; 420 genCtx.blockSizeVertexXYZ[2] = genCtx.q * genCtx.blockSizeElementXYZ[2] + 1; 421 422 genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0]; /* number of elements in x direction, used for global id on element */ 423 genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (nelems * genCtx.q); /* distance between 2 nodes in x direction */ 424 genCtx.NX = (genCtx.q * genCtx.nex + 1); 425 genCtx.xstride = 1; 426 genCtx.ney = genCtx.N * genCtx.B * genCtx.blockSizeElementXYZ[1]; /* number of elements in y direction .... */ 427 genCtx.dy = (genCtx.xyzbounds[3] - genCtx.xyzbounds[2]) / (nelems * genCtx.q); /* distance between 2 nodes in y direction */ 428 genCtx.NY = (genCtx.q * genCtx.ney + 1); 429 genCtx.ystride = genCtx.blockSizeVertexXYZ[0]; 430 genCtx.nez = genCtx.K * genCtx.C * genCtx.blockSizeElementXYZ[2]; /* number of elements in z direction .... */ 431 genCtx.dz = (genCtx.xyzbounds[5] - genCtx.xyzbounds[4]) / (nelems * genCtx.q); /* distance between 2 nodes in z direction */ 432 genCtx.NZ = (genCtx.q * genCtx.nez + 1); 433 genCtx.zstride = genCtx.blockSizeVertexXYZ[0] * genCtx.blockSizeVertexXYZ[1]; 434 break; 435 case 2: 436 genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1; 437 genCtx.blockSizeVertexXYZ[1] = genCtx.q * genCtx.blockSizeElementXYZ[1] + 1; 438 genCtx.blockSizeVertexXYZ[2] = 0; 439 440 genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0]; /* number of elements in x direction, used for global id on element */ 441 genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (genCtx.nex * genCtx.q); /* distance between 2 nodes in x direction */ 442 genCtx.NX = (genCtx.q * genCtx.nex + 1); 443 genCtx.xstride = 1; 444 genCtx.ney = genCtx.N * genCtx.B * genCtx.blockSizeElementXYZ[1]; /* number of elements in y direction .... */ 445 genCtx.dy = (genCtx.xyzbounds[3] - genCtx.xyzbounds[2]) / (nelems * genCtx.q); /* distance between 2 nodes in y direction */ 446 genCtx.NY = (genCtx.q * genCtx.ney + 1); 447 genCtx.ystride = genCtx.blockSizeVertexXYZ[0]; 448 break; 449 case 1: 450 genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 0; 451 genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1; 452 453 genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0]; /* number of elements in x direction, used for global id on element */ 454 genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (nelems * genCtx.q); /* distance between 2 nodes in x direction */ 455 genCtx.NX = (genCtx.q * genCtx.nex + 1); 456 genCtx.xstride = 1; 457 break; 458 } 459 460 /* Lets check for some valid input */ 461 PetscCheck(genCtx.dim >= 1 && genCtx.dim <= 3, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid topological dimension specified: %" PetscInt_FMT ".", genCtx.dim); 462 PetscCheck(genCtx.M * genCtx.N * genCtx.K == nprocs, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid [m, n, k] data: %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ". Product must be equal to global size = %" PetscInt_FMT ".", genCtx.M, 463 genCtx.N, genCtx.K, nprocs); 464 /* validate the bounds data */ 465 PetscCheck(genCtx.xyzbounds[0] < genCtx.xyzbounds[1], PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "X-dim: Left boundary cannot be greater than right. [%G >= %G]", genCtx.xyzbounds[0], genCtx.xyzbounds[1]); 466 PetscCheck(genCtx.dim <= 1 || genCtx.xyzbounds[2] < genCtx.xyzbounds[3], PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Y-dim: Left boundary cannot be greater than right. [%G >= %G]", genCtx.xyzbounds[2], genCtx.xyzbounds[3]); 467 PetscCheck(genCtx.dim <= 2 || genCtx.xyzbounds[4] < genCtx.xyzbounds[5], PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Z-dim: Left boundary cannot be greater than right. [%G >= %G]", genCtx.xyzbounds[4], genCtx.xyzbounds[5]); 468 469 PetscCall(PetscInfo(NULL, "Local elements:= %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.blockSizeElementXYZ[0], genCtx.blockSizeElementXYZ[1], genCtx.blockSizeElementXYZ[2])); 470 PetscCall(PetscInfo(NULL, "Local vertices:= %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.blockSizeVertexXYZ[0], genCtx.blockSizeVertexXYZ[1], genCtx.blockSizeVertexXYZ[2])); 471 PetscCall(PetscInfo(NULL, "Local blocks/processors := %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.A, genCtx.B, genCtx.C)); 472 PetscCall(PetscInfo(NULL, "Local processors := %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.M, genCtx.N, genCtx.K)); 473 PetscCall(PetscInfo(NULL, "Local nexyz:= %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.nex, genCtx.ney, genCtx.nez)); 474 PetscCall(PetscInfo(NULL, "Local delxyz:= %g, %g, %g\n", genCtx.dx, genCtx.dy, genCtx.dz)); 475 PetscCall(PetscInfo(NULL, "Local strides:= %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT "\n", genCtx.xstride, genCtx.ystride, genCtx.zstride)); 476 PetscFunctionReturn(0); 477 } 478 479 /*@C 480 DMMoabCreateBoxMesh - Creates a mesh on the tensor product (box) of intervals with genCtx specified bounds. 481 482 Collective 483 484 Input Parameters: 485 + comm - The communicator for the DM object 486 . dim - The spatial dimension 487 . 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 488 . nele - The number of discrete elements in each direction 489 - user_nghost - The number of ghosted layers needed in the partitioned mesh 490 491 Output Parameter: 492 . dm - The DM object 493 494 Level: beginner 495 496 .seealso: `DMSetType()`, `DMCreate()`, `DMMoabLoadFromFile()` 497 @*/ 498 PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool useSimplex, const PetscReal *bounds, PetscInt nele, PetscInt nghost, DM *dm) 499 { 500 moab::ErrorCode merr; 501 PetscInt a, b, c, n, global_size, global_rank; 502 DM_Moab *dmmoab; 503 moab::Interface *mbImpl; 504 #ifdef MOAB_HAVE_MPI 505 moab::ParallelComm *pcomm; 506 #endif 507 moab::ReadUtilIface *readMeshIface; 508 moab::Range verts, cells, edges, faces, adj, dim3, dim2; 509 DMMoabMeshGeneratorCtx genCtx; 510 const PetscInt npts = nele + 1; /* Number of points in every dimension */ 511 512 moab::Tag global_id_tag, part_tag, geom_tag, mat_tag, dir_tag, neu_tag; 513 moab::Range ownedvtx, ownedelms, localvtxs, localelms; 514 moab::EntityHandle regionset; 515 PetscInt ml = 0, nl = 0, kl = 0; 516 517 PetscFunctionBegin; 518 PetscCheck(dim >= 1 && dim <= 3, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension argument for mesh: dim=[1,3]."); 519 520 PetscCall(PetscLogEventRegister("GenerateMesh", DM_CLASSID, &genCtx.generateMesh)); 521 PetscCall(PetscLogEventRegister("AddVertices", DM_CLASSID, &genCtx.generateVertices)); 522 PetscCall(PetscLogEventRegister("AddElements", DM_CLASSID, &genCtx.generateElements)); 523 PetscCall(PetscLogEventRegister("ParResolve", DM_CLASSID, &genCtx.parResolve)); 524 PetscCall(PetscLogEventBegin(genCtx.generateMesh, 0, 0, 0, 0)); 525 PetscCallMPI(MPI_Comm_size(comm, &global_size)); 526 /* total number of vertices in all dimensions */ 527 n = pow(npts, dim); 528 529 /* do some error checking */ 530 PetscCheck(n >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be >= 2."); 531 PetscCheck(global_size <= n, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of processors must be less than or equal to number of elements."); 532 PetscCheck(nghost >= 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of ghost layers cannot be negative."); 533 534 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 535 PetscCall(DMMoabCreateMoab(comm, NULL, NULL, NULL, dm)); 536 537 /* get all the necessary handles from the private DM object */ 538 dmmoab = (DM_Moab *)(*dm)->data; 539 mbImpl = dmmoab->mbiface; 540 #ifdef MOAB_HAVE_MPI 541 pcomm = dmmoab->pcomm; 542 global_rank = pcomm->rank(); 543 #else 544 global_rank = 0; 545 global_size = 1; 546 #endif 547 global_id_tag = dmmoab->ltog_tag; 548 dmmoab->dim = dim; 549 dmmoab->nghostrings = nghost; 550 dmmoab->refct = 1; 551 552 /* create a file set to associate all entities in current mesh */ 553 merr = mbImpl->create_meshset(moab::MESHSET_SET, dmmoab->fileset); 554 MBERR("Creating file set failed", merr); 555 556 /* No errors yet; proceed with building the mesh */ 557 merr = mbImpl->query_interface(readMeshIface); 558 MBERRNM(merr); 559 560 genCtx.M = genCtx.N = genCtx.K = 1; 561 genCtx.A = genCtx.B = genCtx.C = 1; 562 genCtx.blockSizeElementXYZ[0] = 0; 563 genCtx.blockSizeElementXYZ[1] = 0; 564 genCtx.blockSizeElementXYZ[2] = 0; 565 566 PetscOptionsBegin(comm, "", "DMMoab Creation Options", "DMMOAB"); 567 /* Handle DMMoab spatial resolution */ 568 PetscCall(PetscOptionsInt("-dmb_grid_x", "Number of grid points in x direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[0], &genCtx.blockSizeElementXYZ[0], &genCtx.usrxyzgrid)); 569 if (dim > 1) PetscCall(PetscOptionsInt("-dmb_grid_y", "Number of grid points in y direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[1], &genCtx.blockSizeElementXYZ[1], &genCtx.usrxyzgrid)); 570 if (dim > 2) PetscCall(PetscOptionsInt("-dmb_grid_z", "Number of grid points in z direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[2], &genCtx.blockSizeElementXYZ[2], &genCtx.usrxyzgrid)); 571 572 /* Handle DMMoab parallel distribution */ 573 PetscCall(PetscOptionsInt("-dmb_processors_x", "Number of processors in x direction", "DMMoabSetNumProcs", genCtx.M, &genCtx.M, &genCtx.usrprocgrid)); 574 if (dim > 1) PetscCall(PetscOptionsInt("-dmb_processors_y", "Number of processors in y direction", "DMMoabSetNumProcs", genCtx.N, &genCtx.N, &genCtx.usrprocgrid)); 575 if (dim > 2) PetscCall(PetscOptionsInt("-dmb_processors_z", "Number of processors in z direction", "DMMoabSetNumProcs", genCtx.K, &genCtx.K, &genCtx.usrprocgrid)); 576 577 /* Handle DMMoab block refinement */ 578 PetscCall(PetscOptionsInt("-dmb_refine_x", "Number of refinement blocks in x direction", "DMMoabSetRefinement", genCtx.A, &genCtx.A, &genCtx.usrrefgrid)); 579 if (dim > 1) PetscCall(PetscOptionsInt("-dmb_refine_y", "Number of refinement blocks in y direction", "DMMoabSetRefinement", genCtx.B, &genCtx.B, &genCtx.usrrefgrid)); 580 if (dim > 2) PetscCall(PetscOptionsInt("-dmb_refine_z", "Number of refinement blocks in z direction", "DMMoabSetRefinement", genCtx.C, &genCtx.C, &genCtx.usrrefgrid)); 581 PetscOptionsEnd(); 582 583 PetscCall(DMMBUtil_InitializeOptions(genCtx, dim, useSimplex, global_rank, global_size, bounds, nele)); 584 585 //PetscCheck(nele>=nprocs,PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The dimensional discretization size should be greater or equal to number of processors: %" PetscInt_FMT " < %" PetscInt_FMT,nele,nprocs); 586 587 if (genCtx.adjEnts) genCtx.keep_skins = true; /* do not delete anything - consumes more memory */ 588 589 /* determine m, n, k for processor rank */ 590 ml = nl = kl = 0; 591 switch (genCtx.dim) { 592 case 1: 593 ml = (genCtx.cumfraction); 594 break; 595 case 2: 596 nl = (genCtx.cumfraction); 597 break; 598 default: 599 kl = (genCtx.cumfraction) / genCtx.q / genCtx.blockSizeElementXYZ[2] / genCtx.C; //genCtx.K 600 break; 601 } 602 603 /* 604 * so there are a total of M * A * blockSizeElement elements in x direction (so M * A * blockSizeElement + 1 verts in x direction) 605 * so there are a total of N * B * blockSizeElement elements in y direction (so N * B * blockSizeElement + 1 verts in y direction) 606 * so there are a total of K * C * blockSizeElement elements in z direction (so K * C * blockSizeElement + 1 verts in z direction) 607 608 * there are ( M * A blockSizeElement) * ( N * B * blockSizeElement) * (K * C * blockSizeElement) hexas 609 * there are ( M * A * blockSizeElement + 1) * ( N * B * blockSizeElement + 1) * (K * C * blockSizeElement + 1) vertices 610 * x is the first dimension that varies 611 */ 612 613 /* generate the block at (a, b, c); it will represent a partition , it will get a partition tag */ 614 PetscInt dum_id = -1; 615 merr = mbImpl->tag_get_handle("GLOBAL_ID", 1, moab::MB_TYPE_INTEGER, global_id_tag); 616 MBERR("Getting Global_ID Tag handle failed", merr); 617 618 merr = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, mat_tag); 619 MBERR("Getting Material set Tag handle failed", merr); 620 merr = mbImpl->tag_get_handle(DIRICHLET_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, dir_tag); 621 MBERR("Getting Dirichlet set Tag handle failed", merr); 622 merr = mbImpl->tag_get_handle(NEUMANN_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, neu_tag); 623 MBERR("Getting Neumann set Tag handle failed", merr); 624 625 merr = mbImpl->tag_get_handle("PARALLEL_PARTITION", 1, moab::MB_TYPE_INTEGER, part_tag, moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &dum_id); 626 MBERR("Getting Partition Tag handle failed", merr); 627 628 /* lets create some sets */ 629 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); 630 MBERRNM(merr); 631 merr = mbImpl->create_meshset(moab::MESHSET_SET, regionset); 632 MBERRNM(merr); 633 PetscCall(PetscLogEventEnd(genCtx.generateMesh, 0, 0, 0, 0)); 634 635 for (a = 0; a < (genCtx.dim > 0 ? genCtx.A : genCtx.A); a++) { 636 for (b = 0; b < (genCtx.dim > 1 ? genCtx.B : 1); b++) { 637 for (c = 0; c < (genCtx.dim > 2 ? genCtx.C : 1); c++) { 638 moab::EntityHandle startv; 639 640 PetscCall(PetscLogEventBegin(genCtx.generateVertices, 0, 0, 0, 0)); 641 PetscCall(DMMoab_GenerateVertices_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, verts)); 642 PetscCall(PetscLogEventEnd(genCtx.generateVertices, 0, 0, 0, 0)); 643 644 PetscCall(PetscLogEventBegin(genCtx.generateElements, 0, 0, 0, 0)); 645 PetscCall(DMMoab_GenerateElements_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, cells)); 646 PetscCall(PetscLogEventEnd(genCtx.generateElements, 0, 0, 0, 0)); 647 648 PetscInt part_num = 0; 649 switch (genCtx.dim) { 650 case 3: 651 part_num += (c + kl * genCtx.C) * (genCtx.M * genCtx.A * genCtx.N * genCtx.B); 652 case 2: 653 part_num += (b + nl * genCtx.B) * (genCtx.M * genCtx.A); 654 case 1: 655 part_num += (a + ml * genCtx.A); 656 break; 657 } 658 659 moab::EntityHandle part_set; 660 merr = mbImpl->create_meshset(moab::MESHSET_SET, part_set); 661 MBERR("Can't create mesh set.", merr); 662 663 merr = mbImpl->add_entities(part_set, verts); 664 MBERR("Can't add vertices to set.", merr); 665 merr = mbImpl->add_entities(part_set, cells); 666 MBERR("Can't add entities to set.", merr); 667 merr = mbImpl->add_entities(regionset, cells); 668 MBERR("Can't add entities to set.", merr); 669 670 /* if needed, add all edges and faces */ 671 if (genCtx.adjEnts) { 672 if (genCtx.dim > 1) { 673 merr = mbImpl->get_adjacencies(cells, 1, true, edges, moab::Interface::UNION); 674 MBERR("Can't get edges", merr); 675 merr = mbImpl->add_entities(part_set, edges); 676 MBERR("Can't add edges to partition set.", merr); 677 } 678 if (genCtx.dim > 2) { 679 merr = mbImpl->get_adjacencies(cells, 2, true, faces, moab::Interface::UNION); 680 MBERR("Can't get faces", merr); 681 merr = mbImpl->add_entities(part_set, faces); 682 MBERR("Can't add faces to partition set.", merr); 683 } 684 edges.clear(); 685 faces.clear(); 686 } 687 verts.clear(); 688 cells.clear(); 689 690 merr = mbImpl->tag_set_data(part_tag, &part_set, 1, &part_num); 691 MBERR("Can't set part tag on set", merr); 692 if (dmmoab->fileset) { 693 merr = mbImpl->add_parent_child(dmmoab->fileset, part_set); 694 MBERR("Can't add part set to file set.", merr); 695 merr = mbImpl->unite_meshset(dmmoab->fileset, part_set); 696 MBERRNM(merr); 697 } 698 merr = mbImpl->add_entities(dmmoab->fileset, &part_set, 1); 699 MBERRNM(merr); 700 } 701 } 702 } 703 704 merr = mbImpl->add_parent_child(dmmoab->fileset, regionset); 705 MBERRNM(merr); 706 707 /* Only in parallel: resolve shared entities between processors and exchange ghost layers */ 708 if (global_size > 1) { 709 PetscCall(PetscLogEventBegin(genCtx.parResolve, 0, 0, 0, 0)); 710 711 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, genCtx.dim, cells); 712 MBERR("Can't get all d-dimensional elements.", merr); 713 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 0, verts); 714 MBERR("Can't get all vertices.", merr); 715 716 if (genCtx.A * genCtx.B * genCtx.C != 1) { // merge needed 717 moab::MergeMesh mm(mbImpl); 718 if (genCtx.newMergeMethod) { 719 merr = mm.merge_using_integer_tag(verts, global_id_tag); 720 MBERR("Can't merge with GLOBAL_ID tag", merr); 721 } else { 722 merr = mm.merge_entities(cells, 0.0001); 723 MBERR("Can't merge with coordinates", merr); 724 } 725 } 726 727 #ifdef MOAB_HAVE_MPI 728 /* check the handles */ 729 merr = pcomm->check_all_shared_handles(); 730 MBERRV(mbImpl, merr); 731 732 /* resolve the shared entities by exchanging information to adjacent processors */ 733 merr = pcomm->resolve_shared_ents(dmmoab->fileset, cells, dim, dim - 1, NULL, &global_id_tag); 734 MBERRV(mbImpl, merr); 735 if (dmmoab->fileset) { 736 merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false, &dmmoab->fileset); 737 MBERRV(mbImpl, merr); 738 } else { 739 merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false); 740 MBERRV(mbImpl, merr); 741 } 742 743 /* Reassign global IDs on all entities. */ 744 merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, false, true, false); 745 MBERRNM(merr); 746 #endif 747 748 PetscCall(PetscLogEventEnd(genCtx.parResolve, 0, 0, 0, 0)); 749 } 750 751 if (!genCtx.keep_skins) { // default is to delete the 1- and 2-dimensional entities 752 // delete all quads and edges 753 moab::Range toDelete; 754 if (genCtx.dim > 1) { 755 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 1, toDelete); 756 MBERR("Can't get edges", merr); 757 } 758 759 if (genCtx.dim > 2) { 760 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 2, toDelete); 761 MBERR("Can't get faces", merr); 762 } 763 764 #ifdef MOAB_HAVE_MPI 765 merr = dmmoab->pcomm->delete_entities(toDelete); 766 MBERR("Can't delete entities", merr); 767 #endif 768 } 769 770 /* set geometric dimension tag for regions */ 771 merr = mbImpl->tag_set_data(geom_tag, ®ionset, 1, &dmmoab->dim); 772 MBERRNM(merr); 773 /* set default material ID for regions */ 774 int default_material = 1; 775 merr = mbImpl->tag_set_data(mat_tag, ®ionset, 1, &default_material); 776 MBERRNM(merr); 777 /* 778 int default_dbc = 0; 779 merr = mbImpl->tag_set_data(dir_tag, &vertexset, 1, &default_dbc);MBERRNM(merr); 780 */ 781 PetscFunctionReturn(0); 782 } 783 784 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) 785 { 786 char *ropts; 787 char ropts_par[PETSC_MAX_PATH_LEN], ropts_pargh[PETSC_MAX_PATH_LEN]; 788 char ropts_dbg[PETSC_MAX_PATH_LEN]; 789 790 PetscFunctionBegin; 791 PetscCall(PetscMalloc1(PETSC_MAX_PATH_LEN, &ropts)); 792 PetscCall(PetscMemzero(&ropts_par, PETSC_MAX_PATH_LEN)); 793 PetscCall(PetscMemzero(&ropts_pargh, PETSC_MAX_PATH_LEN)); 794 PetscCall(PetscMemzero(&ropts_dbg, PETSC_MAX_PATH_LEN)); 795 796 /* do parallel read unless using only one processor */ 797 if (numproc > 1) { 798 // PetscCall(PetscSNPrintf(ropts_par, PETSC_MAX_PATH_LEN, "PARALLEL=%s;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=%d.0.1%s;",MoabReadModes[mode],dim,(by_rank ? ";PARTITION_BY_RANK":""))); 799 PetscCall(PetscSNPrintf(ropts_par, PETSC_MAX_PATH_LEN, "PARALLEL=%s;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS;%s", MoabReadModes[mode], (by_rank ? "PARTITION_BY_RANK;" : ""))); 800 if (nghost) PetscCall(PetscSNPrintf(ropts_pargh, PETSC_MAX_PATH_LEN, "PARALLEL_GHOSTS=%" PetscInt_FMT ".0.%" PetscInt_FMT ";", dim, nghost)); 801 } 802 803 if (dbglevel) { 804 if (numproc > 1) { 805 PetscCall(PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%" PetscInt_FMT ";DEBUG_PIO=%" PetscInt_FMT ";", dbglevel, dbglevel)); 806 } else PetscCall(PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%" PetscInt_FMT ";", dbglevel)); 807 } 808 809 PetscCall(PetscSNPrintf(ropts, PETSC_MAX_PATH_LEN, "%s%s%s%s%s", ropts_par, (nghost ? ropts_pargh : ""), ropts_dbg, (extra_opts ? extra_opts : ""), (dm_opts ? dm_opts : ""))); 810 *read_opts = ropts; 811 PetscFunctionReturn(0); 812 } 813 814 /*@C 815 DMMoabLoadFromFile - Creates a DM object by loading the mesh from a user specified file. 816 817 Collective 818 819 Input Parameters: 820 + comm - The communicator for the DM object 821 . dim - The spatial dimension 822 . filename - The name of the mesh file to be loaded 823 - usrreadopts - The options string to read a MOAB mesh. 824 825 Reference (Parallel Mesh Initialization: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html#fivetwo) 826 827 Output Parameter: 828 . dm - The DM object 829 830 Level: beginner 831 832 .seealso: `DMSetType()`, `DMCreate()`, `DMMoabCreateBoxMesh()` 833 @*/ 834 PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm, PetscInt dim, PetscInt nghost, const char *filename, const char *usrreadopts, DM *dm) 835 { 836 moab::ErrorCode merr; 837 PetscInt nprocs; 838 DM_Moab *dmmoab; 839 moab::Interface *mbiface; 840 #ifdef MOAB_HAVE_MPI 841 moab::ParallelComm *pcomm; 842 #endif 843 moab::Range verts, elems; 844 const char *readopts; 845 846 PetscFunctionBegin; 847 PetscValidPointer(dm, 6); 848 849 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 850 PetscCall(DMMoabCreateMoab(comm, NULL, NULL, NULL, dm)); 851 852 /* get all the necessary handles from the private DM object */ 853 dmmoab = (DM_Moab *)(*dm)->data; 854 mbiface = dmmoab->mbiface; 855 #ifdef MOAB_HAVE_MPI 856 pcomm = dmmoab->pcomm; 857 nprocs = pcomm->size(); 858 #else 859 nprocs = 1; 860 #endif 861 /* TODO: Decipher dimension based on the loaded mesh instead of getting from user */ 862 dmmoab->dim = dim; 863 dmmoab->nghostrings = nghost; 864 dmmoab->refct = 1; 865 866 /* create a file set to associate all entities in current mesh */ 867 merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, dmmoab->fileset); 868 MBERR("Creating file set failed", merr); 869 870 /* add mesh loading options specific to the DM */ 871 PetscCall(DMMoab_GetReadOptions_Private(dmmoab->partition_by_rank, nprocs, dim, nghost, dmmoab->read_mode, dmmoab->rw_dbglevel, dmmoab->extra_read_options, usrreadopts, &readopts)); 872 873 PetscCall(PetscInfo(*dm, "Reading file %s with options: %s\n", filename, readopts)); 874 875 /* Load the mesh from a file. */ 876 if (dmmoab->fileset) { 877 merr = mbiface->load_file(filename, &dmmoab->fileset, readopts); 878 MBERRVM(mbiface, "Reading MOAB file failed.", merr); 879 } else { 880 merr = mbiface->load_file(filename, 0, readopts); 881 MBERRVM(mbiface, "Reading MOAB file failed.", merr); 882 } 883 884 #ifdef MOAB_HAVE_MPI 885 /* Reassign global IDs on all entities. */ 886 /* merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, true, true, true);MBERRNM(merr); */ 887 #endif 888 889 /* load the local vertices */ 890 merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true); 891 MBERRNM(merr); 892 /* load the local elements */ 893 merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true); 894 MBERRNM(merr); 895 896 #ifdef MOAB_HAVE_MPI 897 /* Everything is set up, now just do a tag exchange to update tags 898 on all of the ghost vertexes */ 899 merr = pcomm->exchange_tags(dmmoab->ltog_tag, verts); 900 MBERRV(mbiface, merr); 901 merr = pcomm->exchange_tags(dmmoab->ltog_tag, elems); 902 MBERRV(mbiface, merr); 903 merr = pcomm->collective_sync_partition(); 904 MBERR("Collective sync failed", merr); 905 #endif 906 907 PetscCall(PetscInfo(*dm, "MOAB file '%s' was successfully loaded. Found %zu vertices and %zu elements.\n", filename, verts.size(), elems.size())); 908 PetscCall(PetscFree(readopts)); 909 PetscFunctionReturn(0); 910 } 911 912 /*@C 913 DMMoabRenumberMeshEntities - Order and number all entities (vertices->elements) to be contiguously ordered 914 in parallel 915 916 Collective 917 918 Input Parameters: 919 . dm - The DM object 920 921 Level: advanced 922 923 .seealso: `DMSetUp()`, `DMCreate()` 924 @*/ 925 PetscErrorCode DMMoabRenumberMeshEntities(DM dm) 926 { 927 moab::Range verts; 928 929 PetscFunctionBegin; 930 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 931 932 #ifdef MOAB_HAVE_MPI 933 /* Insert new points */ 934 moab::ErrorCode merr; 935 merr = ((DM_Moab *)dm->data)->pcomm->assign_global_ids(((DM_Moab *)dm->data)->fileset, 3, 0, false, true, false); 936 MBERRNM(merr); 937 #endif 938 PetscFunctionReturn(0); 939 } 940