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