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 PetscInt DMMoab_SetTensorElementConnectivity_Private(DMMoabMeshGeneratorCtx& genCtx, PetscInt offset, PetscInt corner, std::vector<PetscInt>& subent_conn, moab::EntityHandle *connectivity) 29 { 30 switch (genCtx.dim) { 31 case 1: 32 subent_conn.resize(2); 33 moab::CN::SubEntityVertexIndices(moab::MBEDGE, 1, 0, subent_conn.data()); 34 connectivity[offset + subent_conn[0]] = corner; 35 connectivity[offset + subent_conn[1]] = corner + 1; 36 break; 37 case 2: 38 subent_conn.resize(4); 39 moab::CN::SubEntityVertexIndices(moab::MBQUAD, 2, 0, subent_conn.data()); 40 connectivity[offset + subent_conn[0]] = corner; 41 connectivity[offset + subent_conn[1]] = corner + 1; 42 connectivity[offset + subent_conn[2]] = corner + 1 + genCtx.ystride; 43 connectivity[offset + subent_conn[3]] = corner + genCtx.ystride; 44 break; 45 case 3: 46 default: 47 subent_conn.resize(8); 48 moab::CN::SubEntityVertexIndices(moab::MBHEX, 3, 0, subent_conn.data()); 49 connectivity[offset + subent_conn[0]] = corner; 50 connectivity[offset + subent_conn[1]] = corner + 1; 51 connectivity[offset + subent_conn[2]] = corner + 1 + genCtx.ystride; 52 connectivity[offset + subent_conn[3]] = corner + genCtx.ystride; 53 connectivity[offset + subent_conn[4]] = corner + genCtx.zstride; 54 connectivity[offset + subent_conn[5]] = corner + 1 + genCtx.zstride; 55 connectivity[offset + subent_conn[6]] = corner + 1 + genCtx.ystride + genCtx.zstride; 56 connectivity[offset + subent_conn[7]] = corner + genCtx.ystride + genCtx.zstride; 57 break; 58 } 59 return subent_conn.size(); 60 } 61 62 PetscInt DMMoab_SetSimplexElementConnectivity_Private(DMMoabMeshGeneratorCtx& genCtx, PetscInt subelem, PetscInt offset, PetscInt corner, std::vector<PetscInt>& subent_conn, moab::EntityHandle *connectivity) 63 { 64 PetscInt A, B, C, D, E, F, G, H, M; 65 const PetscInt trigen_opts = 1; /* 1 - Aligned diagonally to right, 2 - Aligned diagonally to left, 3 - 4 elements per quad */ 66 A = corner; 67 B = corner + 1; 68 switch (genCtx.dim) { 69 case 1: 70 subent_conn.resize(2); /* only linear EDGE supported now */ 71 moab::CN::SubEntityVertexIndices(moab::MBEDGE, 1, 0, subent_conn.data()); 72 connectivity[offset + subent_conn[0]] = A; 73 connectivity[offset + subent_conn[1]] = B; 74 break; 75 case 2: 76 C = corner + 1 + genCtx.ystride; 77 D = corner + genCtx.ystride; 78 M = corner + 0.5; /* technically -- need to modify vertex generation */ 79 subent_conn.resize(3); /* only linear TRI supported */ 80 moab::CN::SubEntityVertexIndices(moab::MBTRI, 2, 0, subent_conn.data()); 81 if (trigen_opts == 1) { 82 if (subelem) { /* 0 1 2 of a QUAD */ 83 connectivity[offset + subent_conn[0]] = B; 84 connectivity[offset + subent_conn[1]] = C; 85 connectivity[offset + subent_conn[2]] = A; 86 } 87 else { /* 2 3 0 of a QUAD */ 88 connectivity[offset + subent_conn[0]] = D; 89 connectivity[offset + subent_conn[1]] = A; 90 connectivity[offset + subent_conn[2]] = C; 91 } 92 } 93 else if (trigen_opts == 2) { 94 if (subelem) { /* 0 1 2 of a QUAD */ 95 connectivity[offset + subent_conn[0]] = A; 96 connectivity[offset + subent_conn[1]] = B; 97 connectivity[offset + subent_conn[2]] = D; 98 } 99 else { /* 2 3 0 of a QUAD */ 100 connectivity[offset + subent_conn[0]] = C; 101 connectivity[offset + subent_conn[1]] = D; 102 connectivity[offset + subent_conn[2]] = B; 103 } 104 } 105 else { 106 switch (subelem) { /* 0 1 2 of a QUAD */ 107 case 0: 108 connectivity[offset + subent_conn[0]] = A; 109 connectivity[offset + subent_conn[1]] = B; 110 connectivity[offset + subent_conn[2]] = M; 111 break; 112 case 1: 113 connectivity[offset + subent_conn[0]] = B; 114 connectivity[offset + subent_conn[1]] = C; 115 connectivity[offset + subent_conn[2]] = M; 116 break; 117 case 2: 118 connectivity[offset + subent_conn[0]] = C; 119 connectivity[offset + subent_conn[1]] = D; 120 connectivity[offset + subent_conn[2]] = M; 121 break; 122 case 3: 123 connectivity[offset + subent_conn[0]] = D; 124 connectivity[offset + subent_conn[1]] = A; 125 connectivity[offset + subent_conn[2]] = M; 126 break; 127 } 128 } 129 break; 130 case 3: 131 default: 132 C = corner + 1 + genCtx.ystride; 133 D = corner + genCtx.ystride; 134 E = corner + genCtx.zstride; 135 F = corner + 1 + genCtx.zstride; 136 G = corner + 1 + genCtx.ystride + genCtx.zstride; 137 H = corner + genCtx.ystride + genCtx.zstride; 138 subent_conn.resize(4); /* only linear TET supported */ 139 moab::CN::SubEntityVertexIndices(moab::MBTET, 3, 0, subent_conn.data()); 140 switch (subelem) { 141 case 0: /* 4 3 7 6 of a HEX */ 142 connectivity[offset + subent_conn[0]] = E; 143 connectivity[offset + subent_conn[1]] = D; 144 connectivity[offset + subent_conn[2]] = H; 145 connectivity[offset + subent_conn[3]] = G; 146 break; 147 case 1: /* 0 1 2 5 of a HEX */ 148 connectivity[offset + subent_conn[0]] = A; 149 connectivity[offset + subent_conn[1]] = B; 150 connectivity[offset + subent_conn[2]] = C; 151 connectivity[offset + subent_conn[3]] = F; 152 break; 153 case 2: /* 0 3 4 5 of a HEX */ 154 connectivity[offset + subent_conn[0]] = A; 155 connectivity[offset + subent_conn[1]] = D; 156 connectivity[offset + subent_conn[2]] = E; 157 connectivity[offset + subent_conn[3]] = F; 158 break; 159 case 3: /* 2 6 3 5 of a HEX */ 160 connectivity[offset + subent_conn[0]] = C; 161 connectivity[offset + subent_conn[1]] = G; 162 connectivity[offset + subent_conn[2]] = D; 163 connectivity[offset + subent_conn[3]] = F; 164 break; 165 case 4: /* 0 2 3 5 of a HEX */ 166 connectivity[offset + subent_conn[0]] = A; 167 connectivity[offset + subent_conn[1]] = C; 168 connectivity[offset + subent_conn[2]] = D; 169 connectivity[offset + subent_conn[3]] = F; 170 break; 171 case 5: /* 3 6 4 5 of a HEX */ 172 connectivity[offset + subent_conn[0]] = D; 173 connectivity[offset + subent_conn[1]] = G; 174 connectivity[offset + subent_conn[2]] = E; 175 connectivity[offset + subent_conn[3]] = F; 176 break; 177 } 178 break; 179 } 180 return subent_conn.size(); 181 } 182 183 std::pair<PetscInt, PetscInt> DMMoab_SetElementConnectivity_Private(DMMoabMeshGeneratorCtx& genCtx, PetscInt offset, PetscInt corner, moab::EntityHandle *connectivity) 184 { 185 PetscInt vcount = 0; 186 PetscInt simplices_per_tensor[4] = {0, 1, 2, 6}; 187 std::vector<PetscInt> subent_conn; /* only linear edge, tri, tet supported now */ 188 subent_conn.reserve(27); 189 PetscInt m, subelem; 190 if (genCtx.simplex) { 191 subelem = simplices_per_tensor[genCtx.dim]; 192 for (m = 0; m < subelem; m++) { 193 vcount = DMMoab_SetSimplexElementConnectivity_Private(genCtx, m, offset, corner, subent_conn, connectivity); 194 offset += vcount; 195 } 196 } 197 else { 198 subelem = 1; 199 vcount = DMMoab_SetTensorElementConnectivity_Private(genCtx, offset, corner, subent_conn, connectivity); 200 } 201 return std::pair<PetscInt, PetscInt>(vcount * subelem, subelem); 202 } 203 204 PetscErrorCode DMMoab_GenerateVertices_Private(moab::Interface *mbImpl, moab::ReadUtilIface *iface, DMMoabMeshGeneratorCtx& genCtx, PetscInt m, PetscInt n, PetscInt k, 205 PetscInt a, PetscInt b, PetscInt c, moab::Tag& global_id_tag, moab::EntityHandle& startv, moab::Range& uverts) 206 { 207 PetscInt x, y, z, ix, nnodes; 208 PetscErrorCode ierr; 209 PetscInt ii, jj, kk; 210 std::vector<PetscReal*> arrays; 211 PetscInt* gids; 212 moab::ErrorCode merr; 213 214 PetscFunctionBegin; 215 /* we will generate (q*block+1)^3 vertices, and block^3 hexas; q is 1 for linear, 2 for quadratic 216 * the global id of the vertices will come from m, n, k, a, b, c 217 * x will vary from m*A*q*block + a*q*block to m*A*q*block+(a+1)*q*block etc. 218 */ 219 nnodes = genCtx.blockSizeVertexXYZ[0] * (genCtx.dim > 1 ? genCtx.blockSizeVertexXYZ[1] * (genCtx.dim > 2 ? genCtx.blockSizeVertexXYZ[2] : 1) : 1); 220 ierr = PetscMalloc1(nnodes, &gids);CHKERRQ(ierr); 221 222 merr = iface->get_node_coords(3, nnodes, 0, startv, arrays);MBERR("Can't get node coords.", merr); 223 224 /* will start with the lower corner: */ 225 /* x = ( m * genCtx.A + a) * genCtx.q * genCtx.blockSizeElementXYZ[0]; */ 226 /* y = ( n * genCtx.B + b) * genCtx.q * genCtx.blockSizeElementXYZ[1]; */ 227 /* z = ( k * genCtx.C + c) * genCtx.q * genCtx.blockSizeElementXYZ[2]; */ 228 229 x = ( m * genCtx.A + a) * genCtx.q; 230 y = ( n * genCtx.B + b) * genCtx.q; 231 z = ( k * genCtx.C + c) * genCtx.q; 232 PetscInfo3(NULL, "Starting offset for coordinates := %d, %d, %d\n", x, y, z); 233 ix = 0; 234 moab::Range verts(startv, startv + nnodes - 1); 235 for (kk = 0; kk < (genCtx.dim > 2 ? genCtx.blockSizeVertexXYZ[2] : 1); kk++) { 236 for (jj = 0; jj < (genCtx.dim > 1 ? genCtx.blockSizeVertexXYZ[1] : 1); jj++) { 237 for (ii = 0; ii < genCtx.blockSizeVertexXYZ[0]; ii++, ix++) { 238 /* set coordinates for the vertices */ 239 arrays[0][ix] = (x + ii) * genCtx.dx + genCtx.xyzbounds[0]; 240 arrays[1][ix] = (y + jj) * genCtx.dy + genCtx.xyzbounds[2]; 241 arrays[2][ix] = (z + kk) * genCtx.dz + genCtx.xyzbounds[4]; 242 PetscInfo3(NULL, "Creating vertex with coordinates := %f, %f, %f\n", arrays[0][ix], arrays[1][ix], arrays[2][ix]); 243 244 /* If we want to set some tags on the vertices -> use the following entity handle definition: 245 moab::EntityHandle v = startv + ix; 246 */ 247 /* compute the global ID for vertex */ 248 gids[ix] = 1 + (x + ii) + (y + jj) * genCtx.NX + (z + kk) * (genCtx.NX * genCtx.NY); 249 } 250 } 251 } 252 /* set global ID data on vertices */ 253 mbImpl->tag_set_data(global_id_tag, verts, &gids[0]); 254 verts.swap(uverts); 255 ierr = PetscFree(gids);CHKERRQ(ierr); 256 PetscFunctionReturn(0); 257 } 258 259 PetscErrorCode DMMoab_GenerateElements_Private(moab::Interface* mbImpl, moab::ReadUtilIface* iface, DMMoabMeshGeneratorCtx& genCtx, PetscInt m, PetscInt n, PetscInt k, 260 PetscInt a, PetscInt b, PetscInt c, moab::Tag& global_id_tag, moab::EntityHandle startv, moab::Range& cells) 261 { 262 moab::ErrorCode merr; 263 PetscInt ix, ie, xe, ye, ze; 264 PetscInt ii, jj, kk, nvperelem; 265 PetscInt simplices_per_tensor[4] = {0, 1, 2, 6}; 266 PetscInt ntensorelems = genCtx.blockSizeElementXYZ[0] * (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] * (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1) : 1); /*pow(genCtx.blockSizeElement,genCtx.dim);*/ 267 PetscInt nelems = ntensorelems; 268 moab::EntityHandle starte; /* connectivity */ 269 moab::EntityHandle* conn; 270 271 PetscFunctionBegin; 272 switch (genCtx.dim) { 273 case 1: 274 nvperelem = 2; 275 merr = iface->get_element_connect(nelems, 2, moab::MBEDGE, 0, starte, conn);MBERR("Can't get EDGE2 element connectivity.", merr); 276 break; 277 case 2: 278 if (genCtx.simplex) { 279 nvperelem = 3; 280 nelems = ntensorelems * simplices_per_tensor[genCtx.dim]; 281 merr = iface->get_element_connect(nelems, 3, moab::MBTRI, 0, starte, conn);MBERR("Can't get TRI3 element connectivity.", merr); 282 } 283 else { 284 nvperelem = 4; 285 merr = iface->get_element_connect(nelems, 4, moab::MBQUAD, 0, starte, conn);MBERR("Can't get QUAD4 element connectivity.", merr); 286 } 287 break; 288 case 3: 289 default: 290 if (genCtx.simplex) { 291 nvperelem = 4; 292 nelems = ntensorelems * simplices_per_tensor[genCtx.dim]; 293 merr = iface->get_element_connect(nelems, 4, moab::MBTET, 0, starte, conn);MBERR("Can't get TET4 element connectivity.", merr); 294 } 295 else { 296 nvperelem = 8; 297 merr = iface->get_element_connect(nelems, 8, moab::MBHEX, 0, starte, conn);MBERR("Can't get HEX8 element connectivity.", merr); 298 } 299 break; 300 } 301 302 ix = ie = 0; /* index now in the elements, for global ids */ 303 304 /* create a temporary range to store local element handles */ 305 moab::Range tmp(starte, starte + nelems - 1); 306 std::vector<PetscInt> gids(nelems); 307 308 /* identify the elements at the lower corner, for their global ids */ 309 xe = m * genCtx.A * genCtx.blockSizeElementXYZ[0] + a * genCtx.blockSizeElementXYZ[0]; 310 ye = (genCtx.dim > 1 ? n * genCtx.B * genCtx.blockSizeElementXYZ[1] + b * genCtx.blockSizeElementXYZ[1] : 0); 311 ze = (genCtx.dim > 2 ? k * genCtx.C * genCtx.blockSizeElementXYZ[2] + c * genCtx.blockSizeElementXYZ[2] : 0); 312 313 /* create owned elements requested by genCtx */ 314 for (kk = 0; kk < (genCtx.dim > 2 ? genCtx.blockSizeElementXYZ[2] : 1); kk++) { 315 for (jj = 0; jj < (genCtx.dim > 1 ? genCtx.blockSizeElementXYZ[1] : 1); jj++) { 316 for (ii = 0; ii < genCtx.blockSizeElementXYZ[0]; ii++) { 317 318 moab::EntityHandle corner = startv + genCtx.q * ii + genCtx.q * jj * genCtx.ystride + genCtx.q * kk * genCtx.zstride; 319 320 std::pair<PetscInt, PetscInt> entoffset = DMMoab_SetElementConnectivity_Private(genCtx, ix, corner, conn); 321 322 for (PetscInt j = 0; j < entoffset.second; j++) { 323 /* The entity handle for the particular element -> if we want to set some tags is 324 moab::EntityHandle eh = starte + ie + j; 325 */ 326 gids[ie + j] = 1 + ((xe + ii) + (ye + jj) * genCtx.nex + (ze + kk) * (genCtx.nex * genCtx.ney)); 327 /* gids[ie+j] = ie + j + ((xe + ii) + (ye + jj) * genCtx.nex + (ze + kk) * (genCtx.nex * genCtx.ney)); */ 328 /* gids[ie+j] = 1 + ie; */ 329 /* ie++; */ 330 } 331 332 ix += entoffset.first; 333 ie += entoffset.second; 334 } 335 } 336 } 337 if (genCtx.adjEnts) { /* we need to update adjacencies now, because some elements are new */ 338 merr = iface->update_adjacencies(starte, nelems, nvperelem, conn);MBERR("Can't update adjacencies", merr); 339 } 340 tmp.swap(cells); 341 merr = mbImpl->tag_set_data(global_id_tag, cells, &gids[0]);MBERR("Can't set global ids to elements.", merr); 342 PetscFunctionReturn(0); 343 } 344 345 PetscErrorCode DMMBUtil_InitializeOptions(DMMoabMeshGeneratorCtx& genCtx, PetscInt dim, PetscBool simplex, PetscInt rank, PetscInt nprocs, const PetscReal* bounds, PetscInt nelems) 346 { 347 PetscFunctionBegin; 348 /* Initialize all genCtx data */ 349 genCtx.dim = dim; 350 genCtx.simplex = simplex; 351 genCtx.newMergeMethod = genCtx.keep_skins = genCtx.adjEnts = true; 352 /* determine other global quantities for the mesh used for nodes increments */ 353 genCtx.q = 1; 354 genCtx.fraction = genCtx.remainder = genCtx.cumfraction = 0; 355 356 if (!genCtx.usrxyzgrid) { /* not overridden by genCtx - assume nele equally and that genCtx wants a uniform cube mesh */ 357 358 genCtx.fraction = nelems / nprocs; /* partition only by the largest dimension */ 359 genCtx.remainder = nelems % nprocs; /* remainder after partition which gets evenly distributed by round-robin */ 360 genCtx.cumfraction = (rank > 0 ? (genCtx.fraction) * (rank) + (rank - 1 < genCtx.remainder ? rank : genCtx.remainder) : 0); 361 if (rank < genCtx.remainder) /* This process gets "fraction+1" elements */ 362 genCtx.fraction++; 363 364 PetscInfo3(NULL, "Fraction = %D, Remainder = %D, Cumulative fraction = %D\n", genCtx.fraction, genCtx.remainder, genCtx.cumfraction); 365 switch (genCtx.dim) { 366 case 1: 367 genCtx.blockSizeElementXYZ[0] = genCtx.fraction; 368 genCtx.blockSizeElementXYZ[1] = 1; 369 genCtx.blockSizeElementXYZ[2] = 1; 370 break; 371 case 2: 372 genCtx.blockSizeElementXYZ[0] = nelems; 373 genCtx.blockSizeElementXYZ[1] = genCtx.fraction; 374 genCtx.blockSizeElementXYZ[2] = 1; 375 break; 376 case 3: 377 default: 378 genCtx.blockSizeElementXYZ[0] = nelems; 379 genCtx.blockSizeElementXYZ[1] = nelems; 380 genCtx.blockSizeElementXYZ[2] = genCtx.fraction; 381 break; 382 } 383 } 384 385 /* partition only by the largest dimension */ 386 /* Total number of local elements := genCtx.blockSizeElementXYZ[0]*(genCtx.dim>1? genCtx.blockSizeElementXYZ[1]*(genCtx.dim>2 ? genCtx.blockSizeElementXYZ[2]:1) :1); */ 387 if (bounds) { 388 for (PetscInt i = 0; i < 6; i++) 389 genCtx.xyzbounds[i] = bounds[i]; 390 } 391 else { 392 genCtx.xyzbounds[0] = genCtx.xyzbounds[2] = genCtx.xyzbounds[4] = 0.0; 393 genCtx.xyzbounds[1] = genCtx.xyzbounds[3] = genCtx.xyzbounds[5] = 1.0; 394 } 395 396 if (!genCtx.usrprocgrid) { 397 switch (genCtx.dim) { 398 case 1: 399 genCtx.M = nprocs; 400 genCtx.N = genCtx.K = 1; 401 break; 402 case 2: 403 genCtx.N = nprocs; 404 genCtx.M = genCtx.K = 1; 405 break; 406 default: 407 genCtx.K = nprocs; 408 genCtx.M = genCtx.N = 1; 409 break; 410 } 411 } 412 413 if (!genCtx.usrrefgrid) { 414 genCtx.A = genCtx.B = genCtx.C = 1; 415 } 416 417 /* more default values */ 418 genCtx.nex = genCtx.ney = genCtx.nez = 0; 419 genCtx.xstride = genCtx.ystride = genCtx.zstride = 0; 420 genCtx.NX = genCtx.NY = genCtx.NZ = 0; 421 genCtx.nex = genCtx.ney = genCtx.nez = 0; 422 genCtx.blockSizeVertexXYZ[0] = genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 1; 423 424 switch (genCtx.dim) { 425 case 3: 426 genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1; 427 genCtx.blockSizeVertexXYZ[1] = genCtx.q * genCtx.blockSizeElementXYZ[1] + 1; 428 genCtx.blockSizeVertexXYZ[2] = genCtx.q * genCtx.blockSizeElementXYZ[2] + 1; 429 430 genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0]; /* number of elements in x direction, used for global id on element */ 431 genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (nelems * genCtx.q); /* distance between 2 nodes in x direction */ 432 genCtx.NX = (genCtx.q * genCtx.nex + 1); 433 genCtx.xstride = 1; 434 genCtx.ney = genCtx.N * genCtx.B * genCtx.blockSizeElementXYZ[1]; /* number of elements in y direction .... */ 435 genCtx.dy = (genCtx.xyzbounds[3] - genCtx.xyzbounds[2]) / (nelems * genCtx.q); /* distance between 2 nodes in y direction */ 436 genCtx.NY = (genCtx.q * genCtx.ney + 1); 437 genCtx.ystride = genCtx.blockSizeVertexXYZ[0]; 438 genCtx.nez = genCtx.K * genCtx.C * genCtx.blockSizeElementXYZ[2]; /* number of elements in z direction .... */ 439 genCtx.dz = (genCtx.xyzbounds[5] - genCtx.xyzbounds[4]) / (nelems * genCtx.q); /* distance between 2 nodes in z direction */ 440 genCtx.NZ = (genCtx.q * genCtx.nez + 1); 441 genCtx.zstride = genCtx.blockSizeVertexXYZ[0] * genCtx.blockSizeVertexXYZ[1]; 442 break; 443 case 2: 444 genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1; 445 genCtx.blockSizeVertexXYZ[1] = genCtx.q * genCtx.blockSizeElementXYZ[1] + 1; 446 genCtx.blockSizeVertexXYZ[2] = 0; 447 448 genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0]; /* number of elements in x direction, used for global id on element */ 449 genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (genCtx.nex * genCtx.q); /* distance between 2 nodes in x direction */ 450 genCtx.NX = (genCtx.q * genCtx.nex + 1); 451 genCtx.xstride = 1; 452 genCtx.ney = genCtx.N * genCtx.B * genCtx.blockSizeElementXYZ[1]; /* number of elements in y direction .... */ 453 genCtx.dy = (genCtx.xyzbounds[3] - genCtx.xyzbounds[2]) / (nelems * genCtx.q); /* distance between 2 nodes in y direction */ 454 genCtx.NY = (genCtx.q * genCtx.ney + 1); 455 genCtx.ystride = genCtx.blockSizeVertexXYZ[0]; 456 break; 457 case 1: 458 genCtx.blockSizeVertexXYZ[1] = genCtx.blockSizeVertexXYZ[2] = 0; 459 genCtx.blockSizeVertexXYZ[0] = genCtx.q * genCtx.blockSizeElementXYZ[0] + 1; 460 461 genCtx.nex = genCtx.M * genCtx.A * genCtx.blockSizeElementXYZ[0]; /* number of elements in x direction, used for global id on element */ 462 genCtx.dx = (genCtx.xyzbounds[1] - genCtx.xyzbounds[0]) / (nelems * genCtx.q); /* distance between 2 nodes in x direction */ 463 genCtx.NX = (genCtx.q * genCtx.nex + 1); 464 genCtx.xstride = 1; 465 break; 466 } 467 468 /* Lets check for some valid input */ 469 if (genCtx.dim < 1 || genCtx.dim > 3) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid topological dimension specified: %d.\n", genCtx.dim); 470 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); 471 /* validate the bounds data */ 472 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]); 473 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]); 474 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]); 475 476 PetscInfo3(NULL, "Local elements:= %d, %d, %d\n", genCtx.blockSizeElementXYZ[0], genCtx.blockSizeElementXYZ[1], genCtx.blockSizeElementXYZ[2]); 477 PetscInfo3(NULL, "Local vertices:= %d, %d, %d\n", genCtx.blockSizeVertexXYZ[0], genCtx.blockSizeVertexXYZ[1], genCtx.blockSizeVertexXYZ[2]); 478 PetscInfo3(NULL, "Local blocks/processors := %d, %d, %d\n", genCtx.A, genCtx.B, genCtx.C); 479 PetscInfo3(NULL, "Local processors := %d, %d, %d\n", genCtx.M, genCtx.N, genCtx.K); 480 PetscInfo3(NULL, "Local nexyz:= %d, %d, %d\n", genCtx.nex, genCtx.ney, genCtx.nez); 481 PetscInfo3(NULL, "Local delxyz:= %g, %g, %g\n", genCtx.dx, genCtx.dy, genCtx.dz); 482 PetscInfo3(NULL, "Local strides:= %d, %d, %d\n", genCtx.xstride, genCtx.ystride, genCtx.zstride); 483 PetscFunctionReturn(0); 484 } 485 486 /*@C 487 DMMoabCreateBoxMesh - Creates a mesh on the tensor product (box) of intervals with genCtx specified bounds. 488 489 Collective 490 491 Input Parameters: 492 + comm - The communicator for the DM object 493 . dim - The spatial dimension 494 . 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 495 . nele - The number of discrete elements in each direction 496 - user_nghost - The number of ghosted layers needed in the partitioned mesh 497 498 Output Parameter: 499 . dm - The DM object 500 501 Level: beginner 502 503 .seealso: DMSetType(), DMCreate(), DMMoabLoadFromFile() 504 @*/ 505 PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool useSimplex, const PetscReal* bounds, PetscInt nele, PetscInt nghost, DM *dm) 506 { 507 PetscErrorCode ierr; 508 moab::ErrorCode merr; 509 PetscInt a, b, c, n, global_size, global_rank; 510 DM_Moab *dmmoab; 511 moab::Interface *mbImpl; 512 #ifdef MOAB_HAVE_MPI 513 moab::ParallelComm *pcomm; 514 #endif 515 moab::ReadUtilIface *readMeshIface; 516 moab::Range verts, cells, edges, faces, adj, dim3, dim2; 517 DMMoabMeshGeneratorCtx genCtx; 518 const PetscInt npts = nele + 1; /* Number of points in every dimension */ 519 520 moab::Tag global_id_tag, part_tag, geom_tag, mat_tag, dir_tag, neu_tag; 521 moab::Range ownedvtx, ownedelms, localvtxs, localelms; 522 moab::EntityHandle regionset; 523 PetscInt ml = 0, nl = 0, kl = 0; 524 525 PetscFunctionBegin; 526 if (dim < 1 || dim > 3) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension argument for mesh: dim=[1,3].\n"); 527 528 ierr = PetscLogEventRegister("GenerateMesh", DM_CLASSID, &genCtx.generateMesh);CHKERRQ(ierr); 529 ierr = PetscLogEventRegister("AddVertices", DM_CLASSID, &genCtx.generateVertices);CHKERRQ(ierr); 530 ierr = PetscLogEventRegister("AddElements", DM_CLASSID, &genCtx.generateElements);CHKERRQ(ierr); 531 ierr = PetscLogEventRegister("ParResolve", DM_CLASSID, &genCtx.parResolve);CHKERRQ(ierr); 532 ierr = PetscLogEventBegin(genCtx.generateMesh, 0, 0, 0, 0);CHKERRQ(ierr); 533 ierr = MPI_Comm_size(comm, &global_size);CHKERRMPI(ierr); 534 /* total number of vertices in all dimensions */ 535 n = pow(npts, dim); 536 537 /* do some error checking */ 538 if (n < 2) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of points must be >= 2.\n"); 539 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"); 540 if (nghost < 0) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Number of ghost layers cannot be negative.\n"); 541 542 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 543 ierr = DMMoabCreateMoab(comm, NULL, NULL, NULL, dm);CHKERRQ(ierr); 544 545 /* get all the necessary handles from the private DM object */ 546 dmmoab = (DM_Moab*)(*dm)->data; 547 mbImpl = dmmoab->mbiface; 548 #ifdef MOAB_HAVE_MPI 549 pcomm = dmmoab->pcomm; 550 global_rank = pcomm->rank(); 551 #else 552 global_rank = 0; 553 global_size = 1; 554 #endif 555 global_id_tag = dmmoab->ltog_tag; 556 dmmoab->dim = dim; 557 dmmoab->nghostrings = nghost; 558 dmmoab->refct = 1; 559 560 /* create a file set to associate all entities in current mesh */ 561 merr = mbImpl->create_meshset(moab::MESHSET_SET, dmmoab->fileset);MBERR("Creating file set failed", merr); 562 563 /* No errors yet; proceed with building the mesh */ 564 merr = mbImpl->query_interface(readMeshIface);MBERRNM(merr); 565 566 genCtx.M = genCtx.N = genCtx.K = 1; 567 genCtx.A = genCtx.B = genCtx.C = 1; 568 genCtx.blockSizeElementXYZ[0] = 0; 569 genCtx.blockSizeElementXYZ[1] = 0; 570 genCtx.blockSizeElementXYZ[2] = 0; 571 572 ierr = PetscOptionsBegin(comm, "", "DMMoab Creation Options", "DMMOAB");CHKERRQ(ierr); 573 /* Handle DMMoab spatial resolution */ 574 ierr = PetscOptionsInt("-dmb_grid_x", "Number of grid points in x direction", "DMMoabSetSizes", genCtx.blockSizeElementXYZ[0], &genCtx.blockSizeElementXYZ[0], &genCtx.usrxyzgrid);CHKERRQ(ierr); 575 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);} 576 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);} 577 578 /* Handle DMMoab parallel distibution */ 579 ierr = PetscOptionsInt("-dmb_processors_x", "Number of processors in x direction", "DMMoabSetNumProcs", genCtx.M, &genCtx.M, &genCtx.usrprocgrid);CHKERRQ(ierr); 580 if (dim > 1) {ierr = PetscOptionsInt("-dmb_processors_y", "Number of processors in y direction", "DMMoabSetNumProcs", genCtx.N, &genCtx.N, &genCtx.usrprocgrid);CHKERRQ(ierr);} 581 if (dim > 2) {ierr = PetscOptionsInt("-dmb_processors_z", "Number of processors in z direction", "DMMoabSetNumProcs", genCtx.K, &genCtx.K, &genCtx.usrprocgrid);CHKERRQ(ierr);} 582 583 /* Handle DMMoab block refinement */ 584 ierr = PetscOptionsInt("-dmb_refine_x", "Number of refinement blocks in x direction", "DMMoabSetRefinement", genCtx.A, &genCtx.A, &genCtx.usrrefgrid);CHKERRQ(ierr); 585 if (dim > 1) {ierr = PetscOptionsInt("-dmb_refine_y", "Number of refinement blocks in y direction", "DMMoabSetRefinement", genCtx.B, &genCtx.B, &genCtx.usrrefgrid);CHKERRQ(ierr);} 586 if (dim > 2) {ierr = PetscOptionsInt("-dmb_refine_z", "Number of refinement blocks in z direction", "DMMoabSetRefinement", genCtx.C, &genCtx.C, &genCtx.usrrefgrid);CHKERRQ(ierr);} 587 PetscOptionsEnd(); 588 589 ierr = DMMBUtil_InitializeOptions(genCtx, dim, useSimplex, global_rank, global_size, bounds, nele);CHKERRQ(ierr); 590 591 //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); 592 593 if (genCtx.adjEnts) genCtx.keep_skins = true; /* do not delete anything - consumes more memory */ 594 595 /* determine m, n, k for processor rank */ 596 ml = nl = kl = 0; 597 switch (genCtx.dim) { 598 case 1: 599 ml = (genCtx.cumfraction); 600 break; 601 case 2: 602 nl = (genCtx.cumfraction); 603 break; 604 default: 605 kl = (genCtx.cumfraction) / genCtx.q / genCtx.blockSizeElementXYZ[2] / genCtx.C; //genCtx.K 606 break; 607 } 608 609 /* 610 * so there are a total of M * A * blockSizeElement elements in x direction (so M * A * blockSizeElement + 1 verts in x direction) 611 * so there are a total of N * B * blockSizeElement elements in y direction (so N * B * blockSizeElement + 1 verts in y direction) 612 * so there are a total of K * C * blockSizeElement elements in z direction (so K * C * blockSizeElement + 1 verts in z direction) 613 614 * there are ( M * A blockSizeElement) * ( N * B * blockSizeElement) * (K * C * blockSizeElement) hexas 615 * there are ( M * A * blockSizeElement + 1) * ( N * B * blockSizeElement + 1) * (K * C * blockSizeElement + 1) vertices 616 * x is the first dimension that varies 617 */ 618 619 /* generate the block at (a, b, c); it will represent a partition , it will get a partition tag */ 620 PetscInt dum_id = -1; 621 merr = mbImpl->tag_get_handle("GLOBAL_ID", 1, moab::MB_TYPE_INTEGER, global_id_tag);MBERR("Getting Global_ID Tag handle failed", merr); 622 623 merr = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, mat_tag);MBERR("Getting Material set Tag handle failed", merr); 624 merr = mbImpl->tag_get_handle(DIRICHLET_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, dir_tag);MBERR("Getting Dirichlet set Tag handle failed", merr); 625 merr = mbImpl->tag_get_handle(NEUMANN_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, neu_tag);MBERR("Getting Neumann set Tag handle failed", merr); 626 627 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); 628 629 /* lets create some sets */ 630 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); 631 merr = mbImpl->create_meshset(moab::MESHSET_SET, regionset);MBERRNM(merr); 632 ierr = PetscLogEventEnd(genCtx.generateMesh, 0, 0, 0, 0);CHKERRQ(ierr); 633 634 for (a = 0; a < (genCtx.dim > 0 ? genCtx.A : genCtx.A); a++) { 635 for (b = 0; b < (genCtx.dim > 1 ? genCtx.B : 1); b++) { 636 for (c = 0; c < (genCtx.dim > 2 ? genCtx.C : 1); c++) { 637 638 moab::EntityHandle startv; 639 640 ierr = PetscLogEventBegin(genCtx.generateVertices, 0, 0, 0, 0);CHKERRQ(ierr); 641 ierr = DMMoab_GenerateVertices_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, verts);CHKERRQ(ierr); 642 ierr = PetscLogEventEnd(genCtx.generateVertices, 0, 0, 0, 0);CHKERRQ(ierr); 643 644 ierr = PetscLogEventBegin(genCtx.generateElements, 0, 0, 0, 0);CHKERRQ(ierr); 645 ierr = DMMoab_GenerateElements_Private(mbImpl, readMeshIface, genCtx, ml, nl, kl, a, b, c, global_id_tag, startv, cells);CHKERRQ(ierr); 646 ierr = PetscLogEventEnd(genCtx.generateElements, 0, 0, 0, 0);CHKERRQ(ierr); 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);MBERR("Can't create mesh set.", merr); 661 662 merr = mbImpl->add_entities(part_set, verts);MBERR("Can't add vertices to set.", merr); 663 merr = mbImpl->add_entities(part_set, cells);MBERR("Can't add entities to set.", merr); 664 merr = mbImpl->add_entities(regionset, cells);MBERR("Can't add entities to set.", merr); 665 666 /* if needed, add all edges and faces */ 667 if (genCtx.adjEnts) 668 { 669 if (genCtx.dim > 1) { 670 merr = mbImpl->get_adjacencies(cells, 1, true, edges, moab::Interface::UNION);MBERR("Can't get edges", merr); 671 merr = mbImpl->add_entities(part_set, edges);MBERR("Can't add edges to partition set.", merr); 672 } 673 if (genCtx.dim > 2) { 674 merr = mbImpl->get_adjacencies(cells, 2, true, faces, moab::Interface::UNION);MBERR("Can't get faces", merr); 675 merr = mbImpl->add_entities(part_set, faces);MBERR("Can't add faces to partition set.", merr); 676 } 677 edges.clear(); 678 faces.clear(); 679 } 680 verts.clear(); cells.clear(); 681 682 merr = mbImpl->tag_set_data(part_tag, &part_set, 1, &part_num);MBERR("Can't set part tag on set", merr); 683 if (dmmoab->fileset) { 684 merr = mbImpl->add_parent_child(dmmoab->fileset, part_set);MBERR("Can't add part set to file set.", merr); 685 merr = mbImpl->unite_meshset(dmmoab->fileset, part_set);MBERRNM(merr); 686 } 687 merr = mbImpl->add_entities(dmmoab->fileset, &part_set, 1);MBERRNM(merr); 688 } 689 } 690 } 691 692 merr = mbImpl->add_parent_child(dmmoab->fileset, regionset);MBERRNM(merr); 693 694 /* Only in parallel: resolve shared entities between processors and exchange ghost layers */ 695 if (global_size > 1) { 696 697 ierr = PetscLogEventBegin(genCtx.parResolve, 0, 0, 0, 0);CHKERRQ(ierr); 698 699 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, genCtx.dim, cells);MBERR("Can't get all d-dimensional elements.", merr); 700 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 0, verts);MBERR("Can't get all vertices.", merr); 701 702 if (genCtx.A * genCtx.B * genCtx.C != 1) { // merge needed 703 moab::MergeMesh mm(mbImpl); 704 if (genCtx.newMergeMethod) { 705 merr = mm.merge_using_integer_tag(verts, global_id_tag);MBERR("Can't merge with GLOBAL_ID tag", merr); 706 } 707 else { 708 merr = mm.merge_entities(cells, 0.0001);MBERR("Can't merge with coordinates", merr); 709 } 710 } 711 712 #ifdef MOAB_HAVE_MPI 713 /* check the handles */ 714 merr = pcomm->check_all_shared_handles();MBERRV(mbImpl, merr); 715 716 /* resolve the shared entities by exchanging information to adjacent processors */ 717 merr = pcomm->resolve_shared_ents(dmmoab->fileset, cells, dim, dim - 1, NULL, &global_id_tag);MBERRV(mbImpl, merr); 718 if (dmmoab->fileset) { 719 merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false, &dmmoab->fileset);MBERRV(mbImpl, merr); 720 } 721 else { 722 merr = pcomm->exchange_ghost_cells(dim, 0, nghost, dim, true, false);MBERRV(mbImpl, merr); 723 } 724 725 /* Reassign global IDs on all entities. */ 726 merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, false, true, false);MBERRNM(merr); 727 #endif 728 729 ierr = PetscLogEventEnd(genCtx.parResolve, 0, 0, 0, 0);CHKERRQ(ierr); 730 } 731 732 if (!genCtx.keep_skins) { // default is to delete the 1- and 2-dimensional entities 733 // delete all quads and edges 734 moab::Range toDelete; 735 if (genCtx.dim > 1) { 736 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 1, toDelete);MBERR("Can't get edges", merr); 737 } 738 739 if (genCtx.dim > 2) { 740 merr = mbImpl->get_entities_by_dimension(dmmoab->fileset, 2, toDelete);MBERR("Can't get faces", merr); 741 } 742 743 #ifdef MOAB_HAVE_MPI 744 merr = dmmoab->pcomm->delete_entities(toDelete) ;MBERR("Can't delete entities", merr); 745 #endif 746 } 747 748 /* set geometric dimension tag for regions */ 749 merr = mbImpl->tag_set_data(geom_tag, ®ionset, 1, &dmmoab->dim);MBERRNM(merr); 750 /* set default material ID for regions */ 751 int default_material = 1; 752 merr = mbImpl->tag_set_data(mat_tag, ®ionset, 1, &default_material);MBERRNM(merr); 753 /* 754 int default_dbc = 0; 755 merr = mbImpl->tag_set_data(dir_tag, &vertexset, 1, &default_dbc);MBERRNM(merr); 756 */ 757 PetscFunctionReturn(0); 758 } 759 760 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) 761 { 762 char *ropts; 763 char ropts_par[PETSC_MAX_PATH_LEN], ropts_pargh[PETSC_MAX_PATH_LEN]; 764 char ropts_dbg[PETSC_MAX_PATH_LEN]; 765 PetscErrorCode ierr; 766 767 PetscFunctionBegin; 768 ierr = PetscMalloc1(PETSC_MAX_PATH_LEN, &ropts);CHKERRQ(ierr); 769 ierr = PetscMemzero(&ropts_par, PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 770 ierr = PetscMemzero(&ropts_pargh, PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 771 ierr = PetscMemzero(&ropts_dbg, PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 772 773 /* do parallel read unless using only one processor */ 774 if (numproc > 1) { 775 // 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); 776 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); 777 if (nghost) { 778 ierr = PetscSNPrintf(ropts_pargh, PETSC_MAX_PATH_LEN, "PARALLEL_GHOSTS=%d.0.%d;", dim, nghost);CHKERRQ(ierr); 779 } 780 } 781 782 if (dbglevel) { 783 if (numproc > 1) { 784 ierr = PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%d;DEBUG_PIO=%d;", dbglevel, dbglevel);CHKERRQ(ierr); 785 } 786 else { 787 ierr = PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "CPUTIME;DEBUG_IO=%d;", dbglevel);CHKERRQ(ierr); 788 } 789 } 790 791 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); 792 *read_opts = ropts; 793 PetscFunctionReturn(0); 794 } 795 796 /*@C 797 DMMoabLoadFromFile - Creates a DM object by loading the mesh from a user specified file. 798 799 Collective 800 801 Input Parameters: 802 + comm - The communicator for the DM object 803 . dim - The spatial dimension 804 . filename - The name of the mesh file to be loaded 805 - usrreadopts - The options string to read a MOAB mesh. 806 807 Reference (Parallel Mesh Initialization: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html#fivetwo) 808 809 Output Parameter: 810 . dm - The DM object 811 812 Level: beginner 813 814 .seealso: DMSetType(), DMCreate(), DMMoabCreateBoxMesh() 815 @*/ 816 PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm, PetscInt dim, PetscInt nghost, const char* filename, const char* usrreadopts, DM *dm) 817 { 818 moab::ErrorCode merr; 819 PetscInt nprocs; 820 DM_Moab *dmmoab; 821 moab::Interface *mbiface; 822 #ifdef MOAB_HAVE_MPI 823 moab::ParallelComm *pcomm; 824 #endif 825 moab::Range verts, elems; 826 const char *readopts; 827 PetscErrorCode ierr; 828 829 PetscFunctionBegin; 830 PetscValidPointer(dm, 6); 831 832 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 833 ierr = DMMoabCreateMoab(comm, NULL, NULL, NULL, dm);CHKERRQ(ierr); 834 835 /* get all the necessary handles from the private DM object */ 836 dmmoab = (DM_Moab*)(*dm)->data; 837 mbiface = dmmoab->mbiface; 838 #ifdef MOAB_HAVE_MPI 839 pcomm = dmmoab->pcomm; 840 nprocs = pcomm->size(); 841 #else 842 nprocs = 1; 843 #endif 844 /* TODO: Decipher dimension based on the loaded mesh instead of getting from user */ 845 dmmoab->dim = dim; 846 dmmoab->nghostrings = nghost; 847 dmmoab->refct = 1; 848 849 /* create a file set to associate all entities in current mesh */ 850 merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, dmmoab->fileset);MBERR("Creating file set failed", merr); 851 852 /* add mesh loading options specific to the DM */ 853 ierr = DMMoab_GetReadOptions_Private(dmmoab->partition_by_rank, nprocs, dim, nghost, dmmoab->read_mode, 854 dmmoab->rw_dbglevel, dmmoab->extra_read_options, usrreadopts, &readopts);CHKERRQ(ierr); 855 856 PetscInfo2(*dm, "Reading file %s with options: %s\n", filename, readopts); 857 858 /* Load the mesh from a file. */ 859 if (dmmoab->fileset) { 860 merr = mbiface->load_file(filename, &dmmoab->fileset, readopts);MBERRVM(mbiface, "Reading MOAB file failed.", merr); 861 } 862 else { 863 merr = mbiface->load_file(filename, 0, readopts);MBERRVM(mbiface, "Reading MOAB file failed.", merr); 864 } 865 866 #ifdef MOAB_HAVE_MPI 867 /* Reassign global IDs on all entities. */ 868 /* merr = pcomm->assign_global_ids(dmmoab->fileset, dim, 1, true, true, true);MBERRNM(merr); */ 869 #endif 870 871 /* load the local vertices */ 872 merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);MBERRNM(merr); 873 /* load the local elements */ 874 merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);MBERRNM(merr); 875 876 #ifdef MOAB_HAVE_MPI 877 /* Everything is set up, now just do a tag exchange to update tags 878 on all of the ghost vertexes */ 879 merr = pcomm->exchange_tags(dmmoab->ltog_tag, verts);MBERRV(mbiface, merr); 880 merr = pcomm->exchange_tags(dmmoab->ltog_tag, elems);MBERRV(mbiface, merr); 881 merr = pcomm->collective_sync_partition();MBERR("Collective sync failed", merr); 882 #endif 883 884 PetscInfo3(*dm, "MOAB file '%s' was successfully loaded. Found %D vertices and %D elements.\n", filename, verts.size(), elems.size()); 885 ierr = PetscFree(readopts);CHKERRQ(ierr); 886 PetscFunctionReturn(0); 887 } 888 889 /*@C 890 DMMoabRenumberMeshEntities - Order and number all entities (vertices->elements) to be contiguously ordered 891 in parallel 892 893 Collective 894 895 Input Parameters: 896 . dm - The DM object 897 898 Level: advanced 899 900 .seealso: DMSetUp(), DMCreate() 901 @*/ 902 PetscErrorCode DMMoabRenumberMeshEntities(DM dm) 903 { 904 moab::Range verts; 905 906 PetscFunctionBegin; 907 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 908 909 #ifdef MOAB_HAVE_MPI 910 /* Insert new points */ 911 moab::ErrorCode merr; 912 merr = ((DM_Moab*) dm->data)->pcomm->assign_global_ids(((DM_Moab*) dm->data)->fileset, 3, 0, false, true, false);MBERRNM(merr); 913 #endif 914 PetscFunctionReturn(0); 915 } 916 917