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/ScdInterface.hpp> 8 #include <moab/CN.hpp> 9 10 11 static PetscErrorCode DMMoabComputeDomainBounds_Private(moab::ParallelComm* pcomm, PetscInt dim, PetscInt neleglob, PetscInt *ise) 12 { 13 PetscInt size,rank; 14 PetscInt fraction,remainder; 15 PetscInt starts[3],sizes[3]; 16 17 PetscFunctionBegin; 18 if(dim<1 && dim>3) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The problem dimension is invalid: %D",dim); 19 20 size=pcomm->size(); 21 rank=pcomm->rank(); 22 fraction=neleglob/size; /* partition only by the largest dimension */ 23 remainder=neleglob%size; /* remainder after partition which gets evenly distributed by round-robin */ 24 25 if(fraction==0) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The leading dimension size should be greater than number of processors: %D < %D",neleglob,size); 26 27 starts[0]=starts[1]=starts[2]=0; /* default dimensional element = 1 */ 28 sizes[0]=sizes[1]=sizes[2]=neleglob; /* default dimensional element = 1 */ 29 30 if(rank < remainder) { /* This process gets "fraction+1" elements */ 31 sizes[dim-1] = (fraction + 1); 32 starts[dim-1] = rank * (fraction+1); 33 } else { /* This process gets "fraction" elements */ 34 sizes[dim-1] = fraction; 35 starts[dim-1] = (remainder*(fraction+1) + fraction*(rank-remainder)); 36 } 37 38 for(int i=dim-1; i>=0; --i) { 39 ise[2*i]=starts[i];ise[2*i+1]=starts[i]+sizes[i]; 40 } 41 PetscFunctionReturn(0); 42 } 43 44 static void DMMoab_SetStructuredCoords_Private(PetscInt i, PetscInt j, PetscInt k, PetscReal hx, PetscReal hy, PetscReal hz, PetscInt vcount, std::vector<double*>& vcoords) 45 { 46 vcoords[0][vcount] = i*hx; 47 vcoords[1][vcount] = j*hy; 48 vcoords[2][vcount] = k*hz; 49 } 50 51 static void DMMoab_SetTensorElementConnectivity_Private(PetscInt dim, moab::EntityType etype, PetscInt offset, PetscInt nele, PetscInt i, PetscInt j, PetscInt k, PetscInt vfirst, moab::EntityHandle *connectivity) 52 { 53 std::vector<int> subent_conn(pow(2,dim)); /* only linear edge, quad, hex supported now */ 54 55 moab::CN::SubEntityVertexIndices(etype, dim, 0, subent_conn.data()); 56 57 switch(dim) { 58 case 1: 59 connectivity[offset+subent_conn[0]] = vfirst+i; 60 connectivity[offset+subent_conn[1]] = vfirst+(i+1); 61 break; 62 case 2: 63 connectivity[offset+subent_conn[0]] = vfirst+i+j*(nele+1); 64 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+j*(nele+1); 65 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(j+1)*(nele+1); 66 connectivity[offset+subent_conn[3]] = vfirst+i+(j+1)*(nele+1); 67 break; 68 case 3: 69 default: 70 connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k); 71 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*k); 72 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k); 73 connectivity[offset+subent_conn[3]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*k); 74 connectivity[offset+subent_conn[4]] = vfirst+i+(nele+1)*(j+(nele+1)*(k+1)); 75 connectivity[offset+subent_conn[5]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1)); 76 connectivity[offset+subent_conn[6]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*(k+1)); 77 connectivity[offset+subent_conn[7]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1)); 78 break; 79 } 80 } 81 82 static void DMMoab_SetSimplexElementConnectivity_Private(PetscInt dim, PetscInt subelem, moab::EntityType etype, PetscInt offset, PetscInt nele, PetscInt i, PetscInt j, PetscInt k, PetscInt vfirst, moab::EntityHandle *connectivity) 83 { 84 std::vector<int> subent_conn(pow(2,dim)); /* only linear edge, quad, hex supported now */ 85 86 moab::CN::SubEntityVertexIndices(etype, dim, 0, subent_conn.data()); 87 88 switch(dim) { 89 case 1: 90 connectivity[offset+subent_conn[0]] = vfirst+i; 91 connectivity[offset+subent_conn[1]] = vfirst+(i+1); 92 break; 93 case 2: 94 if (subelem) { /* 1 2 3 of a QUAD */ 95 connectivity[offset+subent_conn[0]] = vfirst+i+j*(nele+1); 96 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+j*(nele+1); 97 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(j+1)*(nele+1); 98 } 99 else { /* 3 4 1 of a QUAD */ 100 connectivity[offset+subent_conn[0]] = vfirst+(i+1)+(j+1)*(nele+1); 101 connectivity[offset+subent_conn[1]] = vfirst+i+(j+1)*(nele+1); 102 connectivity[offset+subent_conn[2]] = vfirst+i+j*(nele+1); 103 } 104 break; 105 case 3: 106 default: 107 switch(subelem) { 108 case 0: /* 0 1 2 5 of a HEX */ 109 connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k); 110 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*k); 111 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k); 112 connectivity[offset+subent_conn[3]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1)); 113 break; 114 case 1: /* 0 2 7 5 of a HEX */ 115 connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k); 116 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k); 117 connectivity[offset+subent_conn[2]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1)); 118 connectivity[offset+subent_conn[3]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1)); 119 break; 120 case 2: /* 0 2 3 7 of a HEX */ 121 connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k); 122 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k); 123 connectivity[offset+subent_conn[2]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*k); 124 connectivity[offset+subent_conn[3]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1)); 125 break; 126 case 3: /* 0 5 7 4 of a HEX */ 127 connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k); 128 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1)); 129 connectivity[offset+subent_conn[2]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1)); 130 connectivity[offset+subent_conn[3]] = vfirst+i+(nele+1)*(j+(nele+1)*(k+1)); 131 break; 132 case 4: /* 2 7 5 6 of a HEX */ 133 connectivity[offset+subent_conn[0]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k); 134 connectivity[offset+subent_conn[1]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1)); 135 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1)); 136 connectivity[offset+subent_conn[3]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*(k+1)); 137 break; 138 } 139 break; 140 } 141 } 142 143 static void DMMoab_SetElementConnectivity_Private(PetscBool useSimplex, PetscInt dim, moab::EntityType etype, PetscInt *ecount, PetscInt vpere, PetscInt nele, PetscInt i, PetscInt j, PetscInt k, PetscInt vfirst, moab::EntityHandle *connectivity) 144 { 145 PetscInt m,subelem; 146 if (useSimplex) { 147 subelem=(dim==1 ? 1 : (dim==2 ? 2 : 5)); 148 for (m=0; m<subelem; m++) 149 DMMoab_SetSimplexElementConnectivity_Private(dim, m, etype, (*ecount+m)*vpere, nele, i, j, k, vfirst, connectivity); 150 151 } 152 else { 153 subelem=1; 154 DMMoab_SetTensorElementConnectivity_Private(dim, etype, (*ecount)*vpere, nele, i, j, k, vfirst, connectivity); 155 } 156 *ecount+=subelem; 157 } 158 159 160 /*@ 161 DMMoabCreateBoxMesh - Creates a mesh on the tensor product (box) of intervals with user specified bounds. 162 163 Collective on MPI_Comm 164 165 Input Parameters: 166 + comm - The communicator for the DM object 167 . dim - The spatial dimension 168 . 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 169 . nele - The number of discrete elements in each direction 170 . user_nghost - The number of ghosted layers needed in the partitioned mesh 171 172 Output Parameter: 173 . dm - The DM object 174 175 Level: beginner 176 177 .keywords: DM, create 178 .seealso: DMSetType(), DMCreate(), DMMoabLoadFromFile() 179 @*/ 180 PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool useSimplex, const PetscReal* bounds, PetscInt nele, PetscInt user_nghost, DM *dm) 181 { 182 PetscErrorCode ierr; 183 moab::ErrorCode merr; 184 PetscInt i,j,k,n,nprocs,rank=-1,ecount=0,vcount=0; 185 DM_Moab *dmmoab; 186 moab::Interface *mbiface; 187 moab::ParallelComm *pcomm; 188 moab::ReadUtilIface *readMeshIface; 189 190 moab::Tag id_tag=PETSC_NULL,part_tag=PETSC_NULL; 191 moab::Range ownedvtx,ownedelms,localvtxs,localelms; 192 moab::EntityHandle vfirst,efirst,regionset,faceset,edgeset,vtxset,partset; 193 std::vector<double*> vcoords; 194 moab::EntityHandle *connectivity = 0; 195 moab::EntityType etype=moab::MBHEX; 196 PetscInt ise[6]; 197 std::vector<int> vgid; 198 PetscReal xse[6],defbounds[6]; 199 200 /* TODO: Fix nghost > 0 - now relying on exchange_ghost_cells */ 201 const PetscInt nghost=0; 202 203 moab::Tag geom_tag; 204 205 moab::Range adj,dim3,dim2; 206 bool build_adjacencies=false; 207 208 const PetscInt npts=nele+1; /* Number of points in every dimension */ 209 PetscInt vpere=0,locnele=0,locnpts=0,ghnele,ghnpts; /* Number of verts/element, vertices, elements owned by this process */ 210 211 PetscFunctionBegin; 212 if(dim < 1 || dim > 3) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Invalid dimension argument for mesh: dim=[1,3].\n"); 213 214 ierr = MPI_Comm_size(comm, &nprocs);CHKERRQ(ierr); 215 /* total number of vertices in all dimensions */ 216 n=pow(npts,dim); 217 218 /* do some error checking */ 219 if(n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be >= 2.\n"); 220 if(nprocs > n) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of processors must be less than or equal to number of elements.\n"); 221 if(nghost < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of ghost layers cannot be negative.\n"); 222 223 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 224 ierr = DMMoabCreateMoab(comm, NULL, NULL, NULL, NULL, dm);CHKERRQ(ierr); 225 226 /* get all the necessary handles from the private DM object */ 227 dmmoab = (DM_Moab*)(*dm)->data; 228 mbiface = dmmoab->mbiface; 229 pcomm = dmmoab->pcomm; 230 id_tag = dmmoab->ltog_tag; 231 nprocs = pcomm->size(); 232 dmmoab->dim = dim; 233 234 /* create a file set to associate all entities in current mesh */ 235 merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, dmmoab->fileset);MBERR("Creating file set failed", merr); 236 237 /* No errors yet; proceed with building the mesh */ 238 merr = mbiface->query_interface(readMeshIface);MBERRNM(merr); 239 240 ierr = PetscMemzero(ise,sizeof(PetscInt)*6);CHKERRQ(ierr); 241 242 /* call the collective routine that computes the domain bounds for a structured mesh using MOAB */ 243 ierr = DMMoabComputeDomainBounds_Private(pcomm, dim, nele, ise);CHKERRQ(ierr); 244 245 /* set some variables based on dimension */ 246 switch(dim) { 247 case 1: 248 vpere = 2; 249 locnele = (ise[1]-ise[0]); 250 locnpts = (ise[1]-ise[0]+1); 251 ghnele = (nghost > 0 ? (ise[0] > nghost ? 1 : 0) + (ise[1] < nele - nghost ? 1 : 0) : 0 ); 252 ghnpts = (nghost > 0 ? (ise[0] > 0 ? 1 : 0) + (ise[1] < nele ? 1 : 0) : 0); 253 etype = moab::MBEDGE; 254 break; 255 case 2: 256 locnpts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1); 257 ghnpts = (nghost > 0 ? (ise[2] > 0 ? npts : 0) + (ise[3] < nele ? npts : 0) : 0); 258 if (useSimplex) { 259 vpere = 3; 260 locnele = 2*(ise[1]-ise[0])*(ise[3]-ise[2]); 261 ghnele = 2*(nghost > 0 ? (ise[2] > 0 ? nele : 0) + (ise[3] < nele ? nele : 0) : 0); 262 etype = moab::MBTRI; 263 } 264 else { 265 vpere = 4; 266 locnele = (ise[1]-ise[0])*(ise[3]-ise[2]); 267 ghnele = (nghost > 0 ? (ise[2] > 0 ? nele : 0) + (ise[3] < nele ? nele : 0) : 0); 268 etype = moab::MBQUAD; 269 } 270 break; 271 case 3: 272 default: 273 locnpts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1)*(ise[5]-ise[4]+1); 274 ghnpts = (nghost > 0 ? (ise[4] > 0 ? npts*npts : 0) + (ise[5] < nele ? npts*npts : 0) : 0); 275 if (useSimplex) { 276 vpere = 4; 277 locnele = 5*(ise[1]-ise[0])*(ise[3]-ise[2])*(ise[5]-ise[4]); 278 ghnele = 5*(nghost > 0 ? (ise[4] > 0 ? nele*nele : 0) + (ise[5] < nele ? nele*nele : 0) : 0); 279 etype = moab::MBTET; 280 } 281 else { 282 vpere = 8; 283 locnele = (ise[1]-ise[0])*(ise[3]-ise[2])*(ise[5]-ise[4]); 284 ghnele = (nghost > 0 ? (ise[4] > 0 ? nele*nele : 0) + (ise[5] < nele ? nele*nele : 0) : 0); 285 etype = moab::MBHEX; 286 } 287 break; 288 } 289 290 /* we have a domain of size [1,1,1] - now compute local co-ordinate box */ 291 ierr = PetscMemzero(xse,sizeof(PetscReal)*6);CHKERRQ(ierr); 292 for(i=0; i<6; ++i) { 293 xse[i]=(PetscReal)ise[i]/nele; 294 } 295 296 /* Create vertexes and set the coodinate of each vertex */ 297 merr = readMeshIface->get_node_coords(3,locnpts+ghnpts,0,vfirst,vcoords,n);MBERRNM(merr); 298 299 /* Compute the co-ordinates of vertices and global IDs */ 300 vgid.resize(locnpts+ghnpts); 301 vcount=0; 302 if (!bounds) { /* default box mesh is defined on a unit-cube */ 303 defbounds[0]=0.0; defbounds[1]=1.0; 304 defbounds[2]=0.0; defbounds[3]=1.0; 305 defbounds[4]=0.0; defbounds[5]=1.0; 306 bounds=defbounds; 307 } 308 else { 309 /* validate the bounds data */ 310 if(bounds[0] >= bounds[1]) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"X-dim: Left boundary cannot be greater than right. [%G >= %G]\n",bounds[0],bounds[1]); 311 if(dim > 1 && (bounds[2] >= bounds[3])) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Y-dim: Left boundary cannot be greater than right. [%G >= %G]\n",bounds[2],bounds[3]); 312 if(dim > 2 && (bounds[4] >= bounds[5])) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Z-dim: Left boundary cannot be greater than right. [%G >= %G]\n",bounds[4],bounds[5]); 313 } 314 315 const double hx=(bounds[1]-bounds[0])/nele; 316 const double hy=(dim > 1 ? (bounds[3]-bounds[2])/nele : 0.0); 317 const double hz=(dim > 2 ? (bounds[5]-bounds[4])/nele : 0.0); 318 319 /* create all the owned vertices */ 320 for (k = ise[4]; k <= ise[5]; k++) { 321 for (j = ise[2]; j <= ise[3]; j++) { 322 for (i = ise[0]; i <= ise[1]; i++, vcount++) { 323 DMMoab_SetStructuredCoords_Private(i,j,k,hx,hy,hz,vcount,vcoords); 324 vgid[vcount] = (k*npts+j)*npts+i+1; 325 } 326 } 327 } 328 329 /* create ghosted vertices requested by user - below the current plane */ 330 if (ise[2*dim-2] > 0) { 331 for (k = (dim==3?ise[4]-nghost:ise[4]); k <= (dim==3?ise[4]-1:ise[5]); k++) { 332 for (j = (dim==2?ise[2]-nghost:ise[2]); j <= (dim==2?ise[2]-1:ise[3]); j++) { 333 for (i = (dim>1?ise[0]:ise[0]-nghost); i <= (dim>1?ise[1]:ise[0]-1); i++, vcount++) { 334 DMMoab_SetStructuredCoords_Private(i,j,k,hx,hy,hz,vcount,vcoords); 335 vgid[vcount] = (k*npts+j)*npts+i+1; 336 } 337 } 338 } 339 } 340 341 /* create ghosted vertices requested by user - above the current plane */ 342 if (ise[2*dim-1] < nele) { 343 for (k = (dim==3?ise[5]+1:ise[4]); k <= (dim==3?ise[5]+nghost:ise[5]); k++) { 344 for (j = (dim==2?ise[3]+1:ise[2]); j <= (dim==2?ise[3]+nghost:ise[3]); j++) { 345 for (i = (dim>1?ise[0]:ise[1]+1); i <= (dim>1?ise[1]:ise[1]+nghost); i++, vcount++) { 346 DMMoab_SetStructuredCoords_Private(i,j,k,hx,hy,hz,vcount,vcoords); 347 vgid[vcount] = (k*npts+j)*npts+i+1; 348 } 349 } 350 } 351 } 352 353 if (locnpts+ghnpts != vcount) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of vertices! (%D!=%D)",locnpts+ghnpts,vcount); 354 355 merr = mbiface->get_entities_by_type(0,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr); 356 357 /* The global ID tag is applied to each owned 358 vertex. It acts as an global identifier which MOAB uses to 359 assemble the individual pieces of the mesh: 360 Set the global ID indices */ 361 merr = mbiface->tag_set_data(id_tag,ownedvtx,vgid.data());MBERRNM(merr); 362 363 /* Create elements between mesh points using the ReadUtilInterface 364 get the reference to element connectivities for all local elements from the ReadUtilInterface */ 365 merr = readMeshIface->get_element_connect (locnele+ghnele,vpere,etype,1,efirst,connectivity);MBERRNM(merr); 366 367 /* offset appropriately so that only local ID and not global ID numbers are set for connectivity array */ 368 vfirst-=vgid[0]-1; 369 370 /* 3. Loop over elements in 3 nested loops over i, j, k; for each (i,j,k): 371 and then set the connectivity for each element appropriately */ 372 ecount=0; 373 /* create ghosted elements requested by user - below the current plane */ 374 if (ise[2*dim-2] >= nghost) { 375 for (k = (dim==3?ise[4]-nghost:ise[4]); k < (dim==3?ise[4]:std::max(ise[5],1)); k++) { 376 for (j = (dim==2?ise[2]-nghost:ise[2]); j < (dim==2?ise[2]:std::max(ise[3],1)); j++) { 377 for (i = (dim>1?ise[0]:ise[0]-nghost); i < (dim>1?std::max(ise[1],1):ise[0]); i++) { 378 DMMoab_SetElementConnectivity_Private(useSimplex, dim, etype, &ecount, vpere, nele, i, j, k, vfirst, connectivity); 379 } 380 } 381 } 382 } 383 384 /* create owned elements requested by user */ 385 for (k = ise[4]; k < std::max(ise[5],1); k++) { 386 for (j = ise[2]; j < std::max(ise[3],1); j++) { 387 for (i = ise[0]; i < std::max(ise[1],1); i++) { 388 DMMoab_SetElementConnectivity_Private(useSimplex, dim, etype, &ecount, vpere, nele, i, j, k, vfirst, connectivity); 389 } 390 } 391 } 392 393 /* create ghosted elements requested by user - above the current plane */ 394 if (ise[2*dim-1] <= nele-nghost) { 395 for (k = (dim==3?ise[5]:ise[4]); k < (dim==3?ise[5]+nghost:std::max(ise[5],1)); k++) { 396 for (j = (dim==2?ise[3]:ise[2]); j < (dim==2?ise[3]+nghost:std::max(ise[3],1)); j++) { 397 for (i = (dim>1?ise[0]:ise[1]); i < (dim>1?std::max(ise[1],1):ise[1]+nghost); i++) { 398 DMMoab_SetElementConnectivity_Private(useSimplex, dim, etype, &ecount, vpere, nele, i, j, k, vfirst, connectivity); 399 } 400 } 401 } 402 } 403 404 merr = readMeshIface->update_adjacencies(efirst,locnele+ghnele,vpere,connectivity);MBERRNM(merr); 405 406 /* 2. Get the vertices and hexes from moab and check their numbers against I*J*K and (I-1)*(J-1)*(K-1), resp. 407 first '0' specifies "root set", or entire MOAB instance, second the entity dimension being requested */ 408 merr = mbiface->get_entities_by_dimension(0, dim, ownedelms);MBERRNM(merr); 409 410 if (locnele+ghnele != (int) ownedelms.size()) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of elements! (%D!=%D)",locnele+ghnele,ownedelms.size()); 411 else if(locnpts+ghnpts != (int) ownedvtx.size()) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of vertices! (%D!=%D)",locnpts+ghnpts,ownedvtx.size()); 412 else { 413 ierr = PetscInfo2(NULL, "Created %D elements and %D vertices.\n", ownedelms.size(), ownedvtx.size());CHKERRQ(ierr); 414 } 415 416 /* lets create some sets */ 417 merr = mbiface->tag_get_handle(GEOM_DIMENSION_TAG_NAME, 1, moab::MB_TYPE_INTEGER, geom_tag, moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT);MBERRNM(merr); 418 419 merr = mbiface->create_meshset(moab::MESHSET_SET, regionset);MBERRNM(merr); 420 merr = mbiface->add_entities(regionset, ownedelms);MBERRNM(merr); 421 merr = mbiface->tag_set_data(geom_tag, ®ionset, 1, &dmmoab->dim);MBERRNM(merr); 422 merr = mbiface->add_parent_child(dmmoab->fileset,regionset);MBERRNM(merr); 423 merr = mbiface->unite_meshset(dmmoab->fileset, regionset);MBERRNM(merr); 424 425 merr = mbiface->create_meshset(moab::MESHSET_SET, vtxset);MBERRNM(merr); 426 merr = mbiface->add_entities(vtxset, ownedvtx);MBERRNM(merr); 427 merr = mbiface->add_parent_child(dmmoab->fileset,vtxset);MBERRNM(merr); 428 merr = mbiface->unite_meshset(dmmoab->fileset, vtxset);MBERRNM(merr); 429 430 if (build_adjacencies) { 431 // generate all lower dimensional adjacencies 432 merr = mbiface->get_adjacencies( ownedelms, dim-1, true, adj, moab::Interface::UNION );MBERRNM(merr); 433 merr = dmmoab->pcomm->get_part_entities(dim2, dim-1);MBERRNM(merr); 434 adj.merge(dim2); 435 436 /* create face sets */ 437 merr = mbiface->create_meshset(moab::MESHSET_SET, faceset);MBERRNM(merr); 438 merr = mbiface->add_entities(faceset, adj);MBERRNM(merr); 439 merr = mbiface->add_parent_child(dmmoab->fileset,faceset);MBERRNM(merr); 440 i=dim-1; 441 merr = mbiface->tag_set_data(geom_tag, &faceset, 1, &i);MBERRNM(merr); 442 merr = mbiface->unite_meshset(dmmoab->fileset, faceset);MBERRNM(merr); 443 PetscInfo2(NULL, "Found %d %d-Dim quantities.\n", adj.size(), dim-1); 444 445 if (dim > 2) { 446 dim2.clear(); 447 /* create edge sets, if appropriate i.e., if dim=3 */ 448 merr = mbiface->create_meshset(moab::MESHSET_SET, edgeset);MBERRNM(merr); 449 merr = mbiface->get_adjacencies(adj, dim-1, true, dim2, moab::Interface::UNION );MBERRNM(merr); 450 merr = mbiface->add_entities(edgeset, dim2);MBERRNM(merr); 451 merr = mbiface->add_parent_child(dmmoab->fileset,edgeset);MBERRNM(merr); 452 i=dim-2; 453 merr = mbiface->tag_set_data(geom_tag, &edgeset, 1, &i);MBERRNM(merr); 454 merr = mbiface->unite_meshset(dmmoab->fileset, edgeset);MBERRNM(merr); 455 PetscInfo2(NULL, "Found %d %d-Dim quantities.\n", adj.size(), dim-2); 456 } 457 } 458 459 /* check the handles */ 460 merr = pcomm->check_all_shared_handles();MBERRV(mbiface,merr); 461 462 /* resolve the shared entities by exchanging information to adjacent processors */ 463 merr = mbiface->get_entities_by_type(dmmoab->fileset,etype,ownedelms,true);MBERRNM(merr); 464 merr = pcomm->resolve_shared_ents(dmmoab->fileset,ownedelms,dim,dim-1,NULL,&id_tag);MBERRV(mbiface,merr); 465 466 merr = pcomm->exchange_ghost_cells(dim,0,user_nghost,dim,true,false,&dmmoab->fileset);MBERRV(mbiface,merr); 467 468 /* Reassign global IDs on all entities. */ 469 merr = pcomm->assign_global_ids(dmmoab->fileset,dim,1,false,true,false);MBERRNM(merr); 470 471 /* filter based on parallel status */ 472 merr = dmmoab->pcomm->filter_pstatus(ownedvtx,PSTATUS_NOT_OWNED,PSTATUS_NOT,-1,&localvtxs);MBERRNM(merr); 473 merr = dmmoab->pcomm->filter_pstatus(ownedelms,PSTATUS_NOT_OWNED,PSTATUS_NOT,-1,&localelms);MBERRNM(merr); 474 475 /* set the parallel partition tag data */ 476 merr = mbiface->tag_get_handle("PARALLEL_PARTITION", 1, moab::MB_TYPE_INTEGER, 477 part_tag, moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &rank);MBERRNM(merr); 478 rank=pcomm->rank(); 479 merr = mbiface->create_meshset(moab::MESHSET_SET, partset);MBERRNM(merr); 480 merr = mbiface->add_entities(partset, localvtxs);MBERRNM(merr); 481 merr = mbiface->add_entities(partset, localelms);MBERRNM(merr); 482 merr = mbiface->tag_set_data(part_tag, &partset, 1, &rank);MBERRNM(merr); 483 484 /* Everything is set up, now just do a tag exchange to update tags 485 on all of the ghost vertexes */ 486 merr = mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr); 487 merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, ownedelms);MBERRNM(merr); 488 489 merr = pcomm->exchange_tags(id_tag,ownedvtx);MBERRV(mbiface,merr); 490 merr = pcomm->exchange_tags(id_tag,ownedelms);MBERRV(mbiface,merr); 491 PetscFunctionReturn(0); 492 } 493 494 495 PetscErrorCode DMMoab_GetReadOptions_Private(PetscBool by_rank, PetscInt numproc, PetscInt dim, MoabReadMode mode, PetscInt dbglevel, const char* dm_opts, const char* extra_opts, const char** read_opts) 496 { 497 char *ropts; 498 char ropts_par[PETSC_MAX_PATH_LEN]; 499 char ropts_dbg[PETSC_MAX_PATH_LEN]; 500 PetscErrorCode ierr; 501 502 PetscFunctionBegin; 503 ierr = PetscMalloc1(PETSC_MAX_PATH_LEN,&ropts);CHKERRQ(ierr); 504 ierr = PetscMemzero(&ropts_par,PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 505 ierr = PetscMemzero(&ropts_dbg,PETSC_MAX_PATH_LEN);CHKERRQ(ierr); 506 507 /* do parallel read unless using only one processor */ 508 if (numproc > 1) { 509 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); 510 } 511 512 if (dbglevel) { 513 if (numproc>1) { 514 ierr = PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "%sCPUTIME;DEBUG_IO=%d;DEBUG_PIO=%d;",dbglevel,dbglevel);CHKERRQ(ierr); 515 } 516 else { 517 ierr = PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "%sCPUTIME;DEBUG_IO=%d;",dbglevel);CHKERRQ(ierr); 518 } 519 } 520 521 ierr = PetscSNPrintf(ropts, PETSC_MAX_PATH_LEN, "%s%s%s%s",ropts_par,ropts_dbg,(extra_opts?extra_opts:""),(dm_opts?dm_opts:""));CHKERRQ(ierr); 522 *read_opts = ropts; 523 PetscFunctionReturn(0); 524 } 525 526 527 /*@ 528 DMMoabLoadFromFile - Creates a DM object by loading the mesh from a user specified file. 529 530 Collective on MPI_Comm 531 532 Input Parameters: 533 + comm - The communicator for the DM object 534 . dim - The spatial dimension 535 . filename - The name of the mesh file to be loaded 536 . usrreadopts - The options string to read a MOAB mesh. 537 Reference (Parallel Mesh Initialization: http://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html#fivetwo) 538 539 Output Parameter: 540 . dm - The DM object 541 542 Level: beginner 543 544 .keywords: DM, create 545 546 .seealso: DMSetType(), DMCreate(), DMMoabCreateBoxMesh() 547 @*/ 548 PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm,PetscInt dim,const char* filename, const char* usrreadopts, DM *dm) 549 { 550 moab::ErrorCode merr; 551 PetscInt nprocs; 552 DM_Moab *dmmoab; 553 moab::Interface *mbiface; 554 moab::ParallelComm *pcomm; 555 moab::Range verts,elems; 556 const char *readopts; 557 PetscErrorCode ierr; 558 559 PetscFunctionBegin; 560 PetscValidPointer(dm,5); 561 562 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 563 ierr = DMMoabCreateMoab(comm, NULL, NULL, NULL, NULL, dm);CHKERRQ(ierr); 564 565 /* get all the necessary handles from the private DM object */ 566 dmmoab = (DM_Moab*)(*dm)->data; 567 mbiface = dmmoab->mbiface; 568 pcomm = dmmoab->pcomm; 569 nprocs = pcomm->size(); 570 /* TODO: Decipher dimension based on the loaded mesh instead of getting from user */ 571 dmmoab->dim = dim; 572 573 /* create a file set to associate all entities in current mesh */ 574 merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, dmmoab->fileset);MBERR("Creating file set failed", merr); 575 576 /* add mesh loading options specific to the DM */ 577 ierr = DMMoab_GetReadOptions_Private(dmmoab->partition_by_rank, nprocs, dim, dmmoab->read_mode, 578 dmmoab->rw_dbglevel, dmmoab->extra_read_options, usrreadopts, &readopts);CHKERRQ(ierr); 579 580 PetscInfo2(*dm, "Reading file %s with options: %s\n",filename,readopts); 581 582 /* Load the mesh from a file. */ 583 merr = mbiface->load_file(filename, &dmmoab->fileset, readopts);MBERRVM(mbiface,"Reading MOAB file failed.", merr); 584 585 /* Reassign global IDs on all entities. */ 586 merr = pcomm->assign_global_ids(dmmoab->fileset,dim,1,true,true,true);MBERRNM(merr); 587 588 /* load the local vertices */ 589 merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);MBERRNM(merr); 590 /* load the local elements */ 591 merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);MBERRNM(merr); 592 593 /* Everything is set up, now just do a tag exchange to update tags 594 on all of the ghost vertexes */ 595 merr = pcomm->exchange_tags(dmmoab->ltog_tag,verts);MBERRV(mbiface,merr); 596 merr = pcomm->exchange_tags(dmmoab->ltog_tag,elems);MBERRV(mbiface,merr); 597 598 merr = pcomm->exchange_ghost_cells(dim,0,1,0,true,true,&dmmoab->fileset);MBERRV(mbiface,merr); 599 600 merr = pcomm->collective_sync_partition();MBERR("Collective sync failed", merr); 601 602 PetscInfo3(*dm, "MOAB file '%s' was successfully loaded. Found %D vertices and %D elements.\n", filename, verts.size(), elems.size()); 603 ierr = PetscFree(readopts);CHKERRQ(ierr); 604 PetscFunctionReturn(0); 605 } 606 607