Actual source code: dmmbutil.cxx
petsc-3.7.3 2016-08-01
1: #include <petsc/private/dmmbimpl.h> /*I "petscdmmoab.h" I*/
2: #include <petsc/private/vecimpl.h>
4: #include <petscdmmoab.h>
5: #include <MBTagConventions.hpp>
6: #include <moab/ReadUtilIface.hpp>
7: #include <moab/ScdInterface.hpp>
8: #include <moab/CN.hpp>
13: static 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];
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 */
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);
29: starts[0]=starts[1]=starts[2]=0; /* default dimensional element = 1 */
30: sizes[0]=sizes[1]=sizes[2]=neleglob; /* default dimensional element = 1 */
32: if(rank < remainder) { /* This process gets "fraction+1" elements */
33: sizes[dim-1] = (fraction + 1);
34: starts[dim-1] = rank * (fraction+1);
35: } else { /* This process gets "fraction" elements */
36: sizes[dim-1] = fraction;
37: starts[dim-1] = (remainder*(fraction+1) + fraction*(rank-remainder));
38: }
40: for(int i=dim-1; i>=0; --i) {
41: ise[2*i]=starts[i];ise[2*i+1]=starts[i]+sizes[i];
42: }
43: return(0);
44: }
48: static void DMMoab_SetStructuredCoords_Private(PetscInt i, PetscInt j, PetscInt k, PetscReal hx, PetscReal hy, PetscReal hz, PetscInt vcount, std::vector<double*>& vcoords)
49: {
50: vcoords[0][vcount] = i*hx;
51: vcoords[1][vcount] = j*hy;
52: vcoords[2][vcount] = k*hz;
53: }
57: static void DMMoab_SetTensorElementConnectivity_Private(PetscInt dim, moab::EntityType etype, PetscInt offset, PetscInt nele, PetscInt i, PetscInt j, PetscInt k, PetscInt vfirst, moab::EntityHandle *connectivity)
58: {
59: std::vector<int> subent_conn(pow(2,dim)); /* only linear edge, quad, hex supported now */
61: moab::CN::SubEntityVertexIndices(etype, dim, 0, subent_conn.data());
63: switch(dim) {
64: case 1:
65: connectivity[offset+subent_conn[0]] = vfirst+i;
66: connectivity[offset+subent_conn[1]] = vfirst+(i+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: break;
74: case 3:
75: default:
76: connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k);
77: connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*k);
78: connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k);
79: connectivity[offset+subent_conn[3]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*k);
80: connectivity[offset+subent_conn[4]] = vfirst+i+(nele+1)*(j+(nele+1)*(k+1));
81: connectivity[offset+subent_conn[5]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1));
82: connectivity[offset+subent_conn[6]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*(k+1));
83: connectivity[offset+subent_conn[7]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1));
84: break;
85: }
86: }
90: static void DMMoab_SetSimplexElementConnectivity_Private(PetscInt dim, PetscInt subelem, moab::EntityType etype, PetscInt offset, PetscInt nele, PetscInt i, PetscInt j, PetscInt k, PetscInt vfirst, moab::EntityHandle *connectivity)
91: {
92: std::vector<int> subent_conn(pow(2,dim)); /* only linear edge, quad, hex supported now */
94: moab::CN::SubEntityVertexIndices(etype, dim, 0, subent_conn.data());
96: switch(dim) {
97: case 1:
98: connectivity[offset+subent_conn[0]] = vfirst+i;
99: connectivity[offset+subent_conn[1]] = vfirst+(i+1);
100: break;
101: case 2:
102: if (subelem) { /* 1 2 3 of a QUAD */
103: connectivity[offset+subent_conn[0]] = vfirst+i+j*(nele+1);
104: connectivity[offset+subent_conn[1]] = vfirst+(i+1)+j*(nele+1);
105: connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(j+1)*(nele+1);
106: }
107: else { /* 3 4 1 of a QUAD */
108: connectivity[offset+subent_conn[0]] = vfirst+(i+1)+(j+1)*(nele+1);
109: connectivity[offset+subent_conn[1]] = vfirst+i+(j+1)*(nele+1);
110: connectivity[offset+subent_conn[2]] = vfirst+i+j*(nele+1);
111: }
112: break;
113: case 3:
114: default:
115: switch(subelem) {
116: case 0: /* 0 1 2 5 of a HEX */
117: connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k);
118: connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*k);
119: connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k);
120: connectivity[offset+subent_conn[3]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1));
121: break;
122: case 1: /* 0 2 7 5 of a HEX */
123: connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k);
124: connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k);
125: connectivity[offset+subent_conn[2]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1));
126: connectivity[offset+subent_conn[3]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1));
127: break;
128: case 2: /* 0 2 3 7 of a HEX */
129: connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k);
130: connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k);
131: connectivity[offset+subent_conn[2]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*k);
132: connectivity[offset+subent_conn[3]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1));
133: break;
134: case 3: /* 0 5 7 4 of a HEX */
135: connectivity[offset+subent_conn[0]] = vfirst+i+(nele+1)*(j+(nele+1)*k);
136: connectivity[offset+subent_conn[1]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1));
137: connectivity[offset+subent_conn[2]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1));
138: connectivity[offset+subent_conn[3]] = vfirst+i+(nele+1)*(j+(nele+1)*(k+1));
139: break;
140: case 4: /* 2 7 5 6 of a HEX */
141: connectivity[offset+subent_conn[0]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*k);
142: connectivity[offset+subent_conn[1]] = vfirst+i+(nele+1)*((j+1)+(nele+1)*(k+1));
143: connectivity[offset+subent_conn[2]] = vfirst+(i+1)+(nele+1)*(j+(nele+1)*(k+1));
144: connectivity[offset+subent_conn[3]] = vfirst+(i+1)+(nele+1)*((j+1)+(nele+1)*(k+1));
145: break;
146: }
147: break;
148: }
149: }
153: static void DMMoab_SetElementConnectivity_Private(PetscBool useSimplex, PetscInt dim, moab::EntityType etype, PetscInt *ecount, PetscInt vpere, PetscInt nele, PetscInt i, PetscInt j, PetscInt k, PetscInt vfirst, moab::EntityHandle *connectivity)
154: {
155: PetscInt m,subelem;
156: if (useSimplex) {
157: subelem=(dim==1 ? 1 : (dim==2 ? 2 : 5));
158: for (m=0; m<subelem; m++)
159: DMMoab_SetSimplexElementConnectivity_Private(dim, m, etype, (*ecount+m)*vpere, nele, i, j, k, vfirst, connectivity);
161: }
162: else {
163: subelem=1;
164: DMMoab_SetTensorElementConnectivity_Private(dim, etype, (*ecount)*vpere, nele, i, j, k, vfirst, connectivity);
165: }
166: *ecount+=subelem;
167: }
172: /*@
173: DMMoabCreateBoxMesh - Creates a mesh on the tensor product (box) of intervals with user specified bounds.
175: Collective on MPI_Comm
177: Input Parameters:
178: + comm - The communicator for the DM object
179: . dim - The spatial dimension
180: . bounds - The bounds of the box specified with [x-left, x-right, y-bottom, y-top, z-bottom, z-top] depending on the spatial dimension
181: . nele - The number of discrete elements in each direction
182: . user_nghost - The number of ghosted layers needed in the partitioned mesh
184: Output Parameter:
185: . dm - The DM object
187: Level: beginner
189: .keywords: DM, create
190: .seealso: DMSetType(), DMCreate(), DMMoabLoadFromFile()
191: @*/
192: PetscErrorCode DMMoabCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool useSimplex, const PetscReal* bounds, PetscInt nele, PetscInt user_nghost, DM *dm)
193: {
194: PetscErrorCode ierr;
195: moab::ErrorCode merr;
196: PetscInt i,j,k,n,nprocs;
197: DM_Moab *dmmoab;
198: moab::Interface *mbiface;
199: moab::ParallelComm *pcomm;
200: moab::ReadUtilIface* readMeshIface;
202: moab::Tag id_tag=PETSC_NULL;
203: moab::Range ownedvtx,ownedelms;
204: moab::EntityHandle vfirst,efirst,regionset,faceset,edgeset,vtxset;
205: std::vector<double*> vcoords;
206: moab::EntityHandle *connectivity = 0;
207: moab::EntityType etype=moab::MBHEX;
208: PetscInt ise[6];
209: PetscReal xse[6],defbounds[6];
210: /* TODO: Fix nghost > 0 - now relying on exchange_ghost_cells */
211: const PetscInt nghost=0;
213: moab::Tag geom_tag;
215: moab::Range adj,dim3,dim2;
216: bool build_adjacencies=false;
218: const PetscInt npts=nele+1; /* Number of points in every dimension */
219: PetscInt vpere=0,locnele=0,locnpts=0,ghnele,ghnpts; /* Number of verts/element, vertices, elements owned by this process */
222: if(dim < 1 || dim > 3) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Invalid dimension argument for mesh: dim=[1,3].\n");
224: MPI_Comm_size(comm, &nprocs);
225: /* total number of vertices in all dimensions */
226: n=pow(npts,dim);
228: /* do some error checking */
229: if(n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of points must be >= 2.\n");
230: 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");
231: if(nghost < 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Number of ghost layers cannot be negative.\n");
233: /* Create the basic DMMoab object and keep the default parameters created by DM impls */
234: DMMoabCreateMoab(comm, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, dm);
236: /* get all the necessary handles from the private DM object */
237: dmmoab = (DM_Moab*)(*dm)->data;
238: mbiface = dmmoab->mbiface;
239: pcomm = dmmoab->pcomm;
240: id_tag = dmmoab->ltog_tag;
241: nprocs = pcomm->size();
242: dmmoab->dim = dim;
244: /* create a file set to associate all entities in current mesh */
245: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, dmmoab->fileset);MBERR("Creating file set failed", merr);
247: /* No errors yet; proceed with building the mesh */
248: merr = mbiface->query_interface(readMeshIface);MBERRNM(merr);
250: PetscMemzero(ise,sizeof(PetscInt)*6);
252: /* call the collective routine that computes the domain bounds for a structured mesh using MOAB */
253: DMMoabComputeDomainBounds_Private(pcomm, dim, nele, ise);
255: /* set some variables based on dimension */
256: switch(dim) {
257: case 1:
258: vpere = 2;
259: locnele = (ise[1]-ise[0]);
260: locnpts = (ise[1]-ise[0]+1);
261: ghnele = (nghost > 0 ? (ise[0] > nghost ? 1 : 0) + (ise[1] < nele - nghost ? 1 : 0) : 0 );
262: ghnpts = (nghost > 0 ? (ise[0] > 0 ? 1 : 0) + (ise[1] < nele ? 1 : 0) : 0);
263: etype = moab::MBEDGE;
264: break;
265: case 2:
266: locnpts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1);
267: ghnpts = (nghost > 0 ? (ise[2] > 0 ? npts : 0) + (ise[3] < nele ? npts : 0) : 0);
268: if (useSimplex) {
269: vpere = 3;
270: locnele = 2*(ise[1]-ise[0])*(ise[3]-ise[2]);
271: ghnele = 2*(nghost > 0 ? (ise[2] > 0 ? nele : 0) + (ise[3] < nele ? nele : 0) : 0);
272: etype = moab::MBTRI;
273: }
274: else {
275: vpere = 4;
276: locnele = (ise[1]-ise[0])*(ise[3]-ise[2]);
277: ghnele = (nghost > 0 ? (ise[2] > 0 ? nele : 0) + (ise[3] < nele ? nele : 0) : 0);
278: etype = moab::MBQUAD;
279: }
280: break;
281: case 3:
282: default:
283: locnpts = (ise[1]-ise[0]+1)*(ise[3]-ise[2]+1)*(ise[5]-ise[4]+1);
284: ghnpts = (nghost > 0 ? (ise[4] > 0 ? npts*npts : 0) + (ise[5] < nele ? npts*npts : 0) : 0);
285: if (useSimplex) {
286: vpere = 4;
287: locnele = 5*(ise[1]-ise[0])*(ise[3]-ise[2])*(ise[5]-ise[4]);
288: ghnele = 5*(nghost > 0 ? (ise[4] > 0 ? nele*nele : 0) + (ise[5] < nele ? nele*nele : 0) : 0);
289: etype = moab::MBTET;
290: }
291: else {
292: vpere = 8;
293: locnele = (ise[1]-ise[0])*(ise[3]-ise[2])*(ise[5]-ise[4]);
294: ghnele = (nghost > 0 ? (ise[4] > 0 ? nele*nele : 0) + (ise[5] < nele ? nele*nele : 0) : 0);
295: etype = moab::MBHEX;
296: }
297: break;
298: }
300: /* we have a domain of size [1,1,1] - now compute local co-ordinate box */
301: PetscMemzero(xse,sizeof(PetscReal)*6);
302: for(i=0; i<6; ++i) {
303: xse[i]=(PetscReal)ise[i]/nele;
304: }
306: /* Create vertexes and set the coodinate of each vertex */
307: merr = readMeshIface->get_node_coords(3,locnpts+ghnpts,0,vfirst,vcoords,n);MBERRNM(merr);
309: /* Compute the co-ordinates of vertices and global IDs */
310: std::vector<int> vgid(locnpts+ghnpts);
311: int vcount=0;
313: if (!bounds) { /* default box mesh is defined on a unit-cube */
314: defbounds[0]=0.0; defbounds[1]=1.0;
315: defbounds[2]=0.0; defbounds[3]=1.0;
316: defbounds[4]=0.0; defbounds[5]=1.0;
317: bounds=defbounds;
318: }
319: else {
320: /* validate the bounds data */
321: if(bounds[0] >= bounds[1]) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"X-dim: Left boundary cannot be greater than right. [%G >= %G]\n",bounds[0],bounds[1]);
322: if(dim > 1 && (bounds[2] >= bounds[3])) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Y-dim: Left boundary cannot be greater than right. [%G >= %G]\n",bounds[2],bounds[3]);
323: if(dim > 2 && (bounds[4] >= bounds[5])) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Z-dim: Left boundary cannot be greater than right. [%G >= %G]\n",bounds[4],bounds[5]);
324: }
326: const double hx=(bounds[1]-bounds[0])/nele;
327: const double hy=(dim > 1 ? (bounds[3]-bounds[2])/nele : 0.0);
328: const double hz=(dim > 2 ? (bounds[5]-bounds[4])/nele : 0.0);
330: /* create all the owned vertices */
331: for (k = ise[4]; k <= ise[5]; k++) {
332: for (j = ise[2]; j <= ise[3]; j++) {
333: for (i = ise[0]; i <= ise[1]; i++, vcount++) {
334: DMMoab_SetStructuredCoords_Private(i,j,k,hx,hy,hz,vcount,vcoords);
335: vgid[vcount] = (k*npts+j)*npts+i+1;
336: }
337: }
338: }
340: /* create ghosted vertices requested by user - below the current plane */
341: if (ise[2*dim-2] > 0) {
342: for (k = (dim==3?ise[4]-nghost:ise[4]); k <= (dim==3?ise[4]-1:ise[5]); k++) {
343: for (j = (dim==2?ise[2]-nghost:ise[2]); j <= (dim==2?ise[2]-1:ise[3]); j++) {
344: for (i = (dim>1?ise[0]:ise[0]-nghost); i <= (dim>1?ise[1]:ise[0]-1); i++, vcount++) {
345: DMMoab_SetStructuredCoords_Private(i,j,k,hx,hy,hz,vcount,vcoords);
346: vgid[vcount] = (k*npts+j)*npts+i+1;
347: }
348: }
349: }
350: }
352: /* create ghosted vertices requested by user - above the current plane */
353: if (ise[2*dim-1] < nele) {
354: for (k = (dim==3?ise[5]+1:ise[4]); k <= (dim==3?ise[5]+nghost:ise[5]); k++) {
355: for (j = (dim==2?ise[3]+1:ise[2]); j <= (dim==2?ise[3]+nghost:ise[3]); j++) {
356: for (i = (dim>1?ise[0]:ise[1]+1); i <= (dim>1?ise[1]:ise[1]+nghost); i++, vcount++) {
357: DMMoab_SetStructuredCoords_Private(i,j,k,hx,hy,hz,vcount,vcoords);
358: vgid[vcount] = (k*npts+j)*npts+i+1;
359: }
360: }
361: }
362: }
364: if (locnpts+ghnpts != vcount) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of vertices! (%D!=%D)",locnpts+ghnpts,vcount);
366: merr = mbiface->get_entities_by_type(0,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr);
368: /* The global ID tag is applied to each owned
369: vertex. It acts as an global identifier which MOAB uses to
370: assemble the individual pieces of the mesh:
371: Set the global ID indices */
372: merr = mbiface->tag_set_data(id_tag,ownedvtx,vgid.data());MBERRNM(merr);
373:
374: /* Create elements between mesh points using the ReadUtilInterface
375: get the reference to element connectivities for all local elements from the ReadUtilInterface */
376: merr = readMeshIface->get_element_connect (locnele+ghnele,vpere,etype,1,efirst,connectivity);MBERRNM(merr);
378: /* offset appropriately so that only local ID and not global ID numbers are set for connectivity array */
379: vfirst-=vgid[0]-1;
381: /* 3. Loop over elements in 3 nested loops over i, j, k; for each (i,j,k):
382: and then set the connectivity for each element appropriately */
383: int ecount=0;
385: /* create ghosted elements requested by user - below the current plane */
386: if (ise[2*dim-2] >= nghost) {
387: for (k = (dim==3?ise[4]-nghost:ise[4]); k < (dim==3?ise[4]:std::max(ise[5],1)); k++) {
388: for (j = (dim==2?ise[2]-nghost:ise[2]); j < (dim==2?ise[2]:std::max(ise[3],1)); j++) {
389: for (i = (dim>1?ise[0]:ise[0]-nghost); i < (dim>1?std::max(ise[1],1):ise[0]); i++) {
390: DMMoab_SetElementConnectivity_Private(useSimplex, dim, etype, &ecount, vpere, nele, i, j, k, vfirst, connectivity);
391: }
392: }
393: }
394: }
396: /* create owned elements requested by user */
397: for (k = ise[4]; k < std::max(ise[5],1); k++) {
398: for (j = ise[2]; j < std::max(ise[3],1); j++) {
399: for (i = ise[0]; i < std::max(ise[1],1); i++) {
400: DMMoab_SetElementConnectivity_Private(useSimplex, dim, etype, &ecount, vpere, nele, i, j, k, vfirst, connectivity);
401: }
402: }
403: }
405: /* create ghosted elements requested by user - above the current plane */
406: if (ise[2*dim-1] <= nele-nghost) {
407: for (k = (dim==3?ise[5]:ise[4]); k < (dim==3?ise[5]+nghost:std::max(ise[5],1)); k++) {
408: for (j = (dim==2?ise[3]:ise[2]); j < (dim==2?ise[3]+nghost:std::max(ise[3],1)); j++) {
409: for (i = (dim>1?ise[0]:ise[1]); i < (dim>1?std::max(ise[1],1):ise[1]+nghost); i++) {
410: DMMoab_SetElementConnectivity_Private(useSimplex, dim, etype, &ecount, vpere, nele, i, j, k, vfirst, connectivity);
411: }
412: }
413: }
414: }
416: merr = readMeshIface->update_adjacencies(efirst,locnele+ghnele,vpere,connectivity);MBERRNM(merr);
417:
418: /* 2. Get the vertices and hexes from moab and check their numbers against I*J*K and (I-1)*(J-1)*(K-1), resp.
419: first '0' specifies "root set", or entire MOAB instance, second the entity dimension being requested */
420: merr = mbiface->get_entities_by_dimension(0, dim, ownedelms);MBERRNM(merr);
422: if (locnele+ghnele != (int) ownedelms.size()) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of elements! (%D!=%D)",locnele+ghnele,ownedelms.size());
423: else if(locnpts+ghnpts != (int) ownedvtx.size()) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Created the wrong number of vertices! (%D!=%D)",locnpts+ghnpts,ownedvtx.size());
424: else {
425: PetscInfo2(NULL, "Created %D elements and %D vertices.\n", ownedelms.size(), ownedvtx.size());
426: }
427:
428: /* lets create some sets */
429: merr = mbiface->tag_get_handle(GEOM_DIMENSION_TAG_NAME, 1, moab::MB_TYPE_INTEGER, geom_tag, moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT);MBERRNM(merr);
431: merr = mbiface->create_meshset(moab::MESHSET_SET, regionset);MBERRNM(merr);
432: merr = mbiface->add_entities(regionset, ownedelms);MBERRNM(merr);
433: merr = mbiface->tag_set_data(geom_tag, ®ionset, 1, &dmmoab->dim);MBERRNM(merr);
434: merr = mbiface->add_parent_child(dmmoab->fileset,regionset);MBERRNM(merr);
435: merr = mbiface->unite_meshset(dmmoab->fileset, regionset);MBERRNM(merr);
437: merr = mbiface->create_meshset(moab::MESHSET_SET, vtxset);MBERRNM(merr);
438: merr = mbiface->add_entities(vtxset, ownedvtx);MBERRNM(merr);
439: merr = mbiface->add_parent_child(dmmoab->fileset,vtxset);MBERRNM(merr);
440: merr = mbiface->unite_meshset(dmmoab->fileset, vtxset);MBERRNM(merr);
442: if (build_adjacencies) {
443: // generate all lower dimensional adjacencies
444: merr = mbiface->get_adjacencies( ownedelms, dim-1, true, adj, moab::Interface::UNION );MBERRNM(merr);
445: merr = dmmoab->pcomm->get_part_entities(dim2, dim-1);MBERRNM(merr);
446: adj.merge(dim2);
448: /* create face sets */
449: merr = mbiface->create_meshset(moab::MESHSET_SET, faceset);MBERRNM(merr);
450: merr = mbiface->add_entities(faceset, adj);MBERRNM(merr);
451: merr = mbiface->add_parent_child(dmmoab->fileset,faceset);MBERRNM(merr);
452: i=dim-1;
453: merr = mbiface->tag_set_data(geom_tag, &faceset, 1, &i);MBERRNM(merr);
454: merr = mbiface->unite_meshset(dmmoab->fileset, faceset);MBERRNM(merr);
455: PetscInfo2(NULL, "Found %d %d-Dim quantities.\n", adj.size(), dim-1);
457: if (dim > 2) {
458: dim2.clear();
459: /* create edge sets, if appropriate i.e., if dim=3 */
460: merr = mbiface->create_meshset(moab::MESHSET_SET, edgeset);MBERRNM(merr);
461: merr = mbiface->get_adjacencies(adj, dim-1, true, dim2, moab::Interface::UNION );MBERRNM(merr);
462: merr = mbiface->add_entities(edgeset, dim2);MBERRNM(merr);
463: merr = mbiface->add_parent_child(dmmoab->fileset,edgeset);MBERRNM(merr);
464: i=dim-2;
465: merr = mbiface->tag_set_data(geom_tag, &edgeset, 1, &i);MBERRNM(merr);
466: merr = mbiface->unite_meshset(dmmoab->fileset, edgeset);MBERRNM(merr);
467: PetscInfo2(NULL, "Found %d %d-Dim quantities.\n", adj.size(), dim-2);
468: }
469: }
471: /* check the handles */
472: merr = pcomm->check_all_shared_handles();MBERRV(mbiface,merr);
474: /* resolve the shared entities by exchanging information to adjacent processors */
475: merr = mbiface->get_entities_by_type(dmmoab->fileset,etype,ownedelms,true);MBERRNM(merr);
476: merr = pcomm->resolve_shared_ents(dmmoab->fileset,ownedelms,dim,dim-1,NULL,&id_tag);MBERRV(mbiface,merr);
478: merr = pcomm->exchange_ghost_cells(dim,0,user_nghost,dim,true,false,&dmmoab->fileset);MBERRV(mbiface,merr);
480: /* Reassign global IDs on all entities. */
481: merr = pcomm->assign_global_ids(dmmoab->fileset,dim,1,false,true,false);MBERRNM(merr);
483: /* Everything is set up, now just do a tag exchange to update tags
484: on all of the ghost vertexes */
485: merr = mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,ownedvtx,true);MBERRNM(merr);
486: merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, ownedelms);MBERRNM(merr);
487:
488: merr = pcomm->exchange_tags(id_tag,ownedvtx);MBERRV(mbiface,merr);
489: merr = pcomm->exchange_tags(id_tag,ownedelms);MBERRV(mbiface,merr);
490: return(0);
491: }
496: PetscErrorCode DMMoab_GetReadOptions_Private(PetscBool by_rank, PetscInt numproc, PetscInt dim, MoabReadMode mode, PetscInt dbglevel, const char* dm_opts, const char* extra_opts, const char** read_opts)
497: {
498: char *ropts;
499: char ropts_par[PETSC_MAX_PATH_LEN];
500: char ropts_dbg[PETSC_MAX_PATH_LEN];
504: PetscMalloc(PETSC_MAX_PATH_LEN,&ropts);
505: PetscMemzero(&ropts_par,PETSC_MAX_PATH_LEN);
506: PetscMemzero(&ropts_dbg,PETSC_MAX_PATH_LEN);
508: /* do parallel read unless using only one processor */
509: if (numproc > 1) {
510: PetscSNPrintf(ropts_par, PETSC_MAX_PATH_LEN, "PARALLEL=%s;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=%d.0.1%s;",MoabReadModes[mode],dim,(by_rank ? ";PARTITION_BY_RANK":""));
511: }
513: if (dbglevel) {
514: if (numproc>1) {
515: PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "%sCPUTIME;DEBUG_IO=%d;DEBUG_PIO=%d;",dbglevel,dbglevel);
516: }
517: else {
518: PetscSNPrintf(ropts_dbg, PETSC_MAX_PATH_LEN, "%sCPUTIME;DEBUG_IO=%d;",dbglevel);
519: }
520: }
522: PetscSNPrintf(ropts, PETSC_MAX_PATH_LEN, "%s%s%s%s",ropts_par,ropts_dbg,(extra_opts?extra_opts:""),(dm_opts?dm_opts:""));
523: *read_opts = ropts;
524: return(0);
525: }
530: /*@
531: DMMoabLoadFromFile - Creates a DM object by loading the mesh from a user specified file.
533: Collective on MPI_Comm
535: Input Parameters:
536: + comm - The communicator for the DM object
537: . dim - The spatial dimension
538: . filename - The name of the mesh file to be loaded
539: . usrreadopts - The options string to read a MOAB mesh.
540: Reference (Parallel Mesh Initialization: http://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html#fivetwo)
542: Output Parameter:
543: . dm - The DM object
545: Level: beginner
547: .keywords: DM, create
549: .seealso: DMSetType(), DMCreate(), DMMoabCreateBoxMesh()
550: @*/
551: PetscErrorCode DMMoabLoadFromFile(MPI_Comm comm,PetscInt dim,const char* filename, const char* usrreadopts, DM *dm)
552: {
553: moab::ErrorCode merr;
554: PetscInt nprocs;
555: DM_Moab *dmmoab;
556: moab::Interface *mbiface;
557: moab::ParallelComm *pcomm;
558: moab::Range verts,elems;
559: const char *readopts;
565: /* Create the basic DMMoab object and keep the default parameters created by DM impls */
566: DMMoabCreateMoab(comm, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, dm);
568: /* get all the necessary handles from the private DM object */
569: dmmoab = (DM_Moab*)(*dm)->data;
570: mbiface = dmmoab->mbiface;
571: pcomm = dmmoab->pcomm;
572: nprocs = pcomm->size();
573: /* TODO: Decipher dimension based on the loaded mesh instead of getting from user */
574: dmmoab->dim = dim;
576: /* create a file set to associate all entities in current mesh */
577: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, dmmoab->fileset);MBERR("Creating file set failed", merr);
579: /* add mesh loading options specific to the DM */
580: DMMoab_GetReadOptions_Private(dmmoab->partition_by_rank, nprocs, dim, dmmoab->read_mode,
581: dmmoab->rw_dbglevel, dmmoab->extra_read_options, usrreadopts, &readopts);
583: PetscInfo2(*dm, "Reading file %s with options: %s\n",filename,readopts);
585: /* Load the mesh from a file. */
586: merr = mbiface->load_file(filename, &dmmoab->fileset, readopts);MBERRVM(mbiface,"Reading MOAB file failed.", merr);
588: /* Reassign global IDs on all entities. */
589: merr = pcomm->assign_global_ids(dmmoab->fileset,dim,1,true,true,true);MBERRNM(merr);
591: /* load the local vertices */
592: merr = mbiface->get_entities_by_type(dmmoab->fileset, moab::MBVERTEX, verts, true);MBERRNM(merr);
593: /* load the local elements */
594: merr = mbiface->get_entities_by_dimension(dmmoab->fileset, dim, elems, true);MBERRNM(merr);
596: /* Everything is set up, now just do a tag exchange to update tags
597: on all of the ghost vertexes */
598: merr = pcomm->exchange_tags(dmmoab->ltog_tag,verts);MBERRV(mbiface,merr);
599: merr = pcomm->exchange_tags(dmmoab->ltog_tag,elems);MBERRV(mbiface,merr);
601: merr = pcomm->exchange_ghost_cells(dim,0,1,0,true,true,&dmmoab->fileset);MBERRV(mbiface,merr);
603: merr = pcomm->collective_sync_partition();MBERR("Collective sync failed", merr);
605: PetscInfo3(*dm, "MOAB file '%s' was successfully loaded. Found %D vertices and %D elements.\n", filename, verts.size(), elems.size());
606: PetscFree(readopts);
607: return(0);
608: }