Actual source code: dacreate.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/dmdaimpl.h> /*I "petscdmda.h" I*/
6: PetscErrorCode DMSetFromOptions_DA(DM da)
7: {
9: DM_DA *dd = (DM_DA*)da->data;
10: PetscInt refine = 0,maxnlevels = 100,*refx,*refy,*refz,n,i;
11: PetscBool negativeMNP = PETSC_FALSE,bM = PETSC_FALSE,bN = PETSC_FALSE, bP = PETSC_FALSE,flg;
16: if (dd->M < 0) {
17: dd->M = -dd->M;
18: bM = PETSC_TRUE;
19: negativeMNP = PETSC_TRUE;
20: }
21: if (dd->dim > 1 && dd->N < 0) {
22: dd->N = -dd->N;
23: bN = PETSC_TRUE;
24: negativeMNP = PETSC_TRUE;
25: }
26: if (dd->dim > 2 && dd->P < 0) {
27: dd->P = -dd->P;
28: bP = PETSC_TRUE;
29: negativeMNP = PETSC_TRUE;
30: }
32: PetscOptionsHead("DMDA Options");
33: if (bM) {PetscOptionsInt("-da_grid_x","Number of grid points in x direction","DMDASetSizes",dd->M,&dd->M,NULL);}
34: if (bN) {PetscOptionsInt("-da_grid_y","Number of grid points in y direction","DMDASetSizes",dd->N,&dd->N,NULL);}
35: if (bP) {PetscOptionsInt("-da_grid_z","Number of grid points in z direction","DMDASetSizes",dd->P,&dd->P,NULL);}
37: PetscOptionsInt("-da_overlap","Decomposition overlap in all directions","DMDASetOverlap",dd->xol,&dd->xol,&flg);
38: if (flg) {DMDASetOverlap(da,dd->xol,dd->xol,dd->xol);}
39: PetscOptionsInt("-da_overlap_x","Decomposition overlap in x direction","DMDASetOverlap",dd->xol,&dd->xol,NULL);
40: if (dd->dim > 1) {PetscOptionsInt("-da_overlap_y","Decomposition overlap in y direction","DMDASetOverlap",dd->yol,&dd->yol,NULL);}
41: if (dd->dim > 2) {PetscOptionsInt("-da_overlap_z","Decomposition overlap in z direction","DMDASetOverlap",dd->zol,&dd->zol,NULL);}
43: PetscOptionsInt("-da_local_subdomains","","DMDASetNumLocalSubdomains",dd->Nsub,&dd->Nsub,&flg);
44: if (flg) {DMDASetNumLocalSubDomains(da,dd->Nsub);}
46: /* Handle DMDA parallel distibution */
47: PetscOptionsInt("-da_processors_x","Number of processors in x direction","DMDASetNumProcs",dd->m,&dd->m,NULL);
48: if (dd->dim > 1) {PetscOptionsInt("-da_processors_y","Number of processors in y direction","DMDASetNumProcs",dd->n,&dd->n,NULL);}
49: if (dd->dim > 2) {PetscOptionsInt("-da_processors_z","Number of processors in z direction","DMDASetNumProcs",dd->p,&dd->p,NULL);}
50: /* Handle DMDA refinement */
51: PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DMDASetRefinementFactor",dd->refine_x,&dd->refine_x,NULL);
52: if (dd->dim > 1) {PetscOptionsInt("-da_refine_y","Refinement ratio in y direction","DMDASetRefinementFactor",dd->refine_y,&dd->refine_y,NULL);}
53: if (dd->dim > 2) {PetscOptionsInt("-da_refine_z","Refinement ratio in z direction","DMDASetRefinementFactor",dd->refine_z,&dd->refine_z,NULL);}
54: dd->coarsen_x = dd->refine_x; dd->coarsen_y = dd->refine_y; dd->coarsen_z = dd->refine_z;
56: /* Get refinement factors, defaults taken from the coarse DMDA */
57: PetscMalloc3(maxnlevels,PetscInt,&refx,maxnlevels,PetscInt,&refy,maxnlevels,PetscInt,&refz);
58: DMDAGetRefinementFactor(da,&refx[0],&refy[0],&refz[0]);
59: for (i=1; i<maxnlevels; i++) {
60: refx[i] = refx[0];
61: refy[i] = refy[0];
62: refz[i] = refz[0];
63: }
64: n = maxnlevels;
65: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
66: if (da->levelup - da->leveldown >= 0) dd->refine_x = refx[da->levelup - da->leveldown];
67: if (da->levelup - da->leveldown >= 1) dd->coarsen_x = refx[da->levelup - da->leveldown - 1];
68: if (dd->dim > 1) {
69: n = maxnlevels;
70: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
71: if (da->levelup - da->leveldown >= 0) dd->refine_y = refy[da->levelup - da->leveldown];
72: if (da->levelup - da->leveldown >= 1) dd->coarsen_y = refy[da->levelup - da->leveldown - 1];
73: }
74: if (dd->dim > 2) {
75: n = maxnlevels;
76: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);
77: if (da->levelup - da->leveldown >= 0) dd->refine_z = refz[da->levelup - da->leveldown];
78: if (da->levelup - da->leveldown >= 1) dd->coarsen_z = refz[da->levelup - da->leveldown - 1];
79: }
81: if (negativeMNP) {PetscOptionsInt("-da_refine","Uniformly refine DA one or more times","None",refine,&refine,NULL);}
82: PetscOptionsTail();
84: while (refine--) {
85: if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
86: dd->M = dd->refine_x*dd->M;
87: } else {
88: dd->M = 1 + dd->refine_x*(dd->M - 1);
89: }
90: if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
91: dd->N = dd->refine_y*dd->N;
92: } else {
93: dd->N = 1 + dd->refine_y*(dd->N - 1);
94: }
95: if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
96: dd->P = dd->refine_z*dd->P;
97: } else {
98: dd->P = 1 + dd->refine_z*(dd->P - 1);
99: }
100: da->levelup++;
101: if (da->levelup - da->leveldown >= 0) {
102: dd->refine_x = refx[da->levelup - da->leveldown];
103: dd->refine_y = refy[da->levelup - da->leveldown];
104: dd->refine_z = refz[da->levelup - da->leveldown];
105: }
106: if (da->levelup - da->leveldown >= 1) {
107: dd->coarsen_x = refx[da->levelup - da->leveldown - 1];
108: dd->coarsen_y = refy[da->levelup - da->leveldown - 1];
109: dd->coarsen_z = refz[da->levelup - da->leveldown - 1];
110: }
111: }
112: PetscFree3(refx,refy,refz);
113: return(0);
114: }
116: extern PetscErrorCode DMCreateGlobalVector_DA(DM,Vec*);
117: extern PetscErrorCode DMCreateLocalVector_DA(DM,Vec*);
118: extern PetscErrorCode DMGlobalToLocalBegin_DA(DM,Vec,InsertMode,Vec);
119: extern PetscErrorCode DMGlobalToLocalEnd_DA(DM,Vec,InsertMode,Vec);
120: extern PetscErrorCode DMLocalToGlobalBegin_DA(DM,Vec,InsertMode,Vec);
121: extern PetscErrorCode DMLocalToGlobalEnd_DA(DM,Vec,InsertMode,Vec);
122: extern PetscErrorCode DMCreateInterpolation_DA(DM,DM,Mat*,Vec*);
123: extern PetscErrorCode DMCreateColoring_DA(DM,ISColoringType,MatType,ISColoring*);
124: extern PetscErrorCode DMCreateMatrix_DA(DM,MatType,Mat*);
125: extern PetscErrorCode DMCreateCoordinateDM_DA(DM,DM*);
126: extern PetscErrorCode DMRefine_DA(DM,MPI_Comm,DM*);
127: extern PetscErrorCode DMCoarsen_DA(DM,MPI_Comm,DM*);
128: extern PetscErrorCode DMRefineHierarchy_DA(DM,PetscInt,DM[]);
129: extern PetscErrorCode DMCoarsenHierarchy_DA(DM,PetscInt,DM[]);
130: extern PetscErrorCode DMCreateInjection_DA(DM,DM,VecScatter*);
131: extern PetscErrorCode DMCreateAggregates_DA(DM,DM,Mat*);
132: extern PetscErrorCode DMView_DA(DM,PetscViewer);
133: extern PetscErrorCode DMSetUp_DA(DM);
134: extern PetscErrorCode DMDestroy_DA(DM);
135: extern PetscErrorCode DMCreateDomainDecomposition_DA(DM,PetscInt*,char***,IS**,IS**,DM**);
136: extern PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM,PetscInt,DM*,VecScatter**,VecScatter**,VecScatter**);
140: PetscErrorCode DMLoad_DA(DM da,PetscViewer viewer)
141: {
142: PetscErrorCode ierr;
143: PetscInt dim,m,n,p,dof,swidth;
144: DMDAStencilType stencil;
145: DMDABoundaryType bx,by,bz;
146: PetscBool coors;
147: DM dac;
148: Vec c;
151: PetscViewerBinaryRead(viewer,&dim,1,PETSC_INT);
152: PetscViewerBinaryRead(viewer,&m,1,PETSC_INT);
153: PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);
154: PetscViewerBinaryRead(viewer,&p,1,PETSC_INT);
155: PetscViewerBinaryRead(viewer,&dof,1,PETSC_INT);
156: PetscViewerBinaryRead(viewer,&swidth,1,PETSC_INT);
157: PetscViewerBinaryRead(viewer,&bx,1,PETSC_ENUM);
158: PetscViewerBinaryRead(viewer,&by,1,PETSC_ENUM);
159: PetscViewerBinaryRead(viewer,&bz,1,PETSC_ENUM);
160: PetscViewerBinaryRead(viewer,&stencil,1,PETSC_ENUM);
162: DMDASetDim(da, dim);
163: DMDASetSizes(da, m,n,p);
164: DMDASetBoundaryType(da, bx, by, bz);
165: DMDASetDof(da, dof);
166: DMDASetStencilType(da, stencil);
167: DMDASetStencilWidth(da, swidth);
168: DMSetUp(da);
169: PetscViewerBinaryRead(viewer,&coors,1,PETSC_ENUM);
170: if (coors) {
171: DMGetCoordinateDM(da,&dac);
172: DMCreateGlobalVector(dac,&c);
173: VecLoad(c,viewer);
174: DMSetCoordinates(da,c);
175: VecDestroy(&c);
176: }
177: return(0);
178: }
182: PetscErrorCode DMCreateSubDM_DA(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
183: {
185: DM_DA *da = (DM_DA*)dm->data;
188: if (da->dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support only implemented for 2d");
189: if (subdm) {
190: DMDACreate2d(PetscObjectComm((PetscObject)dm),da->bx,da->by,da->stencil_type,da->M,da->N,da->m,da->n,numFields,da->s,da->lx,da->ly,subdm);
191: }
192: if (is) {
193: PetscInt *indices,cnt = 0, dof = da->w,i,j;
194: PetscMalloc(sizeof(PetscInt)*da->Nlocal*numFields/dof,&indices);
195: for (i=da->base/dof; i<(da->base+da->Nlocal)/dof; i++) {
196: for (j=0; j<numFields; j++) {
197: indices[cnt++] = dof*i + fields[j];
198: }
199: }
200: if (cnt != da->Nlocal*numFields/dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not equal expected value");
201: ISCreateGeneral(PetscObjectComm((PetscObject)dm),da->Nlocal*numFields/dof,indices,PETSC_OWN_POINTER,is);
202: }
203: return(0);
204: }
208: PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
209: {
210: PetscInt i;
212: DM_DA *dd = (DM_DA*)dm->data;
213: PetscInt dof = dd->w;
216: if (len) *len = dof;
217: if (islist) {
218: Vec v;
219: PetscInt rstart,n;
221: DMGetGlobalVector(dm,&v);
222: VecGetOwnershipRange(v,&rstart,NULL);
223: VecGetLocalSize(v,&n);
224: DMRestoreGlobalVector(dm,&v);
225: PetscMalloc(dof*sizeof(IS),islist);
226: for (i=0; i<dof; i++) {
227: ISCreateStride(PetscObjectComm((PetscObject)dm),n/dof,rstart+i,dof,&(*islist)[i]);
228: }
229: }
230: if (namelist) {
231: PetscMalloc(dof*sizeof(const char*), namelist);
232: if (dd->fieldname) {
233: for (i=0; i<dof; i++) {
234: PetscStrallocpy(dd->fieldname[i],&(*namelist)[i]);
235: }
236: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently DMDA must have fieldnames");
237: }
238: if (dmlist) {
239: DM da;
241: DMDACreate(PetscObjectComm((PetscObject)dm), &da);
242: DMDASetDim(da, dd->dim);
243: DMDASetSizes(da, dd->M, dd->N, dd->P);
244: DMDASetNumProcs(da, dd->m, dd->n, dd->p);
245: DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);
246: DMDASetDof(da, 1);
247: DMDASetStencilType(da, dd->stencil_type);
248: DMDASetStencilWidth(da, dd->s);
249: DMSetUp(da);
250: PetscMalloc(dof*sizeof(DM),dmlist);
251: for (i=0; i<dof-1; i++) {PetscObjectReference((PetscObject)da);}
252: for (i=0; i<dof; i++) (*dmlist)[i] = da;
253: }
254: return(0);
255: }
257: /*MC
258: DMDA = "da" - A DM object that is used to manage data for a structured grid in 1, 2, or 3 dimensions.
259: In the global representation of the vector each process stores a non-overlapping rectangular (or slab in 3d) portion of the grid points.
260: In the local representation these rectangular regions (slabs) are extended in all directions by a stencil width.
262: 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
263: vertex centered.
266: Level: intermediate
268: .seealso: DMType, DMCOMPOSITE, DMDACreate(), DMCreate(), DMSetType()
269: M*/
274: PETSC_EXTERN PetscErrorCode DMCreate_DA(DM da)
275: {
277: DM_DA *dd;
281: PetscNewLog(da,DM_DA,&dd);
282: da->data = dd;
284: dd->dim = -1;
285: dd->interptype = DMDA_Q1;
286: dd->refine_x = 2;
287: dd->refine_y = 2;
288: dd->refine_z = 2;
289: dd->coarsen_x = 2;
290: dd->coarsen_y = 2;
291: dd->coarsen_z = 2;
292: dd->fieldname = NULL;
293: dd->nlocal = -1;
294: dd->Nlocal = -1;
295: dd->M = -1;
296: dd->N = -1;
297: dd->P = -1;
298: dd->m = -1;
299: dd->n = -1;
300: dd->p = -1;
301: dd->w = -1;
302: dd->s = -1;
304: dd->xs = -1; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1;
305: dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1;
307: dd->Nsub = 1;
308: dd->xol = 0;
309: dd->yol = 0;
310: dd->zol = 0;
311: dd->xo = 0;
312: dd->yo = 0;
313: dd->zo = 0;
314: dd->Mo = -1;
315: dd->No = -1;
316: dd->Po = -1;
318: dd->gtol = NULL;
319: dd->ltog = NULL;
320: dd->ltol = NULL;
321: dd->ao = NULL;
322: dd->base = -1;
323: dd->bx = DMDA_BOUNDARY_NONE;
324: dd->by = DMDA_BOUNDARY_NONE;
325: dd->bz = DMDA_BOUNDARY_NONE;
326: dd->stencil_type = DMDA_STENCIL_BOX;
327: dd->interptype = DMDA_Q1;
328: dd->idx = NULL;
329: dd->Nl = -1;
330: dd->lx = NULL;
331: dd->ly = NULL;
332: dd->lz = NULL;
334: dd->elementtype = DMDA_ELEMENT_Q1;
336: PetscStrallocpy(VECSTANDARD,(char**)&da->vectype);
338: da->ops->globaltolocalbegin = DMGlobalToLocalBegin_DA;
339: da->ops->globaltolocalend = DMGlobalToLocalEnd_DA;
340: da->ops->localtoglobalbegin = DMLocalToGlobalBegin_DA;
341: da->ops->localtoglobalend = DMLocalToGlobalEnd_DA;
342: da->ops->createglobalvector = DMCreateGlobalVector_DA;
343: da->ops->createlocalvector = DMCreateLocalVector_DA;
344: da->ops->createinterpolation = DMCreateInterpolation_DA;
345: da->ops->getcoloring = DMCreateColoring_DA;
346: da->ops->creatematrix = DMCreateMatrix_DA;
347: da->ops->refine = DMRefine_DA;
348: da->ops->coarsen = DMCoarsen_DA;
349: da->ops->refinehierarchy = DMRefineHierarchy_DA;
350: da->ops->coarsenhierarchy = DMCoarsenHierarchy_DA;
351: da->ops->getinjection = DMCreateInjection_DA;
352: da->ops->getaggregates = DMCreateAggregates_DA;
353: da->ops->destroy = DMDestroy_DA;
354: da->ops->view = 0;
355: da->ops->setfromoptions = DMSetFromOptions_DA;
356: da->ops->setup = DMSetUp_DA;
357: da->ops->load = DMLoad_DA;
358: da->ops->createcoordinatedm = DMCreateCoordinateDM_DA;
359: da->ops->createsubdm = DMCreateSubDM_DA;
360: da->ops->createfielddecomposition = DMCreateFieldDecomposition_DA;
361: da->ops->createdomaindecomposition = DMCreateDomainDecomposition_DA;
362: da->ops->createddscatters = DMCreateDomainDecompositionScatters_DA;
363: return(0);
364: }
368: /*@
369: DMDACreate - Creates a DMDA object.
371: Collective on MPI_Comm
373: Input Parameter:
374: . comm - The communicator for the DMDA object
376: Output Parameter:
377: . da - The DMDA object
379: Level: advanced
381: Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist?
383: .keywords: DMDA, create
384: .seealso: DMDASetSizes(), DMDADuplicate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
385: @*/
386: PetscErrorCode DMDACreate(MPI_Comm comm, DM *da)
387: {
392: DMCreate(comm,da);
393: DMSetType(*da,DMDA);
394: return(0);
395: }