Actual source code: dacreate.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/daimpl.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;
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,PETSC_NULL);}
34: if (bN) {PetscOptionsInt("-da_grid_y","Number of grid points in y direction","DMDASetSizes",dd->N,&dd->N,PETSC_NULL);}
35: if (bP) {PetscOptionsInt("-da_grid_z","Number of grid points in z direction","DMDASetSizes",dd->P,&dd->P,PETSC_NULL);}
36: /* Handle DMDA parallel distibution */
37: PetscOptionsInt("-da_processors_x","Number of processors in x direction","DMDASetNumProcs",dd->m,&dd->m,PETSC_NULL);
38: if (dd->dim > 1) {PetscOptionsInt("-da_processors_y","Number of processors in y direction","DMDASetNumProcs",dd->n,&dd->n,PETSC_NULL);}
39: if (dd->dim > 2) {PetscOptionsInt("-da_processors_z","Number of processors in z direction","DMDASetNumProcs",dd->p,&dd->p,PETSC_NULL);}
40: /* Handle DMDA refinement */
41: PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DMDASetRefinementFactor",dd->refine_x,&dd->refine_x,PETSC_NULL);
42: if (dd->dim > 1) {PetscOptionsInt("-da_refine_y","Refinement ratio in y direction","DMDASetRefinementFactor",dd->refine_y,&dd->refine_y,PETSC_NULL);}
43: if (dd->dim > 2) {PetscOptionsInt("-da_refine_z","Refinement ratio in z direction","DMDASetRefinementFactor",dd->refine_z,&dd->refine_z,PETSC_NULL);}
44: dd->coarsen_x = dd->refine_x; dd->coarsen_y = dd->refine_y; dd->coarsen_z = dd->refine_z;
46: /* Get refinement factors, defaults taken from the coarse DMDA */
47: PetscMalloc3(maxnlevels,PetscInt,&refx,maxnlevels,PetscInt,&refy,maxnlevels,PetscInt,&refz);
48: DMDAGetRefinementFactor(da,&refx[0],&refy[0],&refz[0]);
49: for (i=1; i<maxnlevels; i++) {
50: refx[i] = refx[0];
51: refy[i] = refy[0];
52: refz[i] = refz[0];
53: }
54: n = maxnlevels;
55: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,PETSC_NULL);
56: if (da->levelup - da->leveldown >= 0) dd->refine_x = refx[da->levelup - da->leveldown];
57: if (da->levelup - da->leveldown >= 1) dd->coarsen_x = refx[da->levelup - da->leveldown - 1];
58: if (dd->dim > 1) {
59: n = maxnlevels;
60: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,PETSC_NULL);
61: if (da->levelup - da->leveldown >= 0) dd->refine_y = refy[da->levelup - da->leveldown];
62: if (da->levelup - da->leveldown >= 1) dd->coarsen_y = refy[da->levelup - da->leveldown - 1];
63: }
64: if (dd->dim > 2) {
65: n = maxnlevels;
66: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,PETSC_NULL);
67: if (da->levelup - da->leveldown >= 0) dd->refine_z = refz[da->levelup - da->leveldown];
68: if (da->levelup - da->leveldown >= 1) dd->coarsen_z = refz[da->levelup - da->leveldown - 1];
69: }
71: if (negativeMNP) {PetscOptionsInt("-da_refine","Uniformly refine DA one or more times","None",refine,&refine,PETSC_NULL);}
72: PetscOptionsTail();
74: while (refine--) {
75: if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){
76: dd->M = dd->refine_x*dd->M;
77: } else {
78: dd->M = 1 + dd->refine_x*(dd->M - 1);
79: }
80: if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){
81: dd->N = dd->refine_y*dd->N;
82: } else {
83: dd->N = 1 + dd->refine_y*(dd->N - 1);
84: }
85: if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){
86: dd->P = dd->refine_z*dd->P;
87: } else {
88: dd->P = 1 + dd->refine_z*(dd->P - 1);
89: }
90: da->levelup++;
91: if (da->levelup - da->leveldown >= 0) {
92: dd->refine_x = refx[da->levelup - da->leveldown];
93: dd->refine_y = refy[da->levelup - da->leveldown];
94: dd->refine_z = refz[da->levelup - da->leveldown];
95: }
96: if (da->levelup - da->leveldown >= 1) {
97: dd->coarsen_x = refx[da->levelup - da->leveldown - 1];
98: dd->coarsen_y = refy[da->levelup - da->leveldown - 1];
99: dd->coarsen_z = refz[da->levelup - da->leveldown - 1];
100: }
101: }
102: PetscFree3(refx,refy,refz);
103: return(0);
104: }
106: extern PetscErrorCode DMCreateGlobalVector_DA(DM,Vec*);
107: extern PetscErrorCode DMCreateLocalVector_DA(DM,Vec*);
108: extern PetscErrorCode DMGlobalToLocalBegin_DA(DM,Vec,InsertMode,Vec);
109: extern PetscErrorCode DMGlobalToLocalEnd_DA(DM,Vec,InsertMode,Vec);
110: extern PetscErrorCode DMLocalToGlobalBegin_DA(DM,Vec,InsertMode,Vec);
111: extern PetscErrorCode DMLocalToGlobalEnd_DA(DM,Vec,InsertMode,Vec);
112: extern PetscErrorCode DMCreateInterpolation_DA(DM,DM,Mat*,Vec*);
113: extern PetscErrorCode DMCreateColoring_DA(DM,ISColoringType,const MatType,ISColoring*);
114: extern PetscErrorCode DMCreateMatrix_DA(DM,const MatType,Mat*);
115: extern PetscErrorCode DMRefine_DA(DM,MPI_Comm,DM*);
116: extern PetscErrorCode DMCoarsen_DA(DM,MPI_Comm,DM*);
117: extern PetscErrorCode DMRefineHierarchy_DA(DM,PetscInt,DM[]);
118: extern PetscErrorCode DMCoarsenHierarchy_DA(DM,PetscInt,DM[]);
119: extern PetscErrorCode DMCreateInjection_DA(DM,DM,VecScatter*);
120: extern PetscErrorCode DMCreateAggregates_DA(DM,DM,Mat*);
121: extern PetscErrorCode DMView_DA(DM,PetscViewer);
122: extern PetscErrorCode DMSetUp_DA(DM);
123: extern PetscErrorCode DMDestroy_DA(DM);
127: PetscErrorCode DMLoad_DA(DM da,PetscViewer viewer)
128: {
129: PetscErrorCode ierr;
130: PetscInt dim,m,n,p,dof,swidth;
131: DMDAStencilType stencil;
132: DMDABoundaryType bx,by,bz;
133: PetscInt classid = DM_FILE_CLASSID,subclassid = DMDA_FILE_CLASSID;
134: PetscBool coors;
135: DM dac;
136: Vec c;
139: PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);
140: if (classid != DM_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not DM next in file");
141: PetscViewerBinaryRead(viewer,&subclassid,1,PETSC_INT);
142: if (subclassid != DMDA_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not DM DA next in file");
143: PetscViewerBinaryRead(viewer,&dim,1,PETSC_INT);
144: PetscViewerBinaryRead(viewer,&m,1,PETSC_INT);
145: PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);
146: PetscViewerBinaryRead(viewer,&p,1,PETSC_INT);
147: PetscViewerBinaryRead(viewer,&dof,1,PETSC_INT);
148: PetscViewerBinaryRead(viewer,&swidth,1,PETSC_INT);
149: PetscViewerBinaryRead(viewer,&bx,1,PETSC_ENUM);
150: PetscViewerBinaryRead(viewer,&by,1,PETSC_ENUM);
151: PetscViewerBinaryRead(viewer,&bz,1,PETSC_ENUM);
152: PetscViewerBinaryRead(viewer,&stencil,1,PETSC_ENUM);
154: DMDASetDim(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,PETSC_ENUM);
162: if (coors) {
163: DMDAGetCoordinateDA(da,&dac);
164: DMCreateGlobalVector(dac,&c);
165: VecLoad(c,viewer);
166: DMDASetCoordinates(da,c);
167: VecDestroy(&c);
168: }
169: return(0);
170: }
174: PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len,char ***namelist, IS **islist, DM** dmlist)
175: {
176: PetscInt i;
178: DM_DA *dd = (DM_DA*)dm->data;
179: PetscInt dof = dd->w;
182: if(len) *len = dof;
183: if (islist) {
184: Vec v;
185: PetscInt rstart,n;
187: DMGetGlobalVector(dm,&v);
188: VecGetOwnershipRange(v,&rstart,PETSC_NULL);
189: VecGetLocalSize(v,&n);
190: DMRestoreGlobalVector(dm,&v);
191: PetscMalloc(dof*sizeof(IS),islist);
192: for (i=0; i<dof; i++) {
193: ISCreateStride(((PetscObject)dm)->comm,n/dof,rstart+i,dof,&(*islist)[i]);
194: }
195: }
196: if (namelist) {
197: PetscMalloc(dof*sizeof(const char *), namelist);
198: if (dd->fieldname) {
199: for (i=0; i<dof; i++) {
200: PetscStrallocpy(dd->fieldname[i],&(*namelist)[i]);
201: }
202: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently DMDA must have fieldnames");
203: }
204: if (dmlist) {
205: DM da;
207: DMDACreate(((PetscObject)dm)->comm, &da);
208: DMDASetDim(da, dd->dim);
209: DMDASetSizes(da, dd->M, dd->N, dd->P);
210: DMDASetNumProcs(da, dd->m, dd->n, dd->p);
211: DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);
212: DMDASetDof(da, 1);
213: DMDASetStencilType(da, dd->stencil_type);
214: DMDASetStencilWidth(da, dd->s);
215: DMSetUp(da);
216: PetscMalloc(dof*sizeof(DM),dmlist);
217: for (i=0; i<dof-1; i++) {PetscObjectReference((PetscObject)da);}
218: for (i=0; i<dof; i++) (*dmlist)[i] = da;
219: }
221: return(0);
222: }
225: EXTERN_C_BEGIN
228: PetscErrorCode DMCreate_DA(DM da)
229: {
231: DM_DA *dd;
235: PetscNewLog(da,DM_DA,&dd);
236: da->data = dd;
238: dd->dim = -1;
239: dd->interptype = DMDA_Q1;
240: dd->refine_x = 2;
241: dd->refine_y = 2;
242: dd->refine_z = 2;
243: dd->coarsen_x = 2;
244: dd->coarsen_y = 2;
245: dd->coarsen_z = 2;
246: dd->fieldname = PETSC_NULL;
247: dd->nlocal = -1;
248: dd->Nlocal = -1;
249: dd->M = -1;
250: dd->N = -1;
251: dd->P = -1;
252: dd->m = -1;
253: dd->n = -1;
254: dd->p = -1;
255: dd->w = -1;
256: dd->s = -1;
257: dd->xs = -1; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1;
258: dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1;
260: dd->gtol = PETSC_NULL;
261: dd->ltog = PETSC_NULL;
262: dd->ltol = PETSC_NULL;
263: dd->ao = PETSC_NULL;
264: dd->base = -1;
265: dd->bx = DMDA_BOUNDARY_NONE;
266: dd->by = DMDA_BOUNDARY_NONE;
267: dd->bz = DMDA_BOUNDARY_NONE;
268: dd->stencil_type = DMDA_STENCIL_BOX;
269: dd->interptype = DMDA_Q1;
270: dd->idx = PETSC_NULL;
271: dd->Nl = -1;
272: dd->lx = PETSC_NULL;
273: dd->ly = PETSC_NULL;
274: dd->lz = PETSC_NULL;
276: dd->elementtype = DMDA_ELEMENT_Q1;
278: PetscStrallocpy(VECSTANDARD,&da->vectype);
279: da->ops->globaltolocalbegin = DMGlobalToLocalBegin_DA;
280: da->ops->globaltolocalend = DMGlobalToLocalEnd_DA;
281: da->ops->localtoglobalbegin = DMLocalToGlobalBegin_DA;
282: da->ops->localtoglobalend = DMLocalToGlobalEnd_DA;
283: da->ops->createglobalvector = DMCreateGlobalVector_DA;
284: da->ops->createlocalvector = DMCreateLocalVector_DA;
285: da->ops->createinterpolation = DMCreateInterpolation_DA;
286: da->ops->getcoloring = DMCreateColoring_DA;
287: da->ops->creatematrix = DMCreateMatrix_DA;
288: da->ops->refine = DMRefine_DA;
289: da->ops->coarsen = DMCoarsen_DA;
290: da->ops->refinehierarchy = DMRefineHierarchy_DA;
291: da->ops->coarsenhierarchy = DMCoarsenHierarchy_DA;
292: da->ops->getinjection = DMCreateInjection_DA;
293: da->ops->getaggregates = DMCreateAggregates_DA;
294: da->ops->destroy = DMDestroy_DA;
295: da->ops->view = 0;
296: da->ops->setfromoptions = DMSetFromOptions_DA;
297: da->ops->setup = DMSetUp_DA;
298: da->ops->load = DMLoad_DA;
299: da->ops->createfielddecomposition = DMCreateFieldDecomposition_DA;
300: return(0);
301: }
302: EXTERN_C_END
306: /*@
307: DMDACreate - Creates a DMDA object.
309: Collective on MPI_Comm
311: Input Parameter:
312: . comm - The communicator for the DMDA object
314: Output Parameter:
315: . da - The DMDA object
317: Level: advanced
319: Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist?
321: .keywords: DMDA, create
322: .seealso: DMDASetSizes(), DMDADuplicate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
323: @*/
324: PetscErrorCode DMDACreate(MPI_Comm comm, DM *da)
325: {
330: DMCreate(comm,da);
331: DMSetType(*da,DMDA);
332: return(0);
333: }