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 neleadim; 18 PetscInt starts[3],sizes[3]; 19 20 PetscFunctionBegin; 21 if(dim<1 && dim>3) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The problem dimension is invalid: %D",dim); 22 23 size=pcomm->size(); 24 rank=pcomm->rank(); 25 neleadim=(dim==3?neleglob*neleglob:(dim==2?neleglob:1)); 26 fraction=neleglob/size; /* partition only by the largest dimension */ 27 remainder=neleglob%size; /* remainder after partition which gets evenly distributed by round-robin */ 28 29 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); 30 31 PetscPrintf(PETSC_COMM_SELF, "[%D] Input Dim=%D IFR=[%D]; IREM=[%D]; NCOUNT=[%D]\n", rank, dim, fraction, remainder, neleadim); 32 33 starts[0]=starts[1]=starts[2]=0; /* default dimensional element = 1 */ 34 sizes[0]=sizes[1]=sizes[2]=neleglob; /* default dimensional element = 1 */ 35 36 if(rank < remainder) { 37 /* This process gets "fraction+1" elements */ 38 sizes[dim-1] = (fraction + 1); 39 starts[dim-1] = rank * (fraction+1); 40 } else { 41 /* This process gets "fraction" elements */ 42 sizes[dim-1] = fraction; 43 starts[dim-1] = (remainder*(fraction+1) + fraction*(rank-remainder)); 44 } 45 46 for(int i=dim-1; i>=0; --i) { 47 ise[2*i]=starts[i];ise[2*i+1]=starts[i]+sizes[i]; 48 PetscPrintf(PETSC_COMM_SELF, "[%D] Dim=%D ISTART=[%D]; IEND=[%D]; NCOUNT=[%D]\n", rank, i, ise[2*i], ise[2*i+1], sizes[i]); 49 } 50 PetscPrintf(PETSC_COMM_SELF, "[%D] X=[%D, %D]; Y=[%D,%D]; Z=[%D,%D]\n", rank, ise[0], ise[1], ise[2], ise[3], ise[4], ise[5]); 51 52 PetscFunctionReturn(0); 53 } 54 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "DMMoabCreateBoxMesh" 58 PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscInt nele, PetscInt nghost, DM *dm) 59 { 60 PetscErrorCode ierr; 61 moab::ErrorCode merr; 62 PetscInt rank,nprocs; 63 DM_Moab *dmmoab; 64 moab::Interface *mbiface; 65 moab::ParallelComm *pcomm; 66 moab::Tag id_tag=PETSC_NULL; 67 moab::Range range; 68 moab::EntityType etype; 69 moab::ScdInterface *scdiface; 70 PetscInt ise[6]; 71 PetscReal xse[6]; 72 73 // Determine which elements (cells) the current process owns: 74 const PetscInt npts=nele+1; 75 PetscInt my_nele,my_npts; // Number of elements owned by this process 76 PetscInt my_estart; // The starting element for this process 77 PetscInt vpere; 78 79 PetscFunctionBegin; 80 ierr = DMMoabCreateMoab(comm, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, dm);CHKERRQ(ierr); 81 82 dmmoab = (DM_Moab*)(*dm)->data; 83 mbiface = dmmoab->mbiface; 84 pcomm = dmmoab->pcomm; 85 id_tag = dmmoab->ltog_tag; 86 87 nprocs = pcomm->size(); 88 rank = pcomm->rank(); 89 90 // Begin with some error checking: 91 if(npts < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be >= 2"); 92 if(nprocs >= npts) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of processors must be less than number of points"); 93 94 // No errors,proceed with building the mesh: 95 merr = mbiface->query_interface(scdiface);MBERRNM(merr); // get a ScdInterface object through moab instance 96 97 moab::ReadUtilIface* readMeshIface; 98 mbiface->query_interface(readMeshIface); 99 100 ierr = PetscMemzero(ise,sizeof(PetscInt)*6);CHKERRQ(ierr); 101 ierr = DMMoabComputeDomainBounds_Private(pcomm, dim, nele, ise);CHKERRQ(ierr); 102 103 my_estart = ise[2*(dim-1)]; 104 switch(dim) { 105 case 1: 106 vpere = 2; 107 my_nele = (ise[1]-ise[0]); 108 my_npts = (ise[1]-ise[0]+1); 109 etype = moab::MBEDGE; 110 break; 111 case 2: 112 vpere = 4; 113 my_nele = (ise[1]-ise[0])*(ise[3]-ise[2]); 114 my_npts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1); 115 etype = moab::MBQUAD; 116 break; 117 case 3: 118 vpere = 8; 119 my_nele = (ise[1]-ise[0])*(ise[3]-ise[2])*(ise[5]-ise[4]); 120 my_npts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1)*(ise[5]-ise[4]+1); 121 etype = moab::MBHEX; 122 break; 123 } 124 125 /* we have a domain of size [1,1,1] - now compute local co-ordinate box */ 126 ierr = PetscMemzero(xse,sizeof(PetscReal)*6);CHKERRQ(ierr); 127 for(int i=0; i<6; ++i) { 128 xse[i]=(PetscReal)ise[i]/nele; 129 PetscPrintf(PETSC_COMM_SELF, "[%D] Coords %d ; nele = %D; [%D, %G]\n", rank, i, nele, ise[i], xse[i]); 130 } 131 PetscPrintf(PETSC_COMM_SELF, "[%D] Coords X=[%G, %G]; Y=[%G,%G]; Z=[%G,%G]\n", rank, xse[0], xse[1], xse[2], xse[3], xse[4], xse[5]); 132 133 PetscPrintf(PETSC_COMM_SELF, "\n[%D] My start_ele = %D and tot_nele = %D\n", rank,my_estart, my_nele); 134 135 /* 136 // Compute the co-ordinates of vertices 137 std::vector<double> coords(my_npts*vpere*3); // vertices_per_edge = 2, 3 doubles/point 138 std::vector<int> vgid(my_npts); 139 int vcount=0; 140 double hxyz=1.0/nele; 141 for (int k = ise[4]; k <= ise[5]; k++) { 142 for (int j = ise[2]; j <= ise[3]; j++) { 143 for (int i = ise[0]; i <= ise[1]; i++, vcount++) { 144 coords[vcount*3] = i*hxyz; 145 coords[vcount*3+1] = j*hxyz; 146 coords[vcount*3+2] = k*hxyz; 147 vgid[vcount] = (k*nele+j)*nele+i; 148 } 149 } 150 } 151 152 153 // 1. Creates a IxJxK structured mesh, which includes I*J*K vertices and (I-1)*(J-1)*(K-1) hexes. 154 moab::ScdBox *box; 155 moab::HomCoord low(ise[0], ise[2], ise[4]); 156 moab::HomCoord high(ise[1], ise[3], ise[5]); 157 // low.normalize(); high.normalize(); 158 merr = scdiface->construct_box(low, high, // low, high box corners in parametric space 159 coords.data(), vcount, // NULL coords vector and 0 coords (don't specify coords for now) 160 box); // box is the structured box object providing the parametric 161 // structured mesh interface for this tensor grid of elements 162 MBERRNM(merr); 163 */ 164 165 // Create vertexes and set the coodinate of each vertex: 166 moab::EntityHandle vfirst; 167 std::vector<double*> vcoords; 168 const int sequence_size = (my_nele + 2) + 1; 169 merr = readMeshIface->get_node_coords(3,my_npts,0,vfirst,vcoords,sequence_size);MBERRNM(merr); 170 171 // Compute the co-ordinates of vertices and global IDs 172 std::vector<int> vgid(my_npts); 173 int vcount=0; 174 double hxyz=1.0/nele; 175 for (int k = ise[4]; k <= ise[5]; k++) { 176 for (int j = ise[2]; j <= ise[3]; j++) { 177 for (int i = ise[0]; i <= ise[1]; i++, vcount++) { 178 vcoords[0][vcount] = i*hxyz; 179 vcoords[1][vcount] = j*hxyz; 180 vcoords[2][vcount] = k*hxyz; 181 vgid[vcount] = (k*nele+j)*nele+i; 182 } 183 } 184 } 185 186 moab::Range ownedvtx,ownedelms; 187 merr = mbiface->get_entities_by_type(0,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr); 188 189 // Get the global ID tag. The global ID tag is applied to each 190 // vertex. It acts as an global identifier which MOAB uses to 191 // assemble the individual pieces of the mesh: 192 merr = mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);MBERRNM(merr); 193 194 // set the global id for all the owned vertices 195 merr = mbiface->tag_set_data(id_tag,ownedvtx,vgid.data());MBERRNM(merr); 196 197 // Create elements between mesh points. This is done so that VisIt 198 // will interpret the output as a mesh that can be plotted... 199 moab::EntityHandle efirst; 200 moab::EntityHandle *connectivity = 0; 201 std::vector<int> subent_conn(vpere); 202 203 merr = readMeshIface->get_element_connect (my_nele,vpere,etype,1,efirst,connectivity);MBERRNM(merr); 204 205 int ecount=0; 206 for (int k = ise[4]; k < std::max(ise[5],1); k++) { 207 for (int j = ise[2]; j < std::max(ise[3],1); j++) { 208 for (int i = ise[0]; i < std::max(ise[1],1); i++,ecount++) { 209 const int offset = ecount*vpere; 210 moab::CN::SubEntityVertexIndices(etype, dim, 0, subent_conn.data()); 211 212 switch(dim) { 213 case 1: 214 connectivity[offset+subent_conn[0]] = vfirst+i; 215 connectivity[offset+subent_conn[1]] = vfirst+(i+1); 216 break; 217 case 2: 218 connectivity[offset+subent_conn[0]] = vfirst+i+j*(nele+1); 219 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+j*(nele+1); 220 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(j+1)*(nele+1); 221 connectivity[offset+subent_conn[3]] = vfirst+i+(j+1)*(nele+1); 222 break; 223 case 3: 224 connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k); 225 connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*k); 226 connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k); 227 connectivity[offset+subent_conn[3]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*k); 228 connectivity[offset+subent_conn[4]] = vfirst+i+(nele+1)*(j+(nele+1)*(k+1)); 229 connectivity[offset+subent_conn[5]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1)); 230 connectivity[offset+subent_conn[6]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*(k+1)); 231 connectivity[offset+subent_conn[7]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1)); 232 break; 233 } 234 } 235 } 236 } 237 merr = readMeshIface->update_adjacencies(efirst,my_nele,vpere,connectivity);MBERRNM(merr); 238 239 // 2. Get the vertices and hexes from moab and check their numbers against I*J*K and (I-1)*(J-1)*(K-1), resp. 240 // first '0' specifies "root set", or entire MOAB instance, second the entity dimension being requested 241 merr = mbiface->get_entities_by_dimension(0, dim, ownedelms);MBERRNM(merr); 242 243 if (my_nele != (int) ownedelms.size()) 244 SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of elements! (%D!=%D)",my_nele,ownedelms.size()); 245 else if(my_npts != (int) ownedvtx.size()) 246 SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of vertices! (%D!=%D)",my_npts,ownedvtx.size()); 247 else 248 PetscPrintf(PETSC_COMM_WORLD, "Created %D elements and %D vertices.\n", ownedelms.size(), ownedvtx.size()); 249 250 // 3. Loop over elements in 3 nested loops over i, j, k; for each (i,j,k): 251 /* 252 std::vector<double> coords(vpere*3); // vertices_per_edge = 2, 3 doubles/point 253 std::vector<moab::EntityHandle> connect; 254 int count=0; 255 for (int k = ise[4]; k < std::max(ise[5],1); k++) { 256 for (int j = ise[2]; j < std::max(ise[3],1); j++) { 257 for (int i = ise[0]; i < std::max(ise[1],1); i++,count++) { 258 // 3a. Get the element corresponding to (i,j,k) 259 moab::EntityHandle ehandle = box->get_element(i, j, k); 260 if (0 == ehandle) MBERRNM(moab::MB_FAILURE); 261 262 PetscPrintf(PETSC_COMM_SELF, "[%D] element[%D,%D,%D]=%D\n", rank, i,j,k, ehandle); 263 264 // 3b. Get the connectivity of the element 265 merr = mbiface->get_connectivity(&ehandle, 1, connect);MBERRNM(merr); // get the connectivity, in canonical order 266 267 // 3c. Get the coordinates of the vertices comprising that element 268 merr = mbiface->get_coords(connect.data(), connect.size(), coords.data());MBERRNM(merr); // get the coordinates of those vertices 269 270 for (int iv=0; iv<vpere; ++iv) 271 PetscPrintf(PETSC_COMM_SELF, "[%D] \t iv=%D [X,Y,Z]=[%G, %G, %G]\n", rank, iv, coords[iv*3], coords[iv*3+1], coords[iv*3+2]); 272 PetscPrintf(PETSC_COMM_SELF, "\n"); 273 } 274 } 275 } 276 277 merr = readMeshIface->update_adjacencies(box->get_element(ise[0],ise[2],ise[4]),my_nele,vpere,connect.data());MBERRNM(merr); 278 279 280 // 4. Release the structured mesh interface 281 mbiface->release_interface(scdiface); // tell MOAB we're done with the ScdInterface 282 */ 283 284 // The global ID tag is applied to each 285 // vertex. It acts as an global identifier which MOAB uses to 286 // assemble the individual pieces of the mesh: 287 // Set the global ID indices 288 // std::vector<int> global_ids(my_npts); 289 // for (int i = 0; i < my_npts; i++) { 290 // global_ids[i] = i+my_estart; 291 // } 292 293 // set the global id for all the owned vertices 294 // merr = mbiface->tag_set_data(id_tag,ownedvtx,global_ids.data());MBERRNM(merr); 295 296 merr = pcomm->check_all_shared_handles();MBERRNM(merr); 297 298 if (rank) 299 reinterpret_cast<moab::Core*>(mbiface)->print_database(); 300 301 // resolve the shared entities by exchanging information to adjacent processors 302 merr = mbiface->get_entities_by_type(0,etype,ownedelms,true);MBERRNM(merr); 303 merr = pcomm->resolve_shared_ents(0,ownedelms,dim,0);MBERRNM(merr); 304 305 // Reassign global IDs on all entities. 306 merr = pcomm->assign_global_ids(0,dim,0,false,true);MBERRNM(merr); 307 merr = pcomm->exchange_ghost_cells(dim,0,nghost,0,true);MBERRNM(merr); 308 309 // Everything is set up, now just do a tag exchange to update tags 310 // on all of the ghost vertexes: 311 merr = pcomm->exchange_tags(id_tag,ownedvtx);MBERRNM(merr); 312 merr = pcomm->exchange_tags(id_tag,ownedelms);MBERRNM(merr); 313 314 // set the dimension of the mesh 315 merr = mbiface->set_dimension(dim);MBERRNM(merr); 316 317 318 std::stringstream sstr; 319 sstr << "test_" << rank << ".vtk"; 320 mbiface->write_mesh(sstr.str().c_str()); 321 PetscFunctionReturn(0); 322 } 323 324 325 326 327 PetscErrorCode resolve_and_exchange(moab::ParallelComm* mbpc, PetscInt dim) 328 { 329 moab::EntityHandle entity_set; 330 moab::ErrorCode merr; 331 moab::Interface *mbint=mbpc->get_moab(); 332 moab::Range range; 333 moab::Tag tag; 334 PetscInt rank=mbpc->rank(); 335 336 // Create the entity set: 337 merr = mbint->create_meshset(moab::MESHSET_SET, entity_set);MBERRNM(merr); 338 339 // Get a list of elements in the current set: 340 merr = mbint->get_entities_by_dimension(0, dim, range, true);MBERRNM(merr); 341 342 // Add entities to the entity set: 343 merr = mbint->add_entities(entity_set, range);MBERRNM(merr); 344 345 // Add the MATERIAL_SET tag to the entity set: 346 merr = mbint->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, moab::MB_TYPE_INTEGER, tag);MBERRNM(merr); 347 merr = mbint->tag_set_data(tag, &entity_set, 1, &rank);MBERRNM(merr); 348 349 // Set up partition sets. This is where MOAB is actually told what 350 // entities each process owns: 351 merr = mbint->get_entities_by_type_and_tag(0, moab::MBENTITYSET, 352 &tag, NULL, 1, 353 mbpc->partition_sets());MBERRNM(merr); 354 355 // Finally, determine which entites are shared and exchange the 356 // ghosted entities: 357 merr = mbpc->resolve_shared_ents(0, -1, -1);MBERRNM(merr); 358 merr = mbpc->exchange_ghost_cells(-1, 0, 1, 0, true);MBERRNM(merr); 359 PetscFunctionReturn(0); 360 } 361 362 363 #undef __FUNCT__ 364 #define __FUNCT__ "DMMoabLoadFromFile_Private" 365 PetscErrorCode DMMoabLoadFromFile_Private(moab::ParallelComm* pcomm,PetscInt dim,PetscInt npts,PetscInt nghost) 366 { 367 // moab::ErrorCode merr; 368 PetscInt rank,nprocs; 369 // moab::ScdInterface *scdiface; 370 371 /* 372 // Determine which elements (cells) this process owns: 373 const PetscInt nele = npts-1; 374 PetscInt my_nele; // Number of elements owned by this process 375 PetscInt my_estart; // The starting element for this process 376 const PetscInt vertices_per_edge=2; 377 */ 378 PetscFunctionBegin; 379 MPI_Comm_size( PETSC_COMM_WORLD,&nprocs ); 380 MPI_Comm_rank( PETSC_COMM_WORLD,&rank ); 381 382 383 // std::stringstream sstr; 384 // sstr << "test_" << rank << ".vtk"; 385 // mbiface->write_mesh(sstr.str().c_str()); 386 PetscFunctionReturn(0); 387 } 388 389 390 391