xref: /petsc/src/dm/impls/moab/dmmbutil.cxx (revision c8d5365d6c6c3fe4d8ab5d2f7a087b41450baa0c)
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