Actual source code: dacreate.c
2: #include <petsc/private/dmdaimpl.h>
4: PetscErrorCode DMSetFromOptions_DA(PetscOptionItems *PetscOptionsObject,DM da)
5: {
6: DM_DA *dd = (DM_DA*)da->data;
7: PetscInt refine = 0,dim = da->dim,maxnlevels = 100,refx[100],refy[100],refz[100],n,i;
8: PetscBool flg;
16: PetscOptionsHead(PetscOptionsObject,"DMDA Options");
17: PetscOptionsBoundedInt("-da_grid_x","Number of grid points in x direction","DMDASetSizes",dd->M,&dd->M,NULL,1);
18: if (dim > 1) PetscOptionsBoundedInt("-da_grid_y","Number of grid points in y direction","DMDASetSizes",dd->N,&dd->N,NULL,1);
19: if (dim > 2) PetscOptionsBoundedInt("-da_grid_z","Number of grid points in z direction","DMDASetSizes",dd->P,&dd->P,NULL,1);
21: PetscOptionsBoundedInt("-da_overlap","Decomposition overlap in all directions","DMDASetOverlap",dd->xol,&dd->xol,&flg,0);
22: if (flg) DMDASetOverlap(da,dd->xol,dd->xol,dd->xol);
23: PetscOptionsBoundedInt("-da_overlap_x","Decomposition overlap in x direction","DMDASetOverlap",dd->xol,&dd->xol,NULL,0);
24: if (dim > 1) PetscOptionsBoundedInt("-da_overlap_y","Decomposition overlap in y direction","DMDASetOverlap",dd->yol,&dd->yol,NULL,0);
25: if (dim > 2) PetscOptionsBoundedInt("-da_overlap_z","Decomposition overlap in z direction","DMDASetOverlap",dd->zol,&dd->zol,NULL,0);
27: PetscOptionsBoundedInt("-da_local_subdomains","","DMDASetNumLocalSubdomains",dd->Nsub,&dd->Nsub,&flg,PETSC_DECIDE);
28: if (flg) DMDASetNumLocalSubDomains(da,dd->Nsub);
30: /* Handle DMDA parallel distribution */
31: PetscOptionsBoundedInt("-da_processors_x","Number of processors in x direction","DMDASetNumProcs",dd->m,&dd->m,NULL,PETSC_DECIDE);
32: if (dim > 1) PetscOptionsBoundedInt("-da_processors_y","Number of processors in y direction","DMDASetNumProcs",dd->n,&dd->n,NULL,PETSC_DECIDE);
33: if (dim > 2) PetscOptionsBoundedInt("-da_processors_z","Number of processors in z direction","DMDASetNumProcs",dd->p,&dd->p,NULL,PETSC_DECIDE);
34: /* Handle DMDA refinement */
35: PetscOptionsBoundedInt("-da_refine_x","Refinement ratio in x direction","DMDASetRefinementFactor",dd->refine_x,&dd->refine_x,NULL,1);
36: if (dim > 1) PetscOptionsBoundedInt("-da_refine_y","Refinement ratio in y direction","DMDASetRefinementFactor",dd->refine_y,&dd->refine_y,NULL,1);
37: if (dim > 2) PetscOptionsBoundedInt("-da_refine_z","Refinement ratio in z direction","DMDASetRefinementFactor",dd->refine_z,&dd->refine_z,NULL,1);
38: dd->coarsen_x = dd->refine_x; dd->coarsen_y = dd->refine_y; dd->coarsen_z = dd->refine_z;
40: /* Get refinement factors, defaults taken from the coarse DMDA */
41: DMDAGetRefinementFactor(da,&refx[0],&refy[0],&refz[0]);
42: for (i=1; i<maxnlevels; i++) {
43: refx[i] = refx[0];
44: refy[i] = refy[0];
45: refz[i] = refz[0];
46: }
47: n = maxnlevels;
48: PetscOptionsIntArray("-da_refine_hierarchy_x","Refinement factor for each level","None",refx,&n,&flg);
49: if (flg) {
50: dd->refine_x = refx[0];
51: dd->refine_x_hier_n = n;
52: PetscMalloc1(n,&dd->refine_x_hier);
53: PetscArraycpy(dd->refine_x_hier,refx,n);
54: }
55: if (dim > 1) {
56: n = maxnlevels;
57: PetscOptionsIntArray("-da_refine_hierarchy_y","Refinement factor for each level","None",refy,&n,&flg);
58: if (flg) {
59: dd->refine_y = refy[0];
60: dd->refine_y_hier_n = n;
61: PetscMalloc1(n,&dd->refine_y_hier);
62: PetscArraycpy(dd->refine_y_hier,refy,n);
63: }
64: }
65: if (dim > 2) {
66: n = maxnlevels;
67: PetscOptionsIntArray("-da_refine_hierarchy_z","Refinement factor for each level","None",refz,&n,&flg);
68: if (flg) {
69: dd->refine_z = refz[0];
70: dd->refine_z_hier_n = n;
71: PetscMalloc1(n,&dd->refine_z_hier);
72: PetscArraycpy(dd->refine_z_hier,refz,n);
73: }
74: }
76: PetscOptionsBoundedInt("-da_refine","Uniformly refine DA one or more times","None",refine,&refine,NULL,0);
77: PetscOptionsTail();
79: while (refine--) {
80: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
81: dd->M = dd->refine_x*dd->M;
82: } else {
83: dd->M = 1 + dd->refine_x*(dd->M - 1);
84: }
85: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
86: dd->N = dd->refine_y*dd->N;
87: } else {
88: dd->N = 1 + dd->refine_y*(dd->N - 1);
89: }
90: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
91: dd->P = dd->refine_z*dd->P;
92: } else {
93: dd->P = 1 + dd->refine_z*(dd->P - 1);
94: }
95: da->levelup++;
96: if (da->levelup - da->leveldown >= 0) {
97: dd->refine_x = refx[da->levelup - da->leveldown];
98: dd->refine_y = refy[da->levelup - da->leveldown];
99: dd->refine_z = refz[da->levelup - da->leveldown];
100: }
101: if (da->levelup - da->leveldown >= 1) {
102: dd->coarsen_x = refx[da->levelup - da->leveldown - 1];
103: dd->coarsen_y = refy[da->levelup - da->leveldown - 1];
104: dd->coarsen_z = refz[da->levelup - da->leveldown - 1];
105: }
106: }
107: return 0;
108: }
110: extern PetscErrorCode DMCreateGlobalVector_DA(DM,Vec*);
111: extern PetscErrorCode DMCreateLocalVector_DA(DM,Vec*);
112: extern PetscErrorCode DMGlobalToLocalBegin_DA(DM,Vec,InsertMode,Vec);
113: extern PetscErrorCode DMGlobalToLocalEnd_DA(DM,Vec,InsertMode,Vec);
114: extern PetscErrorCode DMLocalToGlobalBegin_DA(DM,Vec,InsertMode,Vec);
115: extern PetscErrorCode DMLocalToGlobalEnd_DA(DM,Vec,InsertMode,Vec);
116: extern PetscErrorCode DMLocalToLocalBegin_DA(DM,Vec,InsertMode,Vec);
117: extern PetscErrorCode DMLocalToLocalEnd_DA(DM,Vec,InsertMode,Vec);
118: extern PetscErrorCode DMCreateInterpolation_DA(DM,DM,Mat*,Vec*);
119: extern PetscErrorCode DMCreateColoring_DA(DM,ISColoringType,ISColoring*);
120: extern PetscErrorCode DMCreateMatrix_DA(DM,Mat*);
121: extern PetscErrorCode DMCreateCoordinateDM_DA(DM,DM*);
122: extern PetscErrorCode DMRefine_DA(DM,MPI_Comm,DM*);
123: extern PetscErrorCode DMCoarsen_DA(DM,MPI_Comm,DM*);
124: extern PetscErrorCode DMRefineHierarchy_DA(DM,PetscInt,DM[]);
125: extern PetscErrorCode DMCoarsenHierarchy_DA(DM,PetscInt,DM[]);
126: extern PetscErrorCode DMCreateInjection_DA(DM,DM,Mat*);
127: extern PetscErrorCode DMView_DA(DM,PetscViewer);
128: extern PetscErrorCode DMSetUp_DA(DM);
129: extern PetscErrorCode DMDestroy_DA(DM);
130: extern PetscErrorCode DMCreateDomainDecomposition_DA(DM,PetscInt*,char***,IS**,IS**,DM**);
131: extern PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM,PetscInt,DM*,VecScatter**,VecScatter**,VecScatter**);
132: PETSC_INTERN PetscErrorCode DMGetCompatibility_DA(DM,DM,PetscBool*,PetscBool*);
134: PetscErrorCode DMLoad_DA(DM da,PetscViewer viewer)
135: {
136: PetscInt dim,m,n,p,dof,swidth;
137: DMDAStencilType stencil;
138: DMBoundaryType bx,by,bz;
139: PetscBool coors;
140: DM dac;
141: Vec c;
143: PetscViewerBinaryRead(viewer,&dim,1,NULL,PETSC_INT);
144: PetscViewerBinaryRead(viewer,&m,1,NULL,PETSC_INT);
145: PetscViewerBinaryRead(viewer,&n,1,NULL,PETSC_INT);
146: PetscViewerBinaryRead(viewer,&p,1,NULL,PETSC_INT);
147: PetscViewerBinaryRead(viewer,&dof,1,NULL,PETSC_INT);
148: PetscViewerBinaryRead(viewer,&swidth,1,NULL,PETSC_INT);
149: PetscViewerBinaryRead(viewer,&bx,1,NULL,PETSC_ENUM);
150: PetscViewerBinaryRead(viewer,&by,1,NULL,PETSC_ENUM);
151: PetscViewerBinaryRead(viewer,&bz,1,NULL,PETSC_ENUM);
152: PetscViewerBinaryRead(viewer,&stencil,1,NULL,PETSC_ENUM);
154: DMSetDimension(da, dim);
155: DMDASetSizes(da, m,n,p);
156: DMDASetBoundaryType(da, bx, by, bz);
157: DMDASetDof(da, dof);
158: DMDASetStencilType(da, stencil);
159: DMDASetStencilWidth(da, swidth);
160: DMSetUp(da);
161: PetscViewerBinaryRead(viewer,&coors,1,NULL,PETSC_ENUM);
162: if (coors) {
163: DMGetCoordinateDM(da,&dac);
164: DMCreateGlobalVector(dac,&c);
165: VecLoad(c,viewer);
166: DMSetCoordinates(da,c);
167: VecDestroy(&c);
168: }
169: return 0;
170: }
172: PetscErrorCode DMCreateSubDM_DA(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
173: {
174: DM_DA *da = (DM_DA*) dm->data;
176: if (subdm) {
177: PetscSF sf;
178: Vec coords;
179: void *ctx;
180: /* Cannot use DMClone since the dof stuff is mixed in. Ugh
181: DMClone(dm, subdm); */
182: DMCreate(PetscObjectComm((PetscObject)dm), subdm);
183: DMGetPointSF(dm, &sf);
184: DMSetPointSF(*subdm, sf);
185: DMGetApplicationContext(dm, &ctx);
186: DMSetApplicationContext(*subdm, ctx);
187: DMGetCoordinatesLocal(dm, &coords);
188: if (coords) {
189: DMSetCoordinatesLocal(*subdm, coords);
190: } else {
191: DMGetCoordinates(dm, &coords);
192: if (coords) DMSetCoordinates(*subdm, coords);
193: }
195: DMSetType(*subdm, DMDA);
196: DMSetDimension(*subdm, dm->dim);
197: DMDASetSizes(*subdm, da->M, da->N, da->P);
198: DMDASetNumProcs(*subdm, da->m, da->n, da->p);
199: DMDASetBoundaryType(*subdm, da->bx, da->by, da->bz);
200: DMDASetDof(*subdm, numFields);
201: DMDASetStencilType(*subdm, da->stencil_type);
202: DMDASetStencilWidth(*subdm, da->s);
203: DMDASetOwnershipRanges(*subdm, da->lx, da->ly, da->lz);
204: }
205: if (is) {
206: PetscInt *indices, cnt = 0, dof = da->w, i, j;
208: PetscMalloc1(da->Nlocal*numFields/dof, &indices);
209: for (i = da->base/dof; i < (da->base+da->Nlocal)/dof; ++i) {
210: for (j = 0; j < numFields; ++j) {
211: indices[cnt++] = dof*i + fields[j];
212: }
213: }
215: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cnt, indices, PETSC_OWN_POINTER, is);
216: }
217: return 0;
218: }
220: PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
221: {
222: PetscInt i;
223: DM_DA *dd = (DM_DA*)dm->data;
224: PetscInt dof = dd->w;
226: if (len) *len = dof;
227: if (islist) {
228: Vec v;
229: PetscInt rstart,n;
231: DMGetGlobalVector(dm,&v);
232: VecGetOwnershipRange(v,&rstart,NULL);
233: VecGetLocalSize(v,&n);
234: DMRestoreGlobalVector(dm,&v);
235: PetscMalloc1(dof,islist);
236: for (i=0; i<dof; i++) {
237: ISCreateStride(PetscObjectComm((PetscObject)dm),n/dof,rstart+i,dof,&(*islist)[i]);
238: }
239: }
240: if (namelist) {
241: PetscMalloc1(dof, namelist);
242: if (dd->fieldname) {
243: for (i=0; i<dof; i++) {
244: PetscStrallocpy(dd->fieldname[i],&(*namelist)[i]);
245: }
246: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently DMDA must have fieldnames");
247: }
248: if (dmlist) {
249: DM da;
251: DMDACreate(PetscObjectComm((PetscObject)dm), &da);
252: DMSetDimension(da, dm->dim);
253: DMDASetSizes(da, dd->M, dd->N, dd->P);
254: DMDASetNumProcs(da, dd->m, dd->n, dd->p);
255: DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);
256: DMDASetDof(da, 1);
257: DMDASetStencilType(da, dd->stencil_type);
258: DMDASetStencilWidth(da, dd->s);
259: DMSetUp(da);
260: PetscMalloc1(dof,dmlist);
261: for (i=0; i<dof-1; i++) PetscObjectReference((PetscObject)da);
262: for (i=0; i<dof; i++) (*dmlist)[i] = da;
263: }
264: return 0;
265: }
267: PetscErrorCode DMClone_DA(DM dm, DM *newdm)
268: {
269: DM_DA *da = (DM_DA *) dm->data;
271: DMSetType(*newdm, DMDA);
272: DMSetDimension(*newdm, dm->dim);
273: DMDASetSizes(*newdm, da->M, da->N, da->P);
274: DMDASetNumProcs(*newdm, da->m, da->n, da->p);
275: DMDASetBoundaryType(*newdm, da->bx, da->by, da->bz);
276: DMDASetDof(*newdm, da->w);
277: DMDASetStencilType(*newdm, da->stencil_type);
278: DMDASetStencilWidth(*newdm, da->s);
279: DMDASetOwnershipRanges(*newdm, da->lx, da->ly, da->lz);
280: DMSetUp(*newdm);
281: return 0;
282: }
284: static PetscErrorCode DMHasCreateInjection_DA(DM dm, PetscBool *flg)
285: {
286: DM_DA *da = (DM_DA *)dm->data;
290: *flg = da->interptype == DMDA_Q1 ? PETSC_TRUE : PETSC_FALSE;
291: return 0;
292: }
294: static PetscErrorCode DMGetDimPoints_DA(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
295: {
296: DMDAGetDepthStratum(dm, dim, pStart, pEnd);
297: return 0;
298: }
300: static PetscErrorCode DMGetNeighbors_DA(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[])
301: {
302: PetscInt dim;
303: DMDAStencilType st;
305: DMDAGetNeighbors(dm,ranks);
306: DMDAGetInfo(dm,&dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&st);
308: switch (dim) {
309: case 1:
310: *nranks = 3;
311: /* if (st == DMDA_STENCIL_STAR) { *nranks = 3; } */
312: break;
313: case 2:
314: *nranks = 9;
315: /* if (st == DMDA_STENCIL_STAR) { *nranks = 5; } */
316: break;
317: case 3:
318: *nranks = 27;
319: /* if (st == DMDA_STENCIL_STAR) { *nranks = 7; } */
320: break;
321: default:
322: break;
323: }
324: return 0;
325: }
327: /*MC
328: DMDA = "da" - A DM object that is used to manage data for a structured grid in 1, 2, or 3 dimensions.
329: In the global representation of the vector each process stores a non-overlapping rectangular (or slab in 3d) portion of the grid points.
330: In the local representation these rectangular regions (slabs) are extended in all directions by a stencil width.
332: The vectors can be thought of as either cell centered or vertex centered on the mesh. But some variables cannot be cell centered and others
333: vertex centered; see the documentation for DMSTAG, a similar DM implementation which supports these staggered grids.
335: Level: intermediate
337: .seealso: DMType, DMCOMPOSITE, DMSTAG, DMDACreate(), DMCreate(), DMSetType()
338: M*/
340: extern PetscErrorCode DMLocatePoints_DA_Regular(DM,Vec,DMPointLocationType,PetscSF);
341: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_DMDA(PetscObject,PetscViewer);
343: PETSC_EXTERN PetscErrorCode DMCreate_DA(DM da)
344: {
345: DM_DA *dd;
348: PetscNewLog(da,&dd);
349: da->data = dd;
351: da->dim = -1;
352: dd->interptype = DMDA_Q1;
353: dd->refine_x = 2;
354: dd->refine_y = 2;
355: dd->refine_z = 2;
356: dd->coarsen_x = 2;
357: dd->coarsen_y = 2;
358: dd->coarsen_z = 2;
359: dd->fieldname = NULL;
360: dd->nlocal = -1;
361: dd->Nlocal = -1;
362: dd->M = -1;
363: dd->N = -1;
364: dd->P = -1;
365: dd->m = -1;
366: dd->n = -1;
367: dd->p = -1;
368: dd->w = -1;
369: dd->s = -1;
371: dd->xs = -1; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1;
372: dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1;
374: dd->Nsub = 1;
375: dd->xol = 0;
376: dd->yol = 0;
377: dd->zol = 0;
378: dd->xo = 0;
379: dd->yo = 0;
380: dd->zo = 0;
381: dd->Mo = -1;
382: dd->No = -1;
383: dd->Po = -1;
385: dd->gtol = NULL;
386: dd->ltol = NULL;
387: dd->ao = NULL;
388: PetscStrallocpy(AOBASIC,(char**)&dd->aotype);
389: dd->base = -1;
390: dd->bx = DM_BOUNDARY_NONE;
391: dd->by = DM_BOUNDARY_NONE;
392: dd->bz = DM_BOUNDARY_NONE;
393: dd->stencil_type = DMDA_STENCIL_BOX;
394: dd->interptype = DMDA_Q1;
395: dd->lx = NULL;
396: dd->ly = NULL;
397: dd->lz = NULL;
399: dd->elementtype = DMDA_ELEMENT_Q1;
401: da->ops->globaltolocalbegin = DMGlobalToLocalBegin_DA;
402: da->ops->globaltolocalend = DMGlobalToLocalEnd_DA;
403: da->ops->localtoglobalbegin = DMLocalToGlobalBegin_DA;
404: da->ops->localtoglobalend = DMLocalToGlobalEnd_DA;
405: da->ops->localtolocalbegin = DMLocalToLocalBegin_DA;
406: da->ops->localtolocalend = DMLocalToLocalEnd_DA;
407: da->ops->createglobalvector = DMCreateGlobalVector_DA;
408: da->ops->createlocalvector = DMCreateLocalVector_DA;
409: da->ops->createinterpolation = DMCreateInterpolation_DA;
410: da->ops->getcoloring = DMCreateColoring_DA;
411: da->ops->creatematrix = DMCreateMatrix_DA;
412: da->ops->refine = DMRefine_DA;
413: da->ops->coarsen = DMCoarsen_DA;
414: da->ops->refinehierarchy = DMRefineHierarchy_DA;
415: da->ops->coarsenhierarchy = DMCoarsenHierarchy_DA;
416: da->ops->createinjection = DMCreateInjection_DA;
417: da->ops->hascreateinjection = DMHasCreateInjection_DA;
418: da->ops->destroy = DMDestroy_DA;
419: da->ops->view = NULL;
420: da->ops->setfromoptions = DMSetFromOptions_DA;
421: da->ops->setup = DMSetUp_DA;
422: da->ops->clone = DMClone_DA;
423: da->ops->load = DMLoad_DA;
424: da->ops->createcoordinatedm = DMCreateCoordinateDM_DA;
425: da->ops->createsubdm = DMCreateSubDM_DA;
426: da->ops->createfielddecomposition = DMCreateFieldDecomposition_DA;
427: da->ops->createdomaindecomposition = DMCreateDomainDecomposition_DA;
428: da->ops->createddscatters = DMCreateDomainDecompositionScatters_DA;
429: da->ops->getdimpoints = DMGetDimPoints_DA;
430: da->ops->getneighbors = DMGetNeighbors_DA;
431: da->ops->locatepoints = DMLocatePoints_DA_Regular;
432: da->ops->getcompatibility = DMGetCompatibility_DA;
433: PetscObjectComposeFunction((PetscObject)da,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_DMDA);
434: return 0;
435: }
437: /*@
438: DMDACreate - Creates a DMDA object.
440: Collective
442: Input Parameter:
443: . comm - The communicator for the DMDA object
445: Output Parameter:
446: . da - The DMDA object
448: Level: advanced
450: Developers Note:
451: Since there exists DMDACreate1/2/3d() should this routine even exist?
453: .seealso: DMDASetSizes(), DMClone(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
454: @*/
455: PetscErrorCode DMDACreate(MPI_Comm comm, DM *da)
456: {
458: DMCreate(comm,da);
459: DMSetType(*da,DMDA);
460: return 0;
461: }