151d15aeeSVijay Mahadevan #include <petsc-private/dmmbimpl.h> /*I "petscdm.h" I*/ 251d15aeeSVijay Mahadevan #include <petsc-private/vecimpl.h> /*I "petscdm.h" I*/ 351d15aeeSVijay Mahadevan 451d15aeeSVijay Mahadevan #include <petscdmmoab.h> 551d15aeeSVijay Mahadevan #include <MBTagConventions.hpp> 651d15aeeSVijay Mahadevan #include <moab/ReadUtilIface.hpp> 751d15aeeSVijay Mahadevan #include <moab/ScdInterface.hpp> 851d15aeeSVijay Mahadevan #include <moab/CN.hpp> 951d15aeeSVijay Mahadevan 1051d15aeeSVijay Mahadevan 1151d15aeeSVijay Mahadevan #undef __FUNCT__ 1251d15aeeSVijay Mahadevan #define __FUNCT__ "DMMoabComputeDomainBounds_Private" 1351d15aeeSVijay Mahadevan PetscErrorCode DMMoabComputeDomainBounds_Private(moab::ParallelComm* pcomm, PetscInt dim, PetscInt neleglob, PetscInt *ise) 1451d15aeeSVijay Mahadevan { 1551d15aeeSVijay Mahadevan PetscInt size,rank; 1651d15aeeSVijay Mahadevan PetscInt fraction,remainder; 1751d15aeeSVijay Mahadevan PetscInt starts[3],sizes[3]; 1851d15aeeSVijay Mahadevan 1951d15aeeSVijay Mahadevan PetscFunctionBegin; 2051d15aeeSVijay Mahadevan if(dim<1 && dim>3) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The problem dimension is invalid: %D",dim); 2151d15aeeSVijay Mahadevan 2251d15aeeSVijay Mahadevan size=pcomm->size(); 2351d15aeeSVijay Mahadevan rank=pcomm->rank(); 2451d15aeeSVijay Mahadevan fraction=neleglob/size; /* partition only by the largest dimension */ 2551d15aeeSVijay Mahadevan remainder=neleglob%size; /* remainder after partition which gets evenly distributed by round-robin */ 2651d15aeeSVijay Mahadevan 2751d15aeeSVijay Mahadevan 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); 2851d15aeeSVijay Mahadevan 2951d15aeeSVijay Mahadevan starts[0]=starts[1]=starts[2]=0; /* default dimensional element = 1 */ 3051d15aeeSVijay Mahadevan sizes[0]=sizes[1]=sizes[2]=neleglob; /* default dimensional element = 1 */ 3151d15aeeSVijay Mahadevan 3251d15aeeSVijay Mahadevan if(rank < remainder) { 3351d15aeeSVijay Mahadevan /* This process gets "fraction+1" elements */ 3451d15aeeSVijay Mahadevan sizes[dim-1] = (fraction + 1); 3551d15aeeSVijay Mahadevan starts[dim-1] = rank * (fraction+1); 3651d15aeeSVijay Mahadevan } else { 3751d15aeeSVijay Mahadevan /* This process gets "fraction" elements */ 3851d15aeeSVijay Mahadevan sizes[dim-1] = fraction; 3951d15aeeSVijay Mahadevan starts[dim-1] = (remainder*(fraction+1) + fraction*(rank-remainder)); 4051d15aeeSVijay Mahadevan } 4151d15aeeSVijay Mahadevan 4251d15aeeSVijay Mahadevan for(int i=dim-1; i>=0; --i) { 4351d15aeeSVijay Mahadevan ise[2*i]=starts[i];ise[2*i+1]=starts[i]+sizes[i]; 4451d15aeeSVijay Mahadevan } 4551d15aeeSVijay Mahadevan 4651d15aeeSVijay Mahadevan PetscFunctionReturn(0); 4751d15aeeSVijay Mahadevan } 4851d15aeeSVijay Mahadevan 49e5136372SVijay Mahadevan static void set_structured_coordinates(PetscInt i, PetscInt j, PetscInt k, PetscReal hxyz, PetscInt vcount, std::vector<double*>& vcoords) 50e5136372SVijay Mahadevan { 51e5136372SVijay Mahadevan vcoords[0][vcount] = i*hxyz; 52e5136372SVijay Mahadevan vcoords[1][vcount] = j*hxyz; 53e5136372SVijay Mahadevan vcoords[2][vcount] = k*hxyz; 54*c8d5365dSVijay Mahadevan // PetscPrintf(PETSC_COMM_SELF, " vertex - %D, %D, %D - [%G, %G]\n", i, j, k, vcoords[0][vcount], vcoords[1][vcount]); 55e5136372SVijay Mahadevan } 56e5136372SVijay Mahadevan 57e5136372SVijay Mahadevan 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) 58e5136372SVijay Mahadevan { 59e5136372SVijay Mahadevan std::vector<int> subent_conn(pow(2,dim)); /* only linear edge, quad, hex supported now */ 60e5136372SVijay Mahadevan 61e5136372SVijay Mahadevan moab::CN::SubEntityVertexIndices(etype, dim, 0, subent_conn.data()); 62e5136372SVijay Mahadevan 63e5136372SVijay Mahadevan switch(dim) { 64e5136372SVijay Mahadevan case 1: 65e5136372SVijay Mahadevan connectivity[offset+subent_conn[0]] = vfirst+i; 66e5136372SVijay Mahadevan connectivity[offset+subent_conn[1]] = vfirst+(i+1); 67e5136372SVijay Mahadevan // 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]]); 68e5136372SVijay Mahadevan break; 69e5136372SVijay Mahadevan case 2: 70e5136372SVijay Mahadevan connectivity[offset+subent_conn[0]] = vfirst+i+j*(nele+1); 71e5136372SVijay Mahadevan connectivity[offset+subent_conn[1]] = vfirst+(i+1)+j*(nele+1); 72e5136372SVijay Mahadevan connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(j+1)*(nele+1); 73e5136372SVijay Mahadevan connectivity[offset+subent_conn[3]] = vfirst+i+(j+1)*(nele+1); 74e5136372SVijay Mahadevan // 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]]); 75e5136372SVijay Mahadevan break; 76e5136372SVijay Mahadevan case 3: 77e5136372SVijay Mahadevan default: 78e5136372SVijay Mahadevan connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k); 79e5136372SVijay Mahadevan connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*k); 80e5136372SVijay Mahadevan connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k); 81e5136372SVijay Mahadevan connectivity[offset+subent_conn[3]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*k); 82e5136372SVijay Mahadevan connectivity[offset+subent_conn[4]] = vfirst+i+(nele+1)*(j+(nele+1)*(k+1)); 83e5136372SVijay Mahadevan connectivity[offset+subent_conn[5]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1)); 84e5136372SVijay Mahadevan connectivity[offset+subent_conn[6]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*(k+1)); 85e5136372SVijay Mahadevan connectivity[offset+subent_conn[7]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1)); 86e5136372SVijay Mahadevan break; 87e5136372SVijay Mahadevan } 88e5136372SVijay Mahadevan 89e5136372SVijay Mahadevan } 9051d15aeeSVijay Mahadevan 9151d15aeeSVijay Mahadevan #undef __FUNCT__ 9251d15aeeSVijay Mahadevan #define __FUNCT__ "DMMoabCreateBoxMesh" 9351d15aeeSVijay Mahadevan PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscInt nele, PetscInt nghost, DM *dm) 9451d15aeeSVijay Mahadevan { 9551d15aeeSVijay Mahadevan PetscErrorCode ierr; 9651d15aeeSVijay Mahadevan moab::ErrorCode merr; 97e5136372SVijay Mahadevan PetscInt i,j,k,n,nprocs,rank; 9851d15aeeSVijay Mahadevan DM_Moab *dmmoab; 9951d15aeeSVijay Mahadevan moab::Interface *mbiface; 10051d15aeeSVijay Mahadevan moab::ParallelComm *pcomm; 101a4717116SVijay Mahadevan moab::ReadUtilIface* readMeshIface; 102a4717116SVijay Mahadevan 10351d15aeeSVijay Mahadevan moab::Tag id_tag=PETSC_NULL; 104a4717116SVijay Mahadevan moab::Range ownedvtx,ownedelms; 105a4717116SVijay Mahadevan moab::EntityHandle vfirst,efirst; 106a4717116SVijay Mahadevan std::vector<double*> vcoords; 107a4717116SVijay Mahadevan moab::EntityHandle *connectivity = 0; 10851d15aeeSVijay Mahadevan moab::EntityType etype; 109*c8d5365dSVijay Mahadevan PetscInt ise[6]; 11051d15aeeSVijay Mahadevan PetscReal xse[6]; 11151d15aeeSVijay Mahadevan 112a4717116SVijay Mahadevan const PetscInt npts=nele+1; /* Number of points in every dimension */ 113e5136372SVijay Mahadevan PetscInt vpere,locnele,locnpts,ghnele,ghnpts; /* Number of verts/element, vertices, elements owned by this process */ 11451d15aeeSVijay Mahadevan 11551d15aeeSVijay Mahadevan PetscFunctionBegin; 116e5136372SVijay Mahadevan if(dim < 1 || dim > 3) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Invalid dimension argument for mesh: dim=[1,3].\n"); 117e5136372SVijay Mahadevan 118*c8d5365dSVijay Mahadevan ierr = MPI_Comm_size(comm, &nprocs);CHKERRQ(ierr); 119e5136372SVijay Mahadevan /* total number of vertices in all dimensions */ 120e5136372SVijay Mahadevan n=pow(npts,dim); 121e5136372SVijay Mahadevan 122e5136372SVijay Mahadevan /* do some error checking */ 123e5136372SVijay Mahadevan if(n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be >= 2.\n"); 124e5136372SVijay Mahadevan 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"); 125e5136372SVijay Mahadevan if(nghost < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of ghost layers cannot be negative.\n"); 126e5136372SVijay Mahadevan 127a4717116SVijay Mahadevan /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 12851d15aeeSVijay Mahadevan ierr = DMMoabCreateMoab(comm, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, dm);CHKERRQ(ierr); 12951d15aeeSVijay Mahadevan 130a4717116SVijay Mahadevan /* get all the necessary handles from the private DM object */ 13151d15aeeSVijay Mahadevan dmmoab = (DM_Moab*)(*dm)->data; 13251d15aeeSVijay Mahadevan mbiface = dmmoab->mbiface; 13351d15aeeSVijay Mahadevan pcomm = dmmoab->pcomm; 13451d15aeeSVijay Mahadevan id_tag = dmmoab->ltog_tag; 13551d15aeeSVijay Mahadevan nprocs = pcomm->size(); 136e5136372SVijay Mahadevan rank = pcomm->rank(); 137*c8d5365dSVijay Mahadevan dmmoab->dim = dim; 13851d15aeeSVijay Mahadevan 139a4717116SVijay Mahadevan /* No errors yet; proceed with building the mesh */ 140a4717116SVijay Mahadevan merr = mbiface->query_interface(readMeshIface);MBERRNM(merr); 14151d15aeeSVijay Mahadevan 14251d15aeeSVijay Mahadevan ierr = PetscMemzero(ise,sizeof(PetscInt)*6);CHKERRQ(ierr); 143a4717116SVijay Mahadevan 144a4717116SVijay Mahadevan /* call the collective routine that computes the domain bounds for a structured mesh using MOAB */ 14551d15aeeSVijay Mahadevan ierr = DMMoabComputeDomainBounds_Private(pcomm, dim, nele, ise);CHKERRQ(ierr); 14651d15aeeSVijay Mahadevan 147a4717116SVijay Mahadevan /* set some variables based on dimension */ 14851d15aeeSVijay Mahadevan switch(dim) { 14951d15aeeSVijay Mahadevan case 1: 15051d15aeeSVijay Mahadevan vpere = 2; 151a4717116SVijay Mahadevan locnele = (ise[1]-ise[0]); 152a4717116SVijay Mahadevan locnpts = (ise[1]-ise[0]+1); 153*c8d5365dSVijay Mahadevan ghnele = (nghost > 0 ? (ise[0] > nghost ? 1 : 0) + (ise[1] < nele - nghost ? 1 : 0) : 0 ); 154*c8d5365dSVijay Mahadevan ghnpts = (nghost > 0 ? (ise[0] > 0 ? 1 : 0) + (ise[1] < nele ? 1 : 0) : 0); 15551d15aeeSVijay Mahadevan etype = moab::MBEDGE; 15651d15aeeSVijay Mahadevan break; 15751d15aeeSVijay Mahadevan case 2: 15851d15aeeSVijay Mahadevan vpere = 4; 159a4717116SVijay Mahadevan locnele = (ise[1]-ise[0])*(ise[3]-ise[2]); 160a4717116SVijay Mahadevan locnpts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1); 161*c8d5365dSVijay Mahadevan ghnele = (nghost > 0 ? (ise[2] > 0 ? nele : 0) + (ise[3] < nele ? nele : 0) : 0); 162*c8d5365dSVijay Mahadevan ghnpts = (nghost > 0 ? (ise[2] > 0 ? npts : 0) + (ise[3] < nele ? npts : 0) : 0); 16351d15aeeSVijay Mahadevan etype = moab::MBQUAD; 16451d15aeeSVijay Mahadevan break; 16551d15aeeSVijay Mahadevan case 3: 16651d15aeeSVijay Mahadevan vpere = 8; 167a4717116SVijay Mahadevan locnele = (ise[1]-ise[0])*(ise[3]-ise[2])*(ise[5]-ise[4]); 168a4717116SVijay Mahadevan locnpts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1)*(ise[5]-ise[4]+1); 169*c8d5365dSVijay Mahadevan ghnele = (nghost > 0 ? (ise[4] > 0 ? nele*nele : 0) + (ise[5] < nele ? nele*nele : 0) : 0); 170*c8d5365dSVijay Mahadevan ghnpts = (nghost > 0 ? (ise[4] > 0 ? npts*npts : 0) + (ise[5] < nele ? npts*npts : 0) : 0); 17151d15aeeSVijay Mahadevan etype = moab::MBHEX; 17251d15aeeSVijay Mahadevan break; 17351d15aeeSVijay Mahadevan } 174*c8d5365dSVijay Mahadevan // 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]); 17551d15aeeSVijay Mahadevan 17651d15aeeSVijay Mahadevan /* we have a domain of size [1,1,1] - now compute local co-ordinate box */ 17751d15aeeSVijay Mahadevan ierr = PetscMemzero(xse,sizeof(PetscReal)*6);CHKERRQ(ierr); 178*c8d5365dSVijay Mahadevan for(i=0; i<6; ++i) { 17951d15aeeSVijay Mahadevan xse[i]=(PetscReal)ise[i]/nele; 18051d15aeeSVijay Mahadevan } 18151d15aeeSVijay Mahadevan 182a4717116SVijay Mahadevan /* Create vertexes and set the coodinate of each vertex */ 183*c8d5365dSVijay Mahadevan merr = readMeshIface->get_node_coords(3,locnpts+ghnpts,0,vfirst,vcoords,n);MBERRNM(merr); 18451d15aeeSVijay Mahadevan 185a4717116SVijay Mahadevan /* Compute the co-ordinates of vertices and global IDs */ 186*c8d5365dSVijay Mahadevan std::vector<int> vgid(locnpts+ghnpts); 18751d15aeeSVijay Mahadevan int vcount=0; 18851d15aeeSVijay Mahadevan double hxyz=1.0/nele; 189e5136372SVijay Mahadevan 190*c8d5365dSVijay Mahadevan /* create all the owned vertices */ 191*c8d5365dSVijay Mahadevan for (k = ise[4]; k <= ise[5]; k++) { 192*c8d5365dSVijay Mahadevan for (j = ise[2]; j <= ise[3]; j++) { 193*c8d5365dSVijay Mahadevan for (i = ise[0]; i <= ise[1]; i++, vcount++) { 194*c8d5365dSVijay Mahadevan // PetscPrintf(PETSC_COMM_SELF, "[%D] creating owned", rank); 195*c8d5365dSVijay Mahadevan set_structured_coordinates(i,j,k,hxyz,vcount,vcoords); 196*c8d5365dSVijay Mahadevan vgid[vcount] = (k*npts+j)*npts+i+1; 197*c8d5365dSVijay Mahadevan } 198*c8d5365dSVijay Mahadevan } 199*c8d5365dSVijay Mahadevan } 200*c8d5365dSVijay Mahadevan 201e5136372SVijay Mahadevan /* create ghosted vertices requested by user - below the current plane */ 202e5136372SVijay Mahadevan if (ise[2*dim-2] > 0) { 203e5136372SVijay Mahadevan for (k = (dim==3?ise[4]-nghost:ise[4]); k <= (dim==3?ise[4]-1:ise[5]); k++) { 204e5136372SVijay Mahadevan for (j = (dim==2?ise[2]-nghost:ise[2]); j <= (dim==2?ise[2]-1:ise[3]); j++) { 205e5136372SVijay Mahadevan for (i = (dim>1?ise[0]:ise[0]-nghost); i <= (dim>1?ise[1]:ise[0]-1); i++, vcount++) { 206*c8d5365dSVijay Mahadevan // PetscPrintf(PETSC_COMM_SELF, "[%D] creating ghost", rank); 207e5136372SVijay Mahadevan set_structured_coordinates(i,j,k,hxyz,vcount,vcoords); 208e5136372SVijay Mahadevan vgid[vcount] = (k*npts+j)*npts+i+1; 209e5136372SVijay Mahadevan } 210e5136372SVijay Mahadevan } 211e5136372SVijay Mahadevan } 212e5136372SVijay Mahadevan } 213e5136372SVijay Mahadevan 214e5136372SVijay Mahadevan /* create ghosted vertices requested by user - above the current plane */ 215e5136372SVijay Mahadevan if (ise[2*dim-1] < nele) { 216e5136372SVijay Mahadevan for (k = (dim==3?ise[5]+1:ise[4]); k <= (dim==3?ise[5]+nghost:ise[5]); k++) { 217e5136372SVijay Mahadevan for (j = (dim==2?ise[3]+1:ise[2]); j <= (dim==2?ise[3]+nghost:ise[3]); j++) { 218e5136372SVijay Mahadevan for (i = (dim>1?ise[0]:ise[1]+1); i <= (dim>1?ise[1]:ise[1]+nghost); i++, vcount++) { 219*c8d5365dSVijay Mahadevan // PetscPrintf(PETSC_COMM_SELF, "[%D] creating ghost", rank); 220e5136372SVijay Mahadevan set_structured_coordinates(i,j,k,hxyz,vcount,vcoords); 221e5136372SVijay Mahadevan vgid[vcount] = (k*npts+j)*npts+i+1; 222e5136372SVijay Mahadevan } 22351d15aeeSVijay Mahadevan } 22451d15aeeSVijay Mahadevan } 22551d15aeeSVijay Mahadevan } 22651d15aeeSVijay Mahadevan 227*c8d5365dSVijay Mahadevan if (locnpts+ghnpts != vcount) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of vertices! (%D!=%D)",locnpts+ghnpts,vcount); 228*c8d5365dSVijay Mahadevan 22951d15aeeSVijay Mahadevan merr = mbiface->get_entities_by_type(0,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr); 230a4717116SVijay Mahadevan merr = mbiface->add_entities (dmmoab->fileset, ownedvtx);MBERRNM(merr); 231*c8d5365dSVijay Mahadevan // merr = mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr); 23251d15aeeSVijay Mahadevan 233a4717116SVijay Mahadevan /* The global ID tag is applied to each owned 234a4717116SVijay Mahadevan vertex. It acts as an global identifier which MOAB uses to 235a4717116SVijay Mahadevan assemble the individual pieces of the mesh: 236a4717116SVijay Mahadevan Set the global ID indices */ 23751d15aeeSVijay Mahadevan merr = mbiface->tag_set_data(id_tag,ownedvtx,vgid.data());MBERRNM(merr); 23851d15aeeSVijay Mahadevan 239a4717116SVijay Mahadevan /* Create elements between mesh points using the ReadUtilInterface 240a4717116SVijay Mahadevan get the reference to element connectivities for all local elements from the ReadUtilInterface */ 241e5136372SVijay Mahadevan merr = readMeshIface->get_element_connect (locnele+ghnele,vpere,etype,1,efirst,connectivity);MBERRNM(merr); 24251d15aeeSVijay Mahadevan 243a4717116SVijay Mahadevan /* offset appropriately so that only local ID and not global ID numbers are set for connectivity array */ 244*c8d5365dSVijay Mahadevan // PetscPrintf(PETSC_COMM_SELF, "[%D] first local handle %D\n", rank, vfirst); 245e5136372SVijay Mahadevan vfirst-=vgid[0]-1; 24651d15aeeSVijay Mahadevan 247a4717116SVijay Mahadevan /* 3. Loop over elements in 3 nested loops over i, j, k; for each (i,j,k): 248a4717116SVijay Mahadevan and then set the connectivity for each element appropriately */ 24951d15aeeSVijay Mahadevan int ecount=0; 25051d15aeeSVijay Mahadevan 251e5136372SVijay Mahadevan /* create ghosted elements requested by user - below the current plane */ 252e5136372SVijay Mahadevan if (ise[2*dim-2] >= nghost) { 253e5136372SVijay Mahadevan for (k = (dim==3?ise[4]-nghost:ise[4]); k < (dim==3?ise[4]:std::max(ise[5],1)); k++) { 254e5136372SVijay Mahadevan for (j = (dim==2?ise[2]-nghost:ise[2]); j < (dim==2?ise[2]:std::max(ise[3],1)); j++) { 255e5136372SVijay Mahadevan for (i = (dim>1?ise[0]:ise[0]-nghost); i < (dim>1?std::max(ise[1],1):ise[0]); i++, ecount++) { 256*c8d5365dSVijay Mahadevan // PetscPrintf(PETSC_COMM_SELF, "[%D] creating ghost element - %D, %D, %D\n", rank, i, j, k); 257*c8d5365dSVijay Mahadevan set_element_connectivity(dim, etype, ecount*vpere, nele, i, j, k, vfirst, connectivity); 25851d15aeeSVijay Mahadevan } 25951d15aeeSVijay Mahadevan } 26051d15aeeSVijay Mahadevan } 26151d15aeeSVijay Mahadevan } 262e5136372SVijay Mahadevan 263e5136372SVijay Mahadevan /* create owned elements requested by user */ 264e5136372SVijay Mahadevan for (k = ise[4]; k < std::max(ise[5],1); k++) { 265e5136372SVijay Mahadevan for (j = ise[2]; j < std::max(ise[3],1); j++) { 266e5136372SVijay Mahadevan for (i = ise[0]; i < std::max(ise[1],1); i++,ecount++) { 267*c8d5365dSVijay Mahadevan // PetscPrintf(PETSC_COMM_SELF, "[%D] creating owned element - %D, %D, %D\n", rank, i, j, k); 268*c8d5365dSVijay Mahadevan set_element_connectivity(dim, etype, ecount*vpere, nele, i, j, k, vfirst, connectivity); 269e5136372SVijay Mahadevan } 270e5136372SVijay Mahadevan } 271e5136372SVijay Mahadevan } 272e5136372SVijay Mahadevan 273e5136372SVijay Mahadevan /* create ghosted elements requested by user - above the current plane */ 274e5136372SVijay Mahadevan if (ise[2*dim-1] <= nele-nghost) { 275e5136372SVijay Mahadevan for (k = (dim==3?ise[5]:ise[4]); k < (dim==3?ise[5]+nghost:std::max(ise[5],1)); k++) { 276e5136372SVijay Mahadevan for (j = (dim==2?ise[3]:ise[2]); j < (dim==2?ise[3]+nghost:std::max(ise[3],1)); j++) { 277e5136372SVijay Mahadevan for (i = (dim>1?ise[0]:ise[1]); i < (dim>1?std::max(ise[1],1):ise[1]+nghost); i++, ecount++) { 278*c8d5365dSVijay Mahadevan // PetscPrintf(PETSC_COMM_SELF, "[%D] creating ghost element - %D, %D, %D\n", rank, i, j, k); 279*c8d5365dSVijay Mahadevan set_element_connectivity(dim, etype, ecount*vpere, nele, i, j, k, vfirst, connectivity); 280e5136372SVijay Mahadevan } 281e5136372SVijay Mahadevan } 282e5136372SVijay Mahadevan } 283e5136372SVijay Mahadevan } 284e5136372SVijay Mahadevan 285e5136372SVijay Mahadevan merr = readMeshIface->update_adjacencies(efirst,locnele+ghnele,vpere,connectivity);MBERRNM(merr); 28651d15aeeSVijay Mahadevan 287a4717116SVijay Mahadevan /* 2. Get the vertices and hexes from moab and check their numbers against I*J*K and (I-1)*(J-1)*(K-1), resp. 288a4717116SVijay Mahadevan first '0' specifies "root set", or entire MOAB instance, second the entity dimension being requested */ 28951d15aeeSVijay Mahadevan merr = mbiface->get_entities_by_dimension(0, dim, ownedelms);MBERRNM(merr); 290a4717116SVijay Mahadevan merr = mbiface->add_entities (dmmoab->fileset, ownedelms);MBERRNM(merr); 291*c8d5365dSVijay Mahadevan // merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, ownedelms);MBERRNM(merr); 29251d15aeeSVijay Mahadevan 293e5136372SVijay Mahadevan // merr = mbiface->unite_meshset(0, dmmoab->fileset);MBERRV(mbiface,merr); 294a4717116SVijay Mahadevan 295e5136372SVijay Mahadevan if (locnele+ghnele != (int) ownedelms.size()) 296e5136372SVijay Mahadevan SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of elements! (%D!=%D)",locnele+ghnele,ownedelms.size()); 297e5136372SVijay Mahadevan else if(locnpts+ghnpts != (int) ownedvtx.size()) 298e5136372SVijay Mahadevan SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of vertices! (%D!=%D)",locnpts+ghnpts,ownedvtx.size()); 29951d15aeeSVijay Mahadevan else 300a4717116SVijay Mahadevan PetscInfo2(*dm, "Created %D elements and %D vertices.\n", ownedelms.size(), ownedvtx.size()); 30151d15aeeSVijay Mahadevan 302a4717116SVijay Mahadevan /* check the handles */ 303a4717116SVijay Mahadevan merr = pcomm->check_all_shared_handles();MBERRV(mbiface,merr); 30451d15aeeSVijay Mahadevan 305*c8d5365dSVijay Mahadevan #if 0 306e5136372SVijay Mahadevan std::stringstream sstr; 307e5136372SVijay Mahadevan sstr << "test_" << pcomm->rank() << ".vtk"; 308e5136372SVijay Mahadevan mbiface->write_mesh(sstr.str().c_str()); 309e5136372SVijay Mahadevan #endif 310e5136372SVijay Mahadevan 311a4717116SVijay Mahadevan /* resolve the shared entities by exchanging information to adjacent processors */ 312a4717116SVijay Mahadevan merr = mbiface->get_entities_by_type(dmmoab->fileset,etype,ownedelms,true);MBERRNM(merr); 313a4717116SVijay Mahadevan merr = pcomm->resolve_shared_ents(dmmoab->fileset,ownedelms,dim,dim-1,&id_tag);MBERRV(mbiface,merr); 31451d15aeeSVijay Mahadevan 315a4717116SVijay Mahadevan /* Reassign global IDs on all entities. */ 316e5136372SVijay Mahadevan merr = pcomm->assign_global_ids(dmmoab->fileset,dim,1,true,true,true);MBERRNM(merr); 31751d15aeeSVijay Mahadevan 318*c8d5365dSVijay Mahadevan merr = pcomm->exchange_ghost_cells(dim,0,1/*nghost*/,dim,true,false,&dmmoab->fileset);MBERRV(mbiface,merr); 319*c8d5365dSVijay Mahadevan 320a4717116SVijay Mahadevan /* Everything is set up, now just do a tag exchange to update tags 321a4717116SVijay Mahadevan on all of the ghost vertexes */ 322*c8d5365dSVijay Mahadevan merr = mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr); 323*c8d5365dSVijay Mahadevan merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, ownedelms);MBERRNM(merr); 324*c8d5365dSVijay Mahadevan 325a4717116SVijay Mahadevan merr = pcomm->exchange_tags(id_tag,ownedvtx);MBERRV(mbiface,merr); 326a4717116SVijay Mahadevan merr = pcomm->exchange_tags(id_tag,ownedelms);MBERRV(mbiface,merr); 32751d15aeeSVijay Mahadevan 328*c8d5365dSVijay Mahadevan // if(rank) { 329*c8d5365dSVijay Mahadevan // ownedvtx.print(0); 330*c8d5365dSVijay Mahadevan // ownedelms.print(0); 331*c8d5365dSVijay Mahadevan // } 33251d15aeeSVijay Mahadevan 33351d15aeeSVijay Mahadevan PetscFunctionReturn(0); 33451d15aeeSVijay Mahadevan } 33551d15aeeSVijay Mahadevan 33651d15aeeSVijay Mahadevan 337a4717116SVijay Mahadevan #undef __FUNCT__ 338a4717116SVijay Mahadevan #define __FUNCT__ "DMMoab_GetReadOptions_Private" 3397023aa44SVijay Mahadevan PetscErrorCode DMMoab_GetReadOptions_Private(PetscBool by_rank, PetscInt numproc, PetscInt dim, MoabReadMode mode, PetscInt dbglevel, const char* extra_opts, const char** read_opts) 34051d15aeeSVijay Mahadevan { 341a4717116SVijay Mahadevan std::ostringstream str; 34251d15aeeSVijay Mahadevan 343a4717116SVijay Mahadevan PetscFunctionBegin; 344a4717116SVijay Mahadevan // do parallel read unless only one processor 345a4717116SVijay Mahadevan if (numproc > 1) { 3467023aa44SVijay Mahadevan str << "PARALLEL=" << mode << ";PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;"; 3477023aa44SVijay Mahadevan str << "PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=" << dim << ".0.1;"; 348a4717116SVijay Mahadevan if (by_rank) 349a4717116SVijay Mahadevan str << "PARTITION_BY_RANK;"; 350a4717116SVijay Mahadevan } 35151d15aeeSVijay Mahadevan 352*c8d5365dSVijay Mahadevan if (dbglevel) { 353*c8d5365dSVijay Mahadevan str << "CPUTIME;DEBUG_IO=" << dbglevel << ";"; 354*c8d5365dSVijay Mahadevan if (numproc>1) 355*c8d5365dSVijay Mahadevan str << "DEBUG_PIO=" << dbglevel << ";"; 356*c8d5365dSVijay Mahadevan } 35751d15aeeSVijay Mahadevan 358*c8d5365dSVijay Mahadevan if (extra_opts) 359*c8d5365dSVijay Mahadevan str << extra_opts; 36051d15aeeSVijay Mahadevan 3617023aa44SVijay Mahadevan *read_opts = str.str().c_str(); 36251d15aeeSVijay Mahadevan PetscFunctionReturn(0); 36351d15aeeSVijay Mahadevan } 36451d15aeeSVijay Mahadevan 36551d15aeeSVijay Mahadevan 36651d15aeeSVijay Mahadevan #undef __FUNCT__ 367a4717116SVijay Mahadevan #define __FUNCT__ "DMMoabLoadFromFile" 368a4717116SVijay Mahadevan PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm,PetscInt dim,const char* filename, const char* usrreadopts, DM *dm) 36951d15aeeSVijay Mahadevan { 370a4717116SVijay Mahadevan moab::ErrorCode merr; 3717023aa44SVijay Mahadevan PetscInt nprocs,dbglevel=0; 3727023aa44SVijay Mahadevan PetscBool part_by_rank=PETSC_FALSE; 373a4717116SVijay Mahadevan DM_Moab *dmmoab; 374a4717116SVijay Mahadevan moab::Interface *mbiface; 375a4717116SVijay Mahadevan moab::ParallelComm *pcomm; 376a4717116SVijay Mahadevan moab::Range verts,elems; 3777023aa44SVijay Mahadevan const char *readopts; 378a4717116SVijay Mahadevan PetscErrorCode ierr; 37951d15aeeSVijay Mahadevan 38051d15aeeSVijay Mahadevan PetscFunctionBegin; 381a4717116SVijay Mahadevan PetscValidPointer(dm,5); 38251d15aeeSVijay Mahadevan 383a4717116SVijay Mahadevan /* Create the basic DMMoab object and keep the default parameters created by DM impls */ 384a4717116SVijay Mahadevan ierr = DMMoabCreateMoab(comm, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, dm);CHKERRQ(ierr); 38551d15aeeSVijay Mahadevan 386a4717116SVijay Mahadevan /* get all the necessary handles from the private DM object */ 387a4717116SVijay Mahadevan dmmoab = (DM_Moab*)(*dm)->data; 388a4717116SVijay Mahadevan mbiface = dmmoab->mbiface; 389a4717116SVijay Mahadevan pcomm = dmmoab->pcomm; 390a4717116SVijay Mahadevan nprocs = pcomm->size(); 391*c8d5365dSVijay Mahadevan dmmoab->dim = dim; 392a4717116SVijay Mahadevan 3937023aa44SVijay Mahadevan /* TODO: Use command-line options to control by_rank, verbosity, MoabReadMode and extra options */ 3947023aa44SVijay Mahadevan ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for reading/writing MOAB based meshes from file", "DMMoab"); 3957023aa44SVijay Mahadevan ierr = PetscOptionsInt("-dmmb_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "dmmbutil.cxx", dbglevel, &dbglevel, NULL);CHKERRQ(ierr); 3967023aa44SVijay Mahadevan 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); 3977023aa44SVijay Mahadevan ierr = PetscOptionsEnd(); 3987023aa44SVijay Mahadevan 399a4717116SVijay Mahadevan /* add mesh loading options specific to the DM */ 4007023aa44SVijay Mahadevan ierr = DMMoab_GetReadOptions_Private(part_by_rank, nprocs, dim, MOAB_PARROPTS_READ_PART, dbglevel, usrreadopts, &readopts);CHKERRQ(ierr); 4017023aa44SVijay Mahadevan 402e5136372SVijay Mahadevan PetscInfo2(*dm, "Reading file %s with options: %s\n",filename,readopts); 403a4717116SVijay Mahadevan 404a4717116SVijay Mahadevan /* Load the mesh from a file. */ 4057023aa44SVijay Mahadevan merr = mbiface->load_file(filename, &dmmoab->fileset, readopts);MBERRVM(mbiface,"Reading MOAB file failed.", merr); 406a4717116SVijay Mahadevan 4076e40195eSVijay Mahadevan /* Reassign global IDs on all entities. */ 408e5136372SVijay Mahadevan merr = pcomm->assign_global_ids(dmmoab->fileset,dim,1,true,true,true);MBERRNM(merr); 409e5136372SVijay Mahadevan 410e5136372SVijay Mahadevan /* load the local vertices */ 411e5136372SVijay Mahadevan merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);MBERRNM(merr); 412e5136372SVijay Mahadevan /* load the local elements */ 413e5136372SVijay Mahadevan merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);MBERRNM(merr); 414e5136372SVijay Mahadevan 415e5136372SVijay Mahadevan /* Everything is set up, now just do a tag exchange to update tags 416e5136372SVijay Mahadevan on all of the ghost vertexes */ 417e5136372SVijay Mahadevan merr = pcomm->exchange_tags(dmmoab->ltog_tag,verts);MBERRV(mbiface,merr); 418e5136372SVijay Mahadevan merr = pcomm->exchange_tags(dmmoab->ltog_tag,elems);MBERRV(mbiface,merr); 419e5136372SVijay Mahadevan 420e5136372SVijay Mahadevan merr = pcomm->exchange_ghost_cells(dim,0,1,0,true,true,&dmmoab->fileset);MBERRV(mbiface,merr); 4216e40195eSVijay Mahadevan 422a4717116SVijay Mahadevan merr = pcomm->collective_sync_partition();MBERR("Collective sync failed", merr); 423a4717116SVijay Mahadevan 424e5136372SVijay Mahadevan PetscInfo3(*dm, "MOAB file '%s' was successfully loaded. Found %D vertices and %D elements.\n", filename, verts.size(), elems.size()); 425e5136372SVijay Mahadevan 426e5136372SVijay Mahadevan #if 1 427e5136372SVijay Mahadevan std::stringstream sstr; 428e5136372SVijay Mahadevan sstr << "test_" << pcomm->rank() << ".vtk"; 429e5136372SVijay Mahadevan mbiface->write_mesh(sstr.str().c_str()); 430e5136372SVijay Mahadevan #endif 431e5136372SVijay Mahadevan 43251d15aeeSVijay Mahadevan PetscFunctionReturn(0); 43351d15aeeSVijay Mahadevan } 43451d15aeeSVijay Mahadevan 435