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 50 #undef __FUNCT__ 51 #define __FUNCT__ "DMMoabCreateBoxMesh" 52 PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscInt nele, PetscInt nghost, DM *dm) 53 { 54 PetscErrorCode ierr; 55 moab::ErrorCode merr; 56 PetscInt nprocs; 57 DM_Moab *dmmoab; 58 moab::Interface *mbiface; 59 moab::ParallelComm *pcomm; 60 moab::ReadUtilIface* readMeshIface; 61 62 moab::Tag id_tag=PETSC_NULL; 63 moab::Range ownedvtx,ownedelms; 64 moab::EntityHandle vfirst,efirst; 65 std::vector<double*> vcoords; 66 moab::EntityHandle *connectivity = 0; 67 std::vector<int> subent_conn; 68 moab::EntityType etype; 69 PetscInt ise[6]; 70 PetscReal xse[6]; 71 72 const PetscInt npts=nele+1; /* Number of points in every dimension */ 73 PetscInt vpere,locnele,locnpts; /* Number of verts/element, vertices, elements owned by this process */ 74 75 PetscFunctionBegin; 76 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 77 ierr = DMMoabCreateMoab(comm, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, dm);CHKERRQ(ierr); 78 79 /* get all the necessary handles from the private DM object */ 80 dmmoab = (DM_Moab*)(*dm)->data; 81 mbiface = dmmoab->mbiface; 82 pcomm = dmmoab->pcomm; 83 id_tag = dmmoab->ltog_tag; 84 nprocs = pcomm->size(); 85 86 /* do some error checking */ 87 if(pow(npts,dim) < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be >= 2"); 88 if(nprocs > pow(nele,dim)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of processors must be less than or equal to number of elements"); 89 90 /* No errors yet; proceed with building the mesh */ 91 merr = mbiface->query_interface(readMeshIface);MBERRNM(merr); 92 93 ierr = PetscMemzero(ise,sizeof(PetscInt)*6);CHKERRQ(ierr); 94 95 /* call the collective routine that computes the domain bounds for a structured mesh using MOAB */ 96 ierr = DMMoabComputeDomainBounds_Private(pcomm, dim, nele, ise);CHKERRQ(ierr); 97 98 /* set some variables based on dimension */ 99 switch(dim) { 100 case 1: 101 vpere = 2; 102 locnele = (ise[1]-ise[0]); 103 locnpts = (ise[1]-ise[0]+1); 104 etype = moab::MBEDGE; 105 break; 106 case 2: 107 vpere = 4; 108 locnele = (ise[1]-ise[0])*(ise[3]-ise[2]); 109 locnpts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1); 110 etype = moab::MBQUAD; 111 break; 112 case 3: 113 vpere = 8; 114 locnele = (ise[1]-ise[0])*(ise[3]-ise[2])*(ise[5]-ise[4]); 115 locnpts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1)*(ise[5]-ise[4]+1); 116 etype = moab::MBHEX; 117 break; 118 } 119 120 /* we have a domain of size [1,1,1] - now compute local co-ordinate box */ 121 ierr = PetscMemzero(xse,sizeof(PetscReal)*6);CHKERRQ(ierr); 122 for(int i=0; i<6; ++i) { 123 xse[i]=(PetscReal)ise[i]/nele; 124 } 125 126 /* Create vertexes and set the coodinate of each vertex */ 127 const int sequence_size = (locnele + vpere) + 1; 128 merr = readMeshIface->get_node_coords(3,locnpts,0,vfirst,vcoords,sequence_size);MBERRNM(merr); 129 130 /* Compute the co-ordinates of vertices and global IDs */ 131 std::vector<int> vgid(locnpts); 132 int vcount=0; 133 double hxyz=1.0/nele; 134 for (int k = ise[4]; k <= ise[5]; k++) { 135 for (int j = ise[2]; j <= ise[3]; j++) { 136 for (int i = ise[0]; i <= ise[1]; i++, vcount++) { 137 vcoords[0][vcount] = i*hxyz; 138 vcoords[1][vcount] = j*hxyz; 139 vcoords[2][vcount] = k*hxyz; 140 vgid[vcount] = (k*nele+j)*nele+i; 141 } 142 } 143 } 144 145 merr = mbiface->get_entities_by_type(0,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr); 146 merr = mbiface->add_entities (dmmoab->fileset, ownedvtx);MBERRNM(merr); 147 148 /* The global ID tag is applied to each owned 149 vertex. It acts as an global identifier which MOAB uses to 150 assemble the individual pieces of the mesh: 151 Set the global ID indices */ 152 merr = mbiface->tag_set_data(id_tag,ownedvtx,vgid.data());MBERRNM(merr); 153 154 /* Create elements between mesh points using the ReadUtilInterface 155 get the reference to element connectivities for all local elements from the ReadUtilInterface */ 156 merr = readMeshIface->get_element_connect (locnele,vpere,etype,1,efirst,connectivity);MBERRNM(merr); 157 158 /* offset appropriately so that only local ID and not global ID numbers are set for connectivity array */ 159 vfirst-=vgid[0]; 160 161 /* 3. Loop over elements in 3 nested loops over i, j, k; for each (i,j,k): 162 and then set the connectivity for each element appropriately */ 163 int ecount=0; 164 subent_conn.resize(vpere); 165 for (int k = ise[4]; k < std::max(ise[5],1); k++) { 166 for (int j = ise[2]; j < std::max(ise[3],1); j++) { 167 for (int i = ise[0]; i < std::max(ise[1],1); i++,ecount++) { 168 const int offset = ecount*vpere; 169 moab::CN::SubEntityVertexIndices(etype, dim, 0, subent_conn.data()); 170 171 switch(dim) { 172 case 1: 173 connectivity[offset+subent_conn[0]] = vfirst+i; 174 connectivity[offset+subent_conn[1]] = vfirst+(i+1); 175 PetscPrintf(PETSC_COMM_WORLD, "ELEMENT[%D,%D,%D]: CONNECTIVITY = %D, %D\n", i,j,k,connectivity[offset+subent_conn[0]], connectivity[offset+subent_conn[1]]); 176 break; 177 case 2: 178 connectivity[offset+subent_conn[0]] = vfirst+i+j*(nele+1); 179 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+j*(nele+1); 180 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(j+1)*(nele+1); 181 connectivity[offset+subent_conn[3]] = vfirst+i+(j+1)*(nele+1); 182 PetscPrintf(PETSC_COMM_WORLD, "ELEMENT[%D,%D,%D]: CONNECTIVITY = %D, %D, %D, %D\n", i,j,k,connectivity[offset+subent_conn[0]], connectivity[offset+subent_conn[1]], connectivity[offset+subent_conn[2]], connectivity[offset+subent_conn[3]]); 183 break; 184 case 3: 185 connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k); 186 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*k); 187 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k); 188 connectivity[offset+subent_conn[3]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*k); 189 connectivity[offset+subent_conn[4]] = vfirst+i+(nele+1)*(j+(nele+1)*(k+1)); 190 connectivity[offset+subent_conn[5]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1)); 191 connectivity[offset+subent_conn[6]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*(k+1)); 192 connectivity[offset+subent_conn[7]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1)); 193 break; 194 } 195 } 196 } 197 } 198 merr = readMeshIface->update_adjacencies(efirst,locnele,vpere,connectivity);MBERRNM(merr); 199 200 /* 2. Get the vertices and hexes from moab and check their numbers against I*J*K and (I-1)*(J-1)*(K-1), resp. 201 first '0' specifies "root set", or entire MOAB instance, second the entity dimension being requested */ 202 merr = mbiface->get_entities_by_dimension(0, dim, ownedelms);MBERRNM(merr); 203 merr = mbiface->add_entities (dmmoab->fileset, ownedelms);MBERRNM(merr); 204 205 // merr = mbiface->unite_meshset(dmmoab->fileset, 0);MBERRV(mbiface,merr); 206 207 if (locnele != (int) ownedelms.size()) 208 SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of elements! (%D!=%D)",locnele,ownedelms.size()); 209 else if(locnpts != (int) ownedvtx.size()) 210 SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of vertices! (%D!=%D)",locnpts,ownedvtx.size()); 211 else 212 PetscInfo2(*dm, "Created %D elements and %D vertices.\n", ownedelms.size(), ownedvtx.size()); 213 214 /* check the handles */ 215 merr = pcomm->check_all_shared_handles();MBERRV(mbiface,merr); 216 217 /* resolve the shared entities by exchanging information to adjacent processors */ 218 merr = mbiface->get_entities_by_type(dmmoab->fileset,etype,ownedelms,true);MBERRNM(merr); 219 merr = pcomm->resolve_shared_ents(dmmoab->fileset,ownedelms,dim,dim-1,&id_tag);MBERRV(mbiface,merr); 220 221 /* Reassign global IDs on all entities. */ 222 merr = pcomm->assign_global_ids(dmmoab->fileset,dim,1,false,true);MBERRNM(merr); 223 merr = pcomm->exchange_ghost_cells(dim,0,nghost,0,true);MBERRV(mbiface,merr); 224 225 /* Everything is set up, now just do a tag exchange to update tags 226 on all of the ghost vertexes */ 227 merr = pcomm->exchange_tags(id_tag,ownedvtx);MBERRV(mbiface,merr); 228 merr = pcomm->exchange_tags(id_tag,ownedelms);MBERRV(mbiface,merr); 229 230 // ierr = DMMoabSetLocalVertices(*dm, &ownedvtx);CHKERRQ(ierr); 231 // ierr = DMMoabSetLocalElements(*dm, &ownedelms);CHKERRQ(ierr); 232 233 merr = mbiface->set_dimension(dim);MBERRNM(merr); 234 235 std::stringstream sstr; 236 sstr << "test_" << pcomm->rank() << ".vtk"; 237 mbiface->write_mesh(sstr.str().c_str()); 238 PetscFunctionReturn(0); 239 } 240 241 242 #undef __FUNCT__ 243 #define __FUNCT__ "DMMoab_GetReadOptions_Private" 244 PetscErrorCode DMMoab_GetReadOptions_Private(PetscBool by_rank, PetscInt numproc, PetscInt dim, MoabReadMode mode, PetscInt dbglevel, const char* extra_opts, const char* read_opts) 245 { 246 std::ostringstream str; 247 248 PetscFunctionBegin; 249 // do parallel read unless only one processor 250 if (numproc > 1) { 251 str << "PARALLEL=" << mode << ";PARTITION=PARALLEL_PARTITION;"; 252 str << "PARTITION_DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=" << dim << ".0.1;"; 253 if (by_rank) 254 str << "PARTITION_BY_RANK;"; 255 } 256 257 if (extra_opts) 258 str << extra_opts << ";"; 259 260 if (dbglevel) 261 str << "DEBUG_IO=" << dbglevel << ";CPUTIME;"; 262 263 read_opts = str.str().c_str(); 264 PetscFunctionReturn(0); 265 } 266 267 268 #undef __FUNCT__ 269 #define __FUNCT__ "DMMoabLoadFromFile" 270 PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm,PetscInt dim,const char* filename, const char* usrreadopts, DM *dm) 271 { 272 moab::ErrorCode merr; 273 PetscInt nprocs; 274 DM_Moab *dmmoab; 275 moab::Interface *mbiface; 276 moab::ParallelComm *pcomm; 277 moab::Range verts,elems; 278 const char* readopts=0; 279 PetscErrorCode ierr; 280 281 PetscFunctionBegin; 282 PetscValidPointer(dm,5); 283 284 /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 285 ierr = DMMoabCreateMoab(comm, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, dm);CHKERRQ(ierr); 286 287 /* get all the necessary handles from the private DM object */ 288 dmmoab = (DM_Moab*)(*dm)->data; 289 mbiface = dmmoab->mbiface; 290 pcomm = dmmoab->pcomm; 291 nprocs = pcomm->size(); 292 293 /* add mesh loading options specific to the DM */ 294 ierr = DMMoab_GetReadOptions_Private(PETSC_TRUE, nprocs, dim, MOAB_PARROPTS_READ_PART, 0, usrreadopts, readopts);CHKERRQ(ierr); 295 296 /* Load the mesh from a file. */ 297 merr = mbiface->load_file(filename, &dmmoab->fileset, ""/*readopts*/);MBERRVM(mbiface,"Reading MOAB file failed.", merr); 298 299 merr = pcomm->collective_sync_partition();MBERR("Collective sync failed", merr); 300 301 /* load the local vertices */ 302 merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);MBERRNM(merr); 303 ierr = DMMoabSetLocalVertices(*dm, &verts);CHKERRQ(ierr); 304 305 /* load the local elements */ 306 merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);MBERRNM(merr); 307 // ierr = DMMoabSetLocalElements(*dm, &elems);CHKERRQ(ierr); 308 309 merr = mbiface->set_dimension(dim);MBERRNM(merr); 310 311 PetscInfo3(*dm, "MOAB file '%s' was successfully loaded. Found %D vertices and %D elements.\n", filename, verts.size(), elems.size()); 312 PetscFunctionReturn(0); 313 } 314 315