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