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