1 #include <petsc-private/dmmbimpl.h> /*I "petscdm.h" I*/ 2 #include <petsc-private/vecimpl.h> /*I "petscdm.h" I*/ 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 #undef __FUNCT__ 12 #define __FUNCT__ "DMMoabComputeDomainBounds_Private" 13 PetscErrorCode DMMoabComputeDomainBounds_Private(moab::ParallelComm* pcomm, PetscInt dim, PetscInt neleglob, PetscInt *ise) 14 { 15 PetscInt size,rank; 16 PetscInt fraction,remainder; 17 PetscInt starts[3],sizes[3]; 18 19 PetscFunctionBegin; 20 if(dim<1 && dim>3) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The problem dimension is invalid: %D",dim); 21 22 size=pcomm->size(); 23 rank=pcomm->rank(); 24 fraction=neleglob/size; /* partition only by the largest dimension */ 25 remainder=neleglob%size; /* remainder after partition which gets evenly distributed by round-robin */ 26 27 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); 28 29 starts[0]=starts[1]=starts[2]=0; /* default dimensional element = 1 */ 30 sizes[0]=sizes[1]=sizes[2]=neleglob; /* default dimensional element = 1 */ 31 32 if(rank < remainder) { 33 /* This process gets "fraction+1" elements */ 34 sizes[dim-1] = (fraction + 1); 35 starts[dim-1] = rank * (fraction+1); 36 } else { 37 /* This process gets "fraction" elements */ 38 sizes[dim-1] = fraction; 39 starts[dim-1] = (remainder*(fraction+1) + fraction*(rank-remainder)); 40 } 41 42 for(int i=dim-1; i>=0; --i) { 43 ise[2*i]=starts[i];ise[2*i+1]=starts[i]+sizes[i]; 44 } 45 46 PetscFunctionReturn(0); 47 } 48 49 static void set_structured_coordinates(PetscInt i, PetscInt j, PetscInt k, PetscReal hxyz, PetscInt vcount, std::vector<double*>& vcoords) 50 { 51 vcoords[0][vcount] = i*hxyz; 52 vcoords[1][vcount] = j*hxyz; 53 vcoords[2][vcount] = k*hxyz; 54 } 55 56 static void set_element_connectivity(PetscInt dim, moab::EntityType etype, PetscInt offset, PetscInt nele, PetscInt i, PetscInt j, PetscInt k, PetscInt vfirst, moab::EntityHandle *connectivity) 57 { 58 std::vector<int> subent_conn(pow(2,dim)); /* only linear edge, quad, hex supported now */ 59 60 moab::CN::SubEntityVertexIndices(etype, dim, 0, subent_conn.data()); 61 62 switch(dim) { 63 case 1: 64 connectivity[offset+subent_conn[0]] = vfirst+i; 65 connectivity[offset+subent_conn[1]] = vfirst+(i+1); 66 break; 67 case 2: 68 connectivity[offset+subent_conn[0]] = vfirst+i+j*(nele+1); 69 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+j*(nele+1); 70 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(j+1)*(nele+1); 71 connectivity[offset+subent_conn[3]] = vfirst+i+(j+1)*(nele+1); 72 break; 73 case 3: 74 default: 75 connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k); 76 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*k); 77 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k); 78 connectivity[offset+subent_conn[3]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*k); 79 connectivity[offset+subent_conn[4]] = vfirst+i+(nele+1)*(j+(nele+1)*(k+1)); 80 connectivity[offset+subent_conn[5]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1)); 81 connectivity[offset+subent_conn[6]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*(k+1)); 82 connectivity[offset+subent_conn[7]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1)); 83 break; 84 } 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "DMMoabCreateBoxMesh" 89 PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscInt nele, PetscInt nghost, DM *dm) 90 { 91 PetscErrorCode ierr; 92 moab::ErrorCode merr; 93 PetscInt i,j,k,n,nprocs,rank; 94 DM_Moab *dmmoab; 95 moab::Interface *mbiface; 96 moab::ParallelComm *pcomm; 97 moab::ReadUtilIface* readMeshIface; 98 99 moab::Tag id_tag=PETSC_NULL; 100 moab::Range ownedvtx,ownedelms; 101 moab::EntityHandle vfirst,efirst; 102 std::vector<double*> vcoords; 103 moab::EntityHandle *connectivity = 0; 104 moab::EntityType etype; 105 PetscInt ise[6]; 106 PetscReal xse[6]; 107 108 const PetscInt npts=nele+1; /* Number of points in every dimension */ 109 PetscInt vpere,locnele,locnpts,ghnele,ghnpts; /* Number of verts/element, vertices, elements owned by this process */ 110 111 PetscFunctionBegin; 112 if(dim < 1 || dim > 3) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Invalid dimension argument for mesh: dim=[1,3].\n"); 113 114 ierr = MPI_Comm_size(comm, &nprocs);CHKERRQ(ierr); 115 /* total number of vertices in all dimensions */ 116 n=pow(npts,dim); 117 118 /* do some error checking */ 119 if(n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be >= 2.\n"); 120 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"); 121 if(nghost < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of ghost layers cannot be negative.\n"); 122 123 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 124 ierr = DMMoabCreateMoab(comm, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, dm);CHKERRQ(ierr); 125 126 /* get all the necessary handles from the private DM object */ 127 dmmoab = (DM_Moab*)(*dm)->data; 128 mbiface = dmmoab->mbiface; 129 pcomm = dmmoab->pcomm; 130 id_tag = dmmoab->ltog_tag; 131 nprocs = pcomm->size(); 132 rank = pcomm->rank(); 133 dmmoab->dim = dim; 134 135 /* No errors yet; proceed with building the mesh */ 136 merr = mbiface->query_interface(readMeshIface);MBERRNM(merr); 137 138 ierr = PetscMemzero(ise,sizeof(PetscInt)*6);CHKERRQ(ierr); 139 140 /* call the collective routine that computes the domain bounds for a structured mesh using MOAB */ 141 ierr = DMMoabComputeDomainBounds_Private(pcomm, dim, nele, ise);CHKERRQ(ierr); 142 143 /* set some variables based on dimension */ 144 switch(dim) { 145 case 1: 146 vpere = 2; 147 locnele = (ise[1]-ise[0]); 148 locnpts = (ise[1]-ise[0]+1); 149 ghnele = (nghost > 0 ? (ise[0] > nghost ? 1 : 0) + (ise[1] < nele - nghost ? 1 : 0) : 0 ); 150 ghnpts = (nghost > 0 ? (ise[0] > 0 ? 1 : 0) + (ise[1] < nele ? 1 : 0) : 0); 151 etype = moab::MBEDGE; 152 break; 153 case 2: 154 vpere = 4; 155 locnele = (ise[1]-ise[0])*(ise[3]-ise[2]); 156 locnpts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1); 157 ghnele = (nghost > 0 ? (ise[2] > 0 ? nele : 0) + (ise[3] < nele ? nele : 0) : 0); 158 ghnpts = (nghost > 0 ? (ise[2] > 0 ? npts : 0) + (ise[3] < nele ? npts : 0) : 0); 159 etype = moab::MBQUAD; 160 break; 161 case 3: 162 vpere = 8; 163 locnele = (ise[1]-ise[0])*(ise[3]-ise[2])*(ise[5]-ise[4]); 164 locnpts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1)*(ise[5]-ise[4]+1); 165 ghnele = (nghost > 0 ? (ise[4] > 0 ? nele*nele : 0) + (ise[5] < nele ? nele*nele : 0) : 0); 166 ghnpts = (nghost > 0 ? (ise[4] > 0 ? npts*npts : 0) + (ise[5] < nele ? npts*npts : 0) : 0); 167 etype = moab::MBHEX; 168 break; 169 } 170 171 /* we have a domain of size [1,1,1] - now compute local co-ordinate box */ 172 ierr = PetscMemzero(xse,sizeof(PetscReal)*6);CHKERRQ(ierr); 173 for(i=0; i<6; ++i) { 174 xse[i]=(PetscReal)ise[i]/nele; 175 } 176 177 /* Create vertexes and set the coodinate of each vertex */ 178 merr = readMeshIface->get_node_coords(3,locnpts+ghnpts,0,vfirst,vcoords,n);MBERRNM(merr); 179 180 /* Compute the co-ordinates of vertices and global IDs */ 181 std::vector<int> vgid(locnpts+ghnpts); 182 int vcount=0; 183 double hxyz=1.0/nele; 184 185 /* create all the owned vertices */ 186 for (k = ise[4]; k <= ise[5]; k++) { 187 for (j = ise[2]; j <= ise[3]; j++) { 188 for (i = ise[0]; i <= ise[1]; i++, vcount++) { 189 set_structured_coordinates(i,j,k,hxyz,vcount,vcoords); 190 vgid[vcount] = (k*npts+j)*npts+i+1; 191 } 192 } 193 } 194 195 /* create ghosted vertices requested by user - below the current plane */ 196 if (ise[2*dim-2] > 0) { 197 for (k = (dim==3?ise[4]-nghost:ise[4]); k <= (dim==3?ise[4]-1:ise[5]); k++) { 198 for (j = (dim==2?ise[2]-nghost:ise[2]); j <= (dim==2?ise[2]-1:ise[3]); j++) { 199 for (i = (dim>1?ise[0]:ise[0]-nghost); i <= (dim>1?ise[1]:ise[0]-1); i++, vcount++) { 200 set_structured_coordinates(i,j,k,hxyz,vcount,vcoords); 201 vgid[vcount] = (k*npts+j)*npts+i+1; 202 } 203 } 204 } 205 } 206 207 /* create ghosted vertices requested by user - above the current plane */ 208 if (ise[2*dim-1] < nele) { 209 for (k = (dim==3?ise[5]+1:ise[4]); k <= (dim==3?ise[5]+nghost:ise[5]); k++) { 210 for (j = (dim==2?ise[3]+1:ise[2]); j <= (dim==2?ise[3]+nghost:ise[3]); j++) { 211 for (i = (dim>1?ise[0]:ise[1]+1); i <= (dim>1?ise[1]:ise[1]+nghost); i++, vcount++) { 212 set_structured_coordinates(i,j,k,hxyz,vcount,vcoords); 213 vgid[vcount] = (k*npts+j)*npts+i+1; 214 } 215 } 216 } 217 } 218 219 if (locnpts+ghnpts != vcount) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of vertices! (%D!=%D)",locnpts+ghnpts,vcount); 220 221 merr = mbiface->get_entities_by_type(0,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr); 222 merr = mbiface->add_entities (dmmoab->fileset, ownedvtx);MBERRNM(merr); 223 224 /* The global ID tag is applied to each owned 225 vertex. It acts as an global identifier which MOAB uses to 226 assemble the individual pieces of the mesh: 227 Set the global ID indices */ 228 merr = mbiface->tag_set_data(id_tag,ownedvtx,vgid.data());MBERRNM(merr); 229 230 /* Create elements between mesh points using the ReadUtilInterface 231 get the reference to element connectivities for all local elements from the ReadUtilInterface */ 232 merr = readMeshIface->get_element_connect (locnele+ghnele,vpere,etype,1,efirst,connectivity);MBERRNM(merr); 233 234 /* offset appropriately so that only local ID and not global ID numbers are set for connectivity array */ 235 // PetscPrintf(PETSC_COMM_SELF, "[%D] first local handle %D\n", rank, vfirst); 236 vfirst-=vgid[0]-1; 237 238 /* 3. Loop over elements in 3 nested loops over i, j, k; for each (i,j,k): 239 and then set the connectivity for each element appropriately */ 240 int ecount=0; 241 242 /* create ghosted elements requested by user - below the current plane */ 243 if (ise[2*dim-2] >= nghost) { 244 for (k = (dim==3?ise[4]-nghost:ise[4]); k < (dim==3?ise[4]:std::max(ise[5],1)); k++) { 245 for (j = (dim==2?ise[2]-nghost:ise[2]); j < (dim==2?ise[2]:std::max(ise[3],1)); j++) { 246 for (i = (dim>1?ise[0]:ise[0]-nghost); i < (dim>1?std::max(ise[1],1):ise[0]); i++, ecount++) { 247 set_element_connectivity(dim, etype, ecount*vpere, nele, i, j, k, vfirst, connectivity); 248 } 249 } 250 } 251 } 252 253 /* create owned elements requested by user */ 254 for (k = ise[4]; k < std::max(ise[5],1); k++) { 255 for (j = ise[2]; j < std::max(ise[3],1); j++) { 256 for (i = ise[0]; i < std::max(ise[1],1); i++,ecount++) { 257 set_element_connectivity(dim, etype, ecount*vpere, nele, i, j, k, vfirst, connectivity); 258 } 259 } 260 } 261 262 /* create ghosted elements requested by user - above the current plane */ 263 if (ise[2*dim-1] <= nele-nghost) { 264 for (k = (dim==3?ise[5]:ise[4]); k < (dim==3?ise[5]+nghost:std::max(ise[5],1)); k++) { 265 for (j = (dim==2?ise[3]:ise[2]); j < (dim==2?ise[3]+nghost:std::max(ise[3],1)); j++) { 266 for (i = (dim>1?ise[0]:ise[1]); i < (dim>1?std::max(ise[1],1):ise[1]+nghost); i++, ecount++) { 267 set_element_connectivity(dim, etype, ecount*vpere, nele, i, j, k, vfirst, connectivity); 268 } 269 } 270 } 271 } 272 273 merr = readMeshIface->update_adjacencies(efirst,locnele+ghnele,vpere,connectivity);MBERRNM(merr); 274 275 /* 2. Get the vertices and hexes from moab and check their numbers against I*J*K and (I-1)*(J-1)*(K-1), resp. 276 first '0' specifies "root set", or entire MOAB instance, second the entity dimension being requested */ 277 merr = mbiface->get_entities_by_dimension(0, dim, ownedelms);MBERRNM(merr); 278 merr = mbiface->add_entities (dmmoab->fileset, ownedelms);MBERRNM(merr); 279 280 if (locnele+ghnele != (int) ownedelms.size()) 281 SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of elements! (%D!=%D)",locnele+ghnele,ownedelms.size()); 282 else if(locnpts+ghnpts != (int) ownedvtx.size()) 283 SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of vertices! (%D!=%D)",locnpts+ghnpts,ownedvtx.size()); 284 else 285 PetscInfo2(*dm, "Created %D elements and %D vertices.\n", ownedelms.size(), ownedvtx.size()); 286 287 /* check the handles */ 288 merr = pcomm->check_all_shared_handles();MBERRV(mbiface,merr); 289 290 #if 0 291 std::stringstream sstr; 292 sstr << "test_" << pcomm->rank() << ".vtk"; 293 mbiface->write_mesh(sstr.str().c_str()); 294 #endif 295 296 /* resolve the shared entities by exchanging information to adjacent processors */ 297 merr = mbiface->get_entities_by_type(dmmoab->fileset,etype,ownedelms,true);MBERRNM(merr); 298 merr = pcomm->resolve_shared_ents(dmmoab->fileset,ownedelms,dim,dim-1,&id_tag);MBERRV(mbiface,merr); 299 300 /* Reassign global IDs on all entities. */ 301 merr = pcomm->assign_global_ids(dmmoab->fileset,dim,1,true,true,true);MBERRNM(merr); 302 303 merr = pcomm->exchange_ghost_cells(dim,0,1/*nghost*/,dim,true,false,&dmmoab->fileset);MBERRV(mbiface,merr); 304 305 /* Everything is set up, now just do a tag exchange to update tags 306 on all of the ghost vertexes */ 307 merr = mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr); 308 merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, ownedelms);MBERRNM(merr); 309 310 merr = pcomm->exchange_tags(id_tag,ownedvtx);MBERRV(mbiface,merr); 311 merr = pcomm->exchange_tags(id_tag,ownedelms);MBERRV(mbiface,merr); 312 PetscFunctionReturn(0); 313 } 314 315 316 #undef __FUNCT__ 317 #define __FUNCT__ "DMMoab_GetReadOptions_Private" 318 PetscErrorCode DMMoab_GetReadOptions_Private(PetscBool by_rank, PetscInt numproc, PetscInt dim, MoabReadMode mode, PetscInt dbglevel, const char* extra_opts, const char** read_opts) 319 { 320 std::ostringstream str; 321 322 PetscFunctionBegin; 323 /* do parallel read unless using only one processor */ 324 if (numproc > 1) { 325 str << "PARALLEL=" << mode << ";PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;"; 326 str << "PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=" << dim << ".0.1;"; 327 if (by_rank) 328 str << "PARTITION_BY_RANK;"; 329 } 330 331 if (dbglevel) { 332 str << "CPUTIME;DEBUG_IO=" << dbglevel << ";"; 333 if (numproc>1) 334 str << "DEBUG_PIO=" << dbglevel << ";"; 335 } 336 337 if (extra_opts) 338 str << extra_opts; 339 340 *read_opts = str.str().c_str(); 341 PetscFunctionReturn(0); 342 } 343 344 345 #undef __FUNCT__ 346 #define __FUNCT__ "DMMoabLoadFromFile" 347 PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm,PetscInt dim,const char* filename, const char* usrreadopts, DM *dm) 348 { 349 moab::ErrorCode merr; 350 PetscInt nprocs,dbglevel=0; 351 PetscBool part_by_rank=PETSC_FALSE; 352 DM_Moab *dmmoab; 353 moab::Interface *mbiface; 354 moab::ParallelComm *pcomm; 355 moab::Range verts,elems; 356 const char *readopts; 357 PetscErrorCode ierr; 358 359 PetscFunctionBegin; 360 PetscValidPointer(dm,5); 361 362 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 363 ierr = DMMoabCreateMoab(comm, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, dm);CHKERRQ(ierr); 364 365 /* get all the necessary handles from the private DM object */ 366 dmmoab = (DM_Moab*)(*dm)->data; 367 mbiface = dmmoab->mbiface; 368 pcomm = dmmoab->pcomm; 369 nprocs = pcomm->size(); 370 dmmoab->dim = dim; 371 372 /* TODO: Use command-line options to control by_rank, verbosity, MoabReadMode and extra options */ 373 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for reading/writing MOAB based meshes from file", "DMMoab"); 374 ierr = PetscOptionsInt("-dmmb_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "dmmbutil.cxx", dbglevel, &dbglevel, NULL);CHKERRQ(ierr); 375 ierr = PetscOptionsBool("-dmmb_part_by_rank", "Use partition by rank when reading MOAB meshes from file", "dmmbutil.cxx", part_by_rank, &part_by_rank, NULL);CHKERRQ(ierr); 376 ierr = PetscOptionsEnd(); 377 378 /* add mesh loading options specific to the DM */ 379 ierr = DMMoab_GetReadOptions_Private(part_by_rank, nprocs, dim, MOAB_PARROPTS_READ_PART, dbglevel, usrreadopts, &readopts);CHKERRQ(ierr); 380 381 PetscInfo2(*dm, "Reading file %s with options: %s\n",filename,readopts); 382 383 /* Load the mesh from a file. */ 384 merr = mbiface->load_file(filename, &dmmoab->fileset, readopts);MBERRVM(mbiface,"Reading MOAB file failed.", merr); 385 386 /* Reassign global IDs on all entities. */ 387 merr = pcomm->assign_global_ids(dmmoab->fileset,dim,1,true,true,true);MBERRNM(merr); 388 389 /* load the local vertices */ 390 merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);MBERRNM(merr); 391 /* load the local elements */ 392 merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);MBERRNM(merr); 393 394 /* Everything is set up, now just do a tag exchange to update tags 395 on all of the ghost vertexes */ 396 merr = pcomm->exchange_tags(dmmoab->ltog_tag,verts);MBERRV(mbiface,merr); 397 merr = pcomm->exchange_tags(dmmoab->ltog_tag,elems);MBERRV(mbiface,merr); 398 399 merr = pcomm->exchange_ghost_cells(dim,0,1,0,true,true,&dmmoab->fileset);MBERRV(mbiface,merr); 400 401 merr = pcomm->collective_sync_partition();MBERR("Collective sync failed", merr); 402 403 PetscInfo3(*dm, "MOAB file '%s' was successfully loaded. Found %D vertices and %D elements.\n", filename, verts.size(), elems.size()); 404 405 #if 1 406 std::stringstream sstr; 407 sstr << "test_" << pcomm->rank() << ".vtk"; 408 mbiface->write_mesh(sstr.str().c_str()); 409 #endif 410 PetscFunctionReturn(0); 411 } 412 413