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