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 26 } DMMoabMeshGeneratorCtx; 27 28 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 } 62 63 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 } 184 185 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 } 206 207 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; 212 PetscErrorCode ierr; 213 PetscInt ii, jj, kk; 214 std::vector<PetscReal*> arrays; 215 PetscInt* gids; 216 moab::ErrorCode merr; 217 218 PetscFunctionBegin; 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 ierr = PetscMalloc1(nnodes, &gids);CHKERRQ(ierr); 225 226 merr = iface->get_node_coords(3, nnodes, 0, startv, arrays);MBERR("Can't get node coords.", merr); 227 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]; */ 232 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]); 247 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 ierr = PetscFree(gids);CHKERRQ(ierr); 260 PetscFunctionReturn(0); 261 } 262 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; 274 275 PetscFunctionBegin; 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 } 305 306 ix = ie = 0; /* index now in the elements, for global ids */ 307 308 /* create a temporary range to store local element handles */ 309 moab::Range tmp(starte, starte + nelems - 1); 310 std::vector<PetscInt> gids(nelems); 311 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); 316 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++) { 321 322 moab::EntityHandle corner = startv + genCtx.q * ii + genCtx.q * jj * genCtx.ystride + genCtx.q * kk * genCtx.zstride; 323 324 std::pair<PetscInt, PetscInt> entoffset = DMMoab_SetElementConnectivity_Private(genCtx, ix, corner, conn); 325 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 } 335 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 PetscFunctionReturn(0); 347 } 348 349 PetscErrorCode DMMBUtil_InitializeOptions(DMMoabMeshGeneratorCtx& genCtx, PetscInt dim, PetscBool simplex, PetscInt rank, PetscInt nprocs, const PetscReal* bounds, PetscInt nelems) 350 { 351 PetscFunctionBegin; 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; 359 360 if (!genCtx.usrxyzgrid) { /* not overridden by genCtx - assume nele equally and that genCtx wants a uniform cube mesh */ 361 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++; 367 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 } 388 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 } 399 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 } 416 417 if (!genCtx.usrrefgrid) { 418 genCtx.A = genCtx.B = genCtx.C = 1; 419 } 420 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; 427 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; 433 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; 451 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; 464 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 } 471 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]); 479 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 PetscFunctionReturn(0); 488 } 489 490 /*@C 491 DMMoabCreateBoxMesh - Creates a mesh on the tensor product (box) of intervals with genCtx specified bounds. 492 493 Collective 494 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 501 502 Output Parameter: 503 . dm - The DM object 504 505 Level: beginner 506 507 .seealso: DMSetType(), DMCreate(), DMMoabLoadFromFile() 508 @*/ 509 PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool useSimplex, const PetscReal* bounds, PetscInt nele, PetscInt nghost, DM *dm) 510 { 511 PetscErrorCode ierr; 512 moab::ErrorCode merr; 513 PetscInt a, b, c, n, global_size, global_rank; 514 DM_Moab *dmmoab; 515 moab::Interface *mbImpl; 516 #ifdef MOAB_HAVE_MPI 517 moab::ParallelComm *pcomm; 518 #endif 519 moab::ReadUtilIface *readMeshIface; 520 moab::Range verts, cells, edges, faces, adj, dim3, dim2; 521 DMMoabMeshGeneratorCtx genCtx; 522 const PetscInt npts = nele + 1; /* Number of points in every dimension */ 523 524 moab::Tag global_id_tag, part_tag, geom_tag, mat_tag, dir_tag, neu_tag; 525 moab::Range ownedvtx, ownedelms, localvtxs, localelms; 526 moab::EntityHandle regionset; 527 PetscInt ml = 0, nl = 0, kl = 0; 528 529 PetscFunctionBegin; 530 if (dim < 1 || dim > 3) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension argument for mesh: dim=[1,3].\n"); 531 532 ierr = PetscLogEventRegister("GenerateMesh", DM_CLASSID, &genCtx.generateMesh);CHKERRQ(ierr); 533 ierr = PetscLogEventRegister("AddVertices", DM_CLASSID, &genCtx.generateVertices);CHKERRQ(ierr); 534 ierr = PetscLogEventRegister("AddElements", DM_CLASSID, &genCtx.generateElements);CHKERRQ(ierr); 535 ierr = PetscLogEventRegister("ParResolve", DM_CLASSID, &genCtx.parResolve);CHKERRQ(ierr); 536 ierr = PetscLogEventBegin(genCtx.generateMesh, 0, 0, 0, 0);CHKERRQ(ierr); 537 ierr = MPI_Comm_size(comm, &global_size);CHKERRQ(ierr); 538 /* total number of vertices in all dimensions */ 539 n = pow(npts, dim); 540 541 /* do some error checking */ 542 if (n < 2) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be >= 2.\n"); 543 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"); 544 if (nghost < 0) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of ghost layers cannot be negative.\n"); 545 546 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 547 ierr = DMMoabCreateMoab(comm, NULL, NULL, NULL, dm);CHKERRQ(ierr); 548 549 /* get all the necessary handles from the private DM object */ 550 dmmoab = (DM_Moab*)(*dm)->data; 551 mbImpl = dmmoab->mbiface; 552 #ifdef MOAB_HAVE_MPI 553 pcomm = dmmoab->pcomm; 554 global_rank = pcomm->rank(); 555 #else 556 global_rank = 0; 557 global_size = 1; 558 #endif 559 global_id_tag = dmmoab->ltog_tag; 560 dmmoab->dim = dim; 561 dmmoab->nghostrings = nghost; 562 dmmoab->refct = 1; 563 564 /* create a file set to associate all entities in current mesh */ 565 merr = mbImpl->create_meshset(moab::MESHSET_SET, dmmoab->fileset);MBERR("Creating file set failed", merr); 566 567 /* No errors yet; proceed with building the mesh */ 568 merr = mbImpl->query_interface(readMeshIface);MBERRNM(merr); 569 570 genCtx.M = genCtx.N = genCtx.K = 1; 571 genCtx.A = genCtx.B = genCtx.C = 1; 572 genCtx.blockSizeElementXYZ[0] = 0; 573 genCtx.blockSizeElementXYZ[1] = 0; 574 genCtx.blockSizeElementXYZ[2] = 0; 575 576 ierr = PetscOptionsBegin(comm, "", "DMMoab Creation Options", "DMMOAB");CHKERRQ(ierr); 577 /* Handle DMMoab spatial resolution */ 578 ierr = PetscOptionsInt("-dmb_grid_x", "Number of grid points in x direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[0], &genCtx.blockSizeElementXYZ[0], &genCtx.usrxyzgrid);CHKERRQ(ierr); 579 if (dim > 1) {ierr = PetscOptionsInt("-dmb_grid_y", "Number of grid points in y direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[1], &genCtx.blockSizeElementXYZ[1], &genCtx.usrxyzgrid);CHKERRQ(ierr);} 580 if (dim > 2) {ierr = PetscOptionsInt("-dmb_grid_z", "Number of grid points in z direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[2], &genCtx.blockSizeElementXYZ[2], &genCtx.usrxyzgrid);CHKERRQ(ierr);} 581 582 /* Handle DMMoab parallel distibution */ 583 ierr = PetscOptionsInt("-dmb_processors_x", "Number of processors in x direction", "DMMoabSetNumProcs", genCtx.M, &genCtx.M, &genCtx.usrprocgrid);CHKERRQ(ierr); 584 if (dim > 1) {ierr = PetscOptionsInt("-dmb_processors_y", "Number of processors in y direction", "DMMoabSetNumProcs", genCtx.N, &genCtx.N, &genCtx.usrprocgrid);CHKERRQ(ierr);} 585 if (dim > 2) {ierr = PetscOptionsInt("-dmb_processors_z", "Number of processors in z direction", "DMMoabSetNumProcs", genCtx.K, &genCtx.K, &genCtx.usrprocgrid);CHKERRQ(ierr);} 586 587 /* Handle DMMoab block refinement */ 588 ierr = PetscOptionsInt("-dmb_refine_x", "Number of refinement blocks in x direction", "DMMoabSetRefinement", genCtx.A, &genCtx.A, &genCtx.usrrefgrid); 589 if (dim > 1) {ierr = PetscOptionsInt("-dmb_refine_y", "Number of refinement blocks in y direction", "DMMoabSetRefinement", genCtx.B, &genCtx.B, &genCtx.usrrefgrid);CHKERRQ(ierr);} 590 if (dim > 2) {ierr = PetscOptionsInt("-dmb_refine_z", "Number of refinement blocks in z direction", "DMMoabSetRefinement", genCtx.C, &genCtx.C, &genCtx.usrrefgrid);CHKERRQ(ierr);} 591 PetscOptionsEnd(); 592 593 ierr = DMMBUtil_InitializeOptions(genCtx, dim, useSimplex, global_rank, global_size, bounds, nele);CHKERRQ(ierr); 594 595 //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); 596 597 if (genCtx.adjEnts) genCtx.keep_skins = true; /* do not delete anything - consumes more memory */ 598 599 /* determine m, n, k for processor rank */ 600 ml = nl = kl = 0; 601 switch (genCtx.dim) { 602 case 1: 603 ml = (genCtx.cumfraction); 604 break; 605 case 2: 606 nl = (genCtx.cumfraction); 607 break; 608 default: 609 kl = (genCtx.cumfraction) / genCtx.q / genCtx.blockSizeElementXYZ[2] / genCtx.C; //genCtx.K 610 break; 611 } 612 613 /* 614 * so there are a total of M * A * blockSizeElement elements in x direction (so M * A * blockSizeElement + 1 verts in x direction) 615 * so there are a total of N * B * blockSizeElement elements in y direction (so N * B * blockSizeElement + 1 verts in y direction) 616 * so there are a total of K * C * blockSizeElement elements in z direction (so K * C * blockSizeElement + 1 verts in z direction) 617 618 * there are ( M * A blockSizeElement ) * ( N * B * blockSizeElement) * (K * C * blockSizeElement ) hexas 619 * there are ( M * A * blockSizeElement + 1) * ( N * B * blockSizeElement + 1 ) * (K * C * blockSizeElement + 1) vertices 620 * x is the first dimension that varies 621 */ 622 623 /* generate the block at (a, b, c); it will represent a partition , it will get a partition tag */ 624 PetscInt dum_id = -1; 625 merr = mbImpl->tag_get_handle("GLOBAL_ID", 1, moab::MB_TYPE_INTEGER, global_id_tag);MBERR("Getting Global_ID Tag handle failed", merr); 626 627 merr = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, mat_tag);MBERR("Getting Material set Tag handle failed", merr); 628 merr = mbImpl->tag_get_handle(DIRICHLET_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, dir_tag);MBERR("Getting Dirichlet set Tag handle failed", merr); 629 merr = mbImpl->tag_get_handle(NEUMANN_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, neu_tag);MBERR("Getting Neumann set Tag handle failed", merr); 630 631 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); 632 633 /* lets create some sets */ 634 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); 635 merr = mbImpl->create_meshset(moab::MESHSET_SET, regionset);MBERRNM(merr); 636 ierr = PetscLogEventEnd(genCtx.generateMesh, 0, 0, 0, 0);CHKERRQ(ierr); 637 638 for (a = 0; a < (genCtx.dim > 0 ? genCtx.A : genCtx.A); a++) { 639 for (b = 0; b < (genCtx.dim > 1 ? genCtx.B : 1); b++) { 640 for (c = 0; c < (genCtx.dim > 2 ? genCtx.C : 1); c++) { 641 642 moab::EntityHandle startv; 643 644 ierr = PetscLogEventBegin(genCtx.generateVertices, 0, 0, 0, 0);CHKERRQ(ierr); 645 ierr = DMMoab_GenerateVertices_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, verts);CHKERRQ(ierr); 646 ierr = PetscLogEventEnd(genCtx.generateVertices, 0, 0, 0, 0);CHKERRQ(ierr); 647 648 ierr = PetscLogEventBegin(genCtx.generateElements, 0, 0, 0, 0);CHKERRQ(ierr); 649 ierr = DMMoab_GenerateElements_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, cells);CHKERRQ(ierr); 650 ierr = PetscLogEventEnd(genCtx.generateElements, 0, 0, 0, 0);CHKERRQ(ierr); 651 652 PetscInt part_num = 0; 653 switch (genCtx.dim) { 654 case 3: 655 part_num += (c + kl * genCtx.C) * (genCtx.M * genCtx.A * genCtx.N * genCtx.B); 656 case 2: 657 part_num += (b + nl * genCtx.B) * (genCtx.M * genCtx.A); 658 case 1: 659 part_num += (a + ml * genCtx.A); 660 break; 661 } 662 663 moab::EntityHandle part_set; 664 merr = mbImpl->create_meshset(moab::MESHSET_SET, part_set);MBERR("Can't create mesh set.", merr); 665 666 merr = mbImpl->add_entities(part_set, verts);MBERR("Can't add vertices to set.", merr); 667 merr = mbImpl->add_entities(part_set, cells);MBERR("Can't add entities to set.", merr); 668 merr = mbImpl->add_entities(regionset, cells);MBERR("Can't add entities to set.", merr); 669 670 /* if needed, add all edges and faces */ 671 if (genCtx.adjEnts) 672 { 673 if (genCtx.dim > 1) { 674 merr = mbImpl->get_adjacencies(cells, 1, true, edges, moab::Interface::UNION);MBERR("Can't get edges", merr); 675 merr = mbImpl->add_entities(part_set, edges);MBERR("Can't add edges to partition set.", merr); 676 } 677 if (genCtx.dim > 2) { 678 merr = mbImpl->get_adjacencies(cells, 2, true, faces, moab::Interface::UNION);MBERR("Can't get faces", merr); 679 merr = mbImpl->add_entities(part_set, faces);MBERR("Can't add faces to partition set.", merr); 680 } 681 edges.clear(); 682 faces.clear(); 683 } 684 verts.clear(); cells.clear(); 685 686 merr = mbImpl->tag_set_data(part_tag, &part_set, 1, &part_num);MBERR("Can't set part tag on set", merr); 687 if (dmmoab->fileset) { 688 merr = mbImpl->add_parent_child(dmmoab->fileset, part_set);MBERR("Can't add part set to file set.", merr); 689 merr = mbImpl->unite_meshset(dmmoab->fileset, part_set);MBERRNM(merr); 690 } 691 merr = mbImpl->add_entities(dmmoab->fileset, &part_set, 1);MBERRNM(merr); 692 } 693 } 694 } 695 696 merr = mbImpl->add_parent_child(dmmoab->fileset, regionset);MBERRNM(merr); 697 698 /* Only in parallel: resolve shared entities between processors and exchange ghost layers */ 699 if (global_size > 1) { 700 701 ierr = PetscLogEventBegin(genCtx.parResolve, 0, 0, 0, 0);CHKERRQ(ierr); 702 703 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, genCtx.dim, cells);MBERR("Can't get all d-dimensional elements.", merr); 704 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 0, verts);MBERR("Can't get all vertices.", merr); 705 706 if (genCtx.A * genCtx.B * genCtx.C != 1) { // merge needed 707 moab::MergeMesh mm(mbImpl); 708 if (genCtx.newMergeMethod) { 709 merr = mm.merge_using_integer_tag(verts, global_id_tag);MBERR("Can't merge with GLOBAL_ID tag", merr); 710 } 711 else { 712 merr = mm.merge_entities(cells, 0.0001);MBERR("Can't merge with coordinates", merr); 713 } 714 } 715 716 #ifdef MOAB_HAVE_MPI 717 /* check the handles */ 718 merr = pcomm->check_all_shared_handles();MBERRV(mbImpl, merr); 719 720 /* resolve the shared entities by exchanging information to adjacent processors */ 721 merr = pcomm->resolve_shared_ents(dmmoab->fileset, cells, dim, dim - 1, NULL, &global_id_tag);MBERRV(mbImpl, merr); 722 if (dmmoab->fileset) { 723 merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false, &dmmoab->fileset);MBERRV(mbImpl, merr); 724 } 725 else { 726 merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false);MBERRV(mbImpl, merr); 727 } 728 729 /* Reassign global IDs on all entities. */ 730 merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, false, true, false);MBERRNM(merr); 731 #endif 732 733 ierr = PetscLogEventEnd(genCtx.parResolve, 0, 0, 0, 0);CHKERRQ(ierr); 734 } 735 736 if (!genCtx.keep_skins) { // default is to delete the 1- and 2-dimensional entities 737 // delete all quads and edges 738 moab::Range toDelete; 739 if (genCtx.dim > 1) { 740 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 1, toDelete);MBERR("Can't get edges", merr); 741 } 742 743 if (genCtx.dim > 2) { 744 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 2, toDelete);MBERR("Can't get faces", merr); 745 } 746 747 #ifdef MOAB_HAVE_MPI 748 merr = dmmoab->pcomm->delete_entities(toDelete) ;MBERR("Can't delete entities", merr); 749 #endif 750 } 751 752 /* set geometric dimension tag for regions */ 753 merr = mbImpl->tag_set_data(geom_tag, ®ionset, 1, &dmmoab->dim);MBERRNM(merr); 754 /* set default material ID for regions */ 755 int default_material = 1; 756 merr = mbImpl->tag_set_data(mat_tag, ®ionset, 1, &default_material);MBERRNM(merr); 757 /* 758 int default_dbc = 0; 759 merr = mbImpl->tag_set_data(dir_tag, &vertexset, 1, &default_dbc);MBERRNM(merr); 760 */ 761 PetscFunctionReturn(0); 762 } 763 764 765 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) 766 { 767 char *ropts; 768 char ropts_par[PETSC_MAX_PATH_LEN], ropts_pargh[PETSC_MAX_PATH_LEN]; 769 char ropts_dbg[PETSC_MAX_PATH_LEN]; 770 PetscErrorCode ierr; 771 772 PetscFunctionBegin; 773 ierr = PetscMalloc1(PETSC_MAX_PATH_LEN, &ropts);CHKERRQ(ierr); 774 ierr = PetscMemzero(&ropts_par, PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 775 ierr = PetscMemzero(&ropts_pargh, PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 776 ierr = PetscMemzero(&ropts_dbg, PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 777 778 /* do parallel read unless using only one processor */ 779 if (numproc > 1) { 780 // ierr = 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":""));CHKERRQ(ierr); 781 ierr = 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;" : ""));CHKERRQ(ierr); 782 if (nghost) { 783 ierr = PetscSNPrintf(ropts_pargh, PETSC_MAX_PATH_LEN, "PARALLEL_GHOSTS=%d.0.%d;", dim, nghost);CHKERRQ(ierr); 784 } 785 } 786 787 if (dbglevel) { 788 if (numproc > 1) { 789 ierr = PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%d;DEBUG_PIO=%d;", dbglevel, dbglevel);CHKERRQ(ierr); 790 } 791 else { 792 ierr = PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%d;", dbglevel);CHKERRQ(ierr); 793 } 794 } 795 796 ierr = 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 : ""));CHKERRQ(ierr); 797 *read_opts = ropts; 798 PetscFunctionReturn(0); 799 } 800 801 802 /*@C 803 DMMoabLoadFromFile - Creates a DM object by loading the mesh from a user specified file. 804 805 Collective 806 807 Input Parameters: 808 + comm - The communicator for the DM object 809 . dim - The spatial dimension 810 . filename - The name of the mesh file to be loaded 811 - usrreadopts - The options string to read a MOAB mesh. 812 813 Reference (Parallel Mesh Initialization: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html#fivetwo) 814 815 Output Parameter: 816 . dm - The DM object 817 818 Level: beginner 819 820 .seealso: DMSetType(), DMCreate(), DMMoabCreateBoxMesh() 821 @*/ 822 PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm, PetscInt dim, PetscInt nghost, const char* filename, const char* usrreadopts, DM *dm) 823 { 824 moab::ErrorCode merr; 825 PetscInt nprocs; 826 DM_Moab *dmmoab; 827 moab::Interface *mbiface; 828 #ifdef MOAB_HAVE_MPI 829 moab::ParallelComm *pcomm; 830 #endif 831 moab::Range verts, elems; 832 const char *readopts; 833 PetscErrorCode ierr; 834 835 PetscFunctionBegin; 836 PetscValidPointer(dm, 6); 837 838 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 839 ierr = DMMoabCreateMoab(comm, NULL, NULL, NULL, dm);CHKERRQ(ierr); 840 841 /* get all the necessary handles from the private DM object */ 842 dmmoab = (DM_Moab*)(*dm)->data; 843 mbiface = dmmoab->mbiface; 844 #ifdef MOAB_HAVE_MPI 845 pcomm = dmmoab->pcomm; 846 nprocs = pcomm->size(); 847 #else 848 nprocs = 1; 849 #endif 850 /* TODO: Decipher dimension based on the loaded mesh instead of getting from user */ 851 dmmoab->dim = dim; 852 dmmoab->nghostrings = nghost; 853 dmmoab->refct = 1; 854 855 /* create a file set to associate all entities in current mesh */ 856 merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, dmmoab->fileset);MBERR("Creating file set failed", merr); 857 858 /* add mesh loading options specific to the DM */ 859 ierr = DMMoab_GetReadOptions_Private(dmmoab->partition_by_rank, nprocs, dim, nghost, dmmoab->read_mode, 860 dmmoab->rw_dbglevel, dmmoab->extra_read_options, usrreadopts, &readopts);CHKERRQ(ierr); 861 862 PetscInfo2(*dm, "Reading file %s with options: %s\n", filename, readopts); 863 864 /* Load the mesh from a file. */ 865 if (dmmoab->fileset) { 866 merr = mbiface->load_file(filename, &dmmoab->fileset, readopts);MBERRVM(mbiface, "Reading MOAB file failed.", merr); 867 } 868 else { 869 merr = mbiface->load_file(filename, 0, readopts);MBERRVM(mbiface, "Reading MOAB file failed.", merr); 870 } 871 872 #ifdef MOAB_HAVE_MPI 873 /* Reassign global IDs on all entities. */ 874 /* merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, true, true, true);MBERRNM(merr); */ 875 #endif 876 877 /* load the local vertices */ 878 merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);MBERRNM(merr); 879 /* load the local elements */ 880 merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);MBERRNM(merr); 881 882 #ifdef MOAB_HAVE_MPI 883 /* Everything is set up, now just do a tag exchange to update tags 884 on all of the ghost vertexes */ 885 merr = pcomm->exchange_tags(dmmoab->ltog_tag, verts);MBERRV(mbiface, merr); 886 merr = pcomm->exchange_tags(dmmoab->ltog_tag, elems);MBERRV(mbiface, merr); 887 merr = pcomm->collective_sync_partition();MBERR("Collective sync failed", merr); 888 #endif 889 890 PetscInfo3(*dm, "MOAB file '%s' was successfully loaded. Found %D vertices and %D elements.\n", filename, verts.size(), elems.size()); 891 ierr = PetscFree(readopts);CHKERRQ(ierr); 892 PetscFunctionReturn(0); 893 } 894 895 896 /*@C 897 DMMoabRenumberMeshEntities - Order and number all entities (vertices->elements) to be contiguously ordered 898 in parallel 899 900 Collective 901 902 Input Parameters: 903 . dm - The DM object 904 905 Level: advanced 906 907 .seealso: DMSetUp(), DMCreate() 908 @*/ 909 PetscErrorCode DMMoabRenumberMeshEntities(DM dm) 910 { 911 moab::Range verts; 912 913 PetscFunctionBegin; 914 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 915 916 #ifdef MOAB_HAVE_MPI 917 /* Insert new points */ 918 moab::ErrorCode merr; 919 merr = ((DM_Moab*) dm->data)->pcomm->assign_global_ids(((DM_Moab*) dm->data)->fileset, 3, 0, false, true, false);MBERRNM(merr); 920 #endif 921 PetscFunctionReturn(0); 922 } 923 924