xref: /petsc/src/dm/impls/da/da.c (revision 8865f1ea4ea560cd84ab8db62e98b7095cdff96f)
1b45d2f2cSJed Brown #include <petsc-private/daimpl.h>    /*I   "petscdmda.h"   I*/
247c6ae99SBarry Smith 
347c6ae99SBarry Smith 
447c6ae99SBarry Smith #undef __FUNCT__
5aa219208SBarry Smith #define __FUNCT__ "DMDASetDim"
647c6ae99SBarry Smith /*@
7aa219208SBarry Smith   DMDASetDim - Sets the dimension
847c6ae99SBarry Smith 
9aa219208SBarry Smith   Collective on DMDA
1047c6ae99SBarry Smith 
1147c6ae99SBarry Smith   Input Parameters:
12aa219208SBarry Smith + da - the DMDA
1347c6ae99SBarry Smith - dim - the dimension (or PETSC_DECIDE)
1447c6ae99SBarry Smith 
1547c6ae99SBarry Smith   Level: intermediate
1647c6ae99SBarry Smith 
17aa219208SBarry Smith .seealso: DaGetDim(), DMDASetSizes()
1847c6ae99SBarry Smith @*/
197087cfbeSBarry Smith PetscErrorCode  DMDASetDim(DM da, PetscInt dim)
2047c6ae99SBarry Smith {
2147c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
2247c6ae99SBarry Smith 
2347c6ae99SBarry Smith   PetscFunctionBegin;
2447c6ae99SBarry Smith   PetscValidHeaderSpecific(da, DM_CLASSID, 1);
25aa219208SBarry Smith   if (dd->dim > 0 && dim != dd->dim) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change DMDA dim from %D after it was set to %D",dd->dim,dim);
2647c6ae99SBarry Smith   dd->dim = dim;
2747c6ae99SBarry Smith   PetscFunctionReturn(0);
2847c6ae99SBarry Smith }
2947c6ae99SBarry Smith 
3047c6ae99SBarry Smith #undef __FUNCT__
31aa219208SBarry Smith #define __FUNCT__ "DMDASetSizes"
3247c6ae99SBarry Smith /*@
33aa219208SBarry Smith   DMDASetSizes - Sets the global sizes
3447c6ae99SBarry Smith 
35aa219208SBarry Smith   Logically Collective on DMDA
3647c6ae99SBarry Smith 
3747c6ae99SBarry Smith   Input Parameters:
38aa219208SBarry Smith + da - the DMDA
3947c6ae99SBarry Smith . M - the global X size (or PETSC_DECIDE)
4047c6ae99SBarry Smith . N - the global Y size (or PETSC_DECIDE)
4147c6ae99SBarry Smith - P - the global Z size (or PETSC_DECIDE)
4247c6ae99SBarry Smith 
4347c6ae99SBarry Smith   Level: intermediate
4447c6ae99SBarry Smith 
45aa219208SBarry Smith .seealso: DMDAGetSize(), PetscSplitOwnership()
4647c6ae99SBarry Smith @*/
477087cfbeSBarry Smith PetscErrorCode  DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
4847c6ae99SBarry Smith {
4947c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
5047c6ae99SBarry Smith 
5147c6ae99SBarry Smith   PetscFunctionBegin;
5247c6ae99SBarry Smith   PetscValidHeaderSpecific(da, DM_CLASSID, 1);
5347c6ae99SBarry Smith   PetscValidLogicalCollectiveInt(da,M,2);
5447c6ae99SBarry Smith   PetscValidLogicalCollectiveInt(da,N,3);
5547c6ae99SBarry Smith   PetscValidLogicalCollectiveInt(da,P,4);
56cb630486SJed Brown   if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
5747c6ae99SBarry Smith 
5847c6ae99SBarry Smith   dd->M = M;
5947c6ae99SBarry Smith   dd->N = N;
6047c6ae99SBarry Smith   dd->P = P;
6147c6ae99SBarry Smith   PetscFunctionReturn(0);
6247c6ae99SBarry Smith }
6347c6ae99SBarry Smith 
6447c6ae99SBarry Smith #undef __FUNCT__
65aa219208SBarry Smith #define __FUNCT__ "DMDASetNumProcs"
6647c6ae99SBarry Smith /*@
67aa219208SBarry Smith   DMDASetNumProcs - Sets the number of processes in each dimension
6847c6ae99SBarry Smith 
69aa219208SBarry Smith   Logically Collective on DMDA
7047c6ae99SBarry Smith 
7147c6ae99SBarry Smith   Input Parameters:
72aa219208SBarry Smith + da - the DMDA
7347c6ae99SBarry Smith . m - the number of X procs (or PETSC_DECIDE)
7447c6ae99SBarry Smith . n - the number of Y procs (or PETSC_DECIDE)
7547c6ae99SBarry Smith - p - the number of Z procs (or PETSC_DECIDE)
7647c6ae99SBarry Smith 
7747c6ae99SBarry Smith   Level: intermediate
7847c6ae99SBarry Smith 
79aa219208SBarry Smith .seealso: DMDASetSizes(), DMDAGetSize(), PetscSplitOwnership()
8047c6ae99SBarry Smith @*/
817087cfbeSBarry Smith PetscErrorCode  DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
8247c6ae99SBarry Smith {
8347c6ae99SBarry Smith   DM_DA          *dd = (DM_DA*)da->data;
84e3f3e4b6SBarry Smith   PetscErrorCode ierr;
8547c6ae99SBarry Smith 
8647c6ae99SBarry Smith   PetscFunctionBegin;
8747c6ae99SBarry Smith   PetscValidHeaderSpecific(da, DM_CLASSID, 1);
8847c6ae99SBarry Smith   PetscValidLogicalCollectiveInt(da,m,2);
8947c6ae99SBarry Smith   PetscValidLogicalCollectiveInt(da,n,3);
9047c6ae99SBarry Smith   PetscValidLogicalCollectiveInt(da,p,4);
91cb630486SJed Brown   if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
9247c6ae99SBarry Smith   dd->m = m;
9347c6ae99SBarry Smith   dd->n = n;
9447c6ae99SBarry Smith   dd->p = p;
95e3f3e4b6SBarry Smith   if (dd->dim == 2) {
96e3f3e4b6SBarry Smith     PetscMPIInt size;
97e3f3e4b6SBarry Smith     ierr = MPI_Comm_size(((PetscObject)da)->comm,&size);CHKERRQ(ierr);
98e3f3e4b6SBarry Smith     if ((dd->m > 0) && (dd->n < 0)) {
99e3f3e4b6SBarry Smith       dd->n = size/dd->m;
100e3f3e4b6SBarry Smith       if (dd->n*dd->m != size) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"%D processes in X direction not divisible into comm size %d",m,size);
101e3f3e4b6SBarry Smith     }
102e3f3e4b6SBarry Smith     if ((dd->n > 0) && (dd->m < 0)) {
103e3f3e4b6SBarry Smith       dd->m = size/dd->n;
104e3f3e4b6SBarry Smith       if (dd->n*dd->m != size) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"%D processes in Y direction not divisible into comm size %d",n,size);
105e3f3e4b6SBarry Smith     }
106e3f3e4b6SBarry Smith   }
10747c6ae99SBarry Smith   PetscFunctionReturn(0);
10847c6ae99SBarry Smith }
10947c6ae99SBarry Smith 
11047c6ae99SBarry Smith #undef __FUNCT__
1113e24977eSEthan Coon #define __FUNCT__ "DMDASetBoundaryType"
11247c6ae99SBarry Smith /*@
1131321219cSEthan Coon   DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries.
11447c6ae99SBarry Smith 
11547c6ae99SBarry Smith   Not collective
11647c6ae99SBarry Smith 
11747c6ae99SBarry Smith   Input Parameter:
118aa219208SBarry Smith + da    - The DMDA
1191321219cSEthan Coon - bx,by,bz - One of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC
12047c6ae99SBarry Smith 
12147c6ae99SBarry Smith   Level: intermediate
12247c6ae99SBarry Smith 
12347c6ae99SBarry Smith .keywords:  distributed array, periodicity
12496e147daSBarry Smith .seealso: DMDACreate(), DMDestroy(), DMDA, DMDABoundaryType
12547c6ae99SBarry Smith @*/
1261321219cSEthan Coon PetscErrorCode  DMDASetBoundaryType(DM da,DMDABoundaryType bx,DMDABoundaryType by,DMDABoundaryType bz)
12747c6ae99SBarry Smith {
12847c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
12947c6ae99SBarry Smith 
13047c6ae99SBarry Smith   PetscFunctionBegin;
13147c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
1325a43f728SLisandro Dalcin   PetscValidLogicalCollectiveEnum(da,bx,2);
1335a43f728SLisandro Dalcin   PetscValidLogicalCollectiveEnum(da,by,3);
1345a43f728SLisandro Dalcin   PetscValidLogicalCollectiveEnum(da,bz,4);
135cb630486SJed Brown   if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1361321219cSEthan Coon   dd->bx = bx;
1371321219cSEthan Coon   dd->by = by;
1381321219cSEthan Coon   dd->bz = bz;
13947c6ae99SBarry Smith   PetscFunctionReturn(0);
14047c6ae99SBarry Smith }
14147c6ae99SBarry Smith 
14247c6ae99SBarry Smith #undef __FUNCT__
143aa219208SBarry Smith #define __FUNCT__ "DMDASetDof"
14447c6ae99SBarry Smith /*@
145aa219208SBarry Smith   DMDASetDof - Sets the number of degrees of freedom per vertex
14647c6ae99SBarry Smith 
14747c6ae99SBarry Smith   Not collective
14847c6ae99SBarry Smith 
14947c6ae99SBarry Smith   Input Parameter:
150aa219208SBarry Smith + da  - The DMDA
15147c6ae99SBarry Smith - dof - Number of degrees of freedom
15247c6ae99SBarry Smith 
15347c6ae99SBarry Smith   Level: intermediate
15447c6ae99SBarry Smith 
15547c6ae99SBarry Smith .keywords:  distributed array, degrees of freedom
15696e147daSBarry Smith .seealso: DMDACreate(), DMDestroy(), DMDA
15747c6ae99SBarry Smith @*/
15854cfb0beSLisandro Dalcin PetscErrorCode  DMDASetDof(DM da, PetscInt dof)
15947c6ae99SBarry Smith {
16047c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
16147c6ae99SBarry Smith 
16247c6ae99SBarry Smith   PetscFunctionBegin;
16347c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
16454cfb0beSLisandro Dalcin   PetscValidLogicalCollectiveInt(da,dof,2);
165cb630486SJed Brown   if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
16647c6ae99SBarry Smith   dd->w  = dof;
1671411c6eeSJed Brown   da->bs = dof;
16847c6ae99SBarry Smith   PetscFunctionReturn(0);
16947c6ae99SBarry Smith }
17047c6ae99SBarry Smith 
17147c6ae99SBarry Smith #undef __FUNCT__
17288661749SPeter Brune #define __FUNCT__ "DMDASetOverlap"
17388661749SPeter Brune /*@
17488661749SPeter Brune   DMDASetOverlap - Sets the size of the per-processor overlap.
17588661749SPeter Brune 
17688661749SPeter Brune   Not collective
17788661749SPeter Brune 
17888661749SPeter Brune   Input Parameter:
17988661749SPeter Brune + da  - The DMDA
18088661749SPeter Brune - dof - Number of degrees of freedom
18188661749SPeter Brune 
18288661749SPeter Brune   Level: intermediate
18388661749SPeter Brune 
18488661749SPeter Brune .keywords:  distributed array, degrees of freedom
18588661749SPeter Brune .seealso: DMDACreate(), DMDestroy(), DMDA
18688661749SPeter Brune @*/
18788661749SPeter Brune PetscErrorCode  DMDASetOverlap(DM da, PetscInt overlap)
18888661749SPeter Brune {
18988661749SPeter Brune   DM_DA *dd = (DM_DA*)da->data;
19088661749SPeter Brune 
19188661749SPeter Brune   PetscFunctionBegin;
19288661749SPeter Brune   PetscValidHeaderSpecific(da,DM_CLASSID,1);
19388661749SPeter Brune   PetscValidLogicalCollectiveInt(da,overlap,2);
19488661749SPeter Brune   if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
19588661749SPeter Brune   dd->overlap = overlap;
19688661749SPeter Brune   PetscFunctionReturn(0);
19788661749SPeter Brune }
19888661749SPeter Brune 
199d886c4f4SPeter Brune #undef __FUNCT__
200d886c4f4SPeter Brune #define __FUNCT__ "DMDASetOffset"
201d886c4f4SPeter Brune /*@
202d886c4f4SPeter Brune   DMDASetOffset - Sets the index offset of the DA.
203d886c4f4SPeter Brune 
204d886c4f4SPeter Brune   Collective on DA
205d886c4f4SPeter Brune 
206d886c4f4SPeter Brune   Input Parameter:
207d886c4f4SPeter Brune + da  - The DMDA
208d886c4f4SPeter Brune . xo  - The offset in the x direction
209d886c4f4SPeter Brune . yo  - The offset in the y direction
210d886c4f4SPeter Brune - zo  - The offset in the z direction
211d886c4f4SPeter Brune 
212d886c4f4SPeter Brune   Level: intermediate
213d886c4f4SPeter Brune 
214d886c4f4SPeter Brune   Notes: This is used primarily to overlap a computation on a local DA with that on a global DA without
215d886c4f4SPeter Brune   changing boundary conditions or subdomain features that depend upon the global offsets.
216d886c4f4SPeter Brune 
217d886c4f4SPeter Brune .keywords:  distributed array, degrees of freedom
218d886c4f4SPeter Brune .seealso: DMDAGetOffset(), DMDAVecGetArray()
219d886c4f4SPeter Brune @*/
22095c13181SPeter Brune PetscErrorCode  DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
221d886c4f4SPeter Brune {
22295c13181SPeter Brune   PetscErrorCode ierr;
223d886c4f4SPeter Brune   DM_DA          *dd = (DM_DA*)da->data;
224d886c4f4SPeter Brune 
225d886c4f4SPeter Brune   PetscFunctionBegin;
226d886c4f4SPeter Brune   PetscValidHeaderSpecific(da,DM_CLASSID,1);
227d886c4f4SPeter Brune   PetscValidLogicalCollectiveInt(da,xo,2);
22895c13181SPeter Brune   PetscValidLogicalCollectiveInt(da,yo,3);
22995c13181SPeter Brune   PetscValidLogicalCollectiveInt(da,zo,4);
23095c13181SPeter Brune   PetscValidLogicalCollectiveInt(da,Mo,5);
23195c13181SPeter Brune   PetscValidLogicalCollectiveInt(da,No,6);
23295c13181SPeter Brune   PetscValidLogicalCollectiveInt(da,Po,7);
233d886c4f4SPeter Brune   dd->xo = xo;
234d886c4f4SPeter Brune   dd->yo = yo;
235d886c4f4SPeter Brune   dd->zo = zo;
23695c13181SPeter Brune   dd->Mo = Mo;
23795c13181SPeter Brune   dd->No = No;
23895c13181SPeter Brune   dd->Po = Po;
23995c13181SPeter Brune 
24095c13181SPeter Brune   if (da->coordinateDM) {
24195c13181SPeter Brune     ierr = DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);CHKERRQ(ierr);
24295c13181SPeter Brune   }
243d886c4f4SPeter Brune   PetscFunctionReturn(0);
244d886c4f4SPeter Brune }
245d886c4f4SPeter Brune 
246d886c4f4SPeter Brune #undef __FUNCT__
247d886c4f4SPeter Brune #define __FUNCT__ "DMDAGetOffset"
248d886c4f4SPeter Brune /*@
249d886c4f4SPeter Brune   DMDAGetOffset - Gets the index offset of the DA.
250d886c4f4SPeter Brune 
251d886c4f4SPeter Brune   Not collective
252d886c4f4SPeter Brune 
253d886c4f4SPeter Brune   Input Parameter:
254d886c4f4SPeter Brune . da  - The DMDA
255d886c4f4SPeter Brune 
256d886c4f4SPeter Brune   Output Parameters:
257d886c4f4SPeter Brune + xo  - The offset in the x direction
258d886c4f4SPeter Brune . yo  - The offset in the y direction
25995c13181SPeter Brune . zo  - The offset in the z direction
26095c13181SPeter Brune . Mo  - The global size in the x direction
26195c13181SPeter Brune . No  - The global size in the y direction
26295c13181SPeter Brune - Po  - The global size in the z direction
263d886c4f4SPeter Brune 
264d886c4f4SPeter Brune   Level: intermediate
265d886c4f4SPeter Brune 
266d886c4f4SPeter Brune .keywords:  distributed array, degrees of freedom
267d886c4f4SPeter Brune .seealso: DMDASetOffset(), DMDAVecGetArray()
268d886c4f4SPeter Brune @*/
26995c13181SPeter Brune PetscErrorCode  DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
270d886c4f4SPeter Brune {
271d886c4f4SPeter Brune   DM_DA *dd = (DM_DA*)da->data;
272d886c4f4SPeter Brune 
273d886c4f4SPeter Brune   PetscFunctionBegin;
274d886c4f4SPeter Brune   PetscValidHeaderSpecific(da,DM_CLASSID,1);
275d886c4f4SPeter Brune   if (xo) *xo = dd->xo;
276d886c4f4SPeter Brune   if (yo) *yo = dd->yo;
277d886c4f4SPeter Brune   if (zo) *zo = dd->zo;
27895c13181SPeter Brune   if (Mo) *Mo = dd->Mo;
27995c13181SPeter Brune   if (No) *No = dd->No;
28095c13181SPeter Brune   if (Po) *Po = dd->Po;
281d886c4f4SPeter Brune   PetscFunctionReturn(0);
282d886c4f4SPeter Brune }
283d886c4f4SPeter Brune 
28488661749SPeter Brune 
28588661749SPeter Brune #undef __FUNCT__
286aa219208SBarry Smith #define __FUNCT__ "DMDASetStencilType"
28747c6ae99SBarry Smith /*@
288aa219208SBarry Smith   DMDASetStencilType - Sets the type of the communication stencil
28947c6ae99SBarry Smith 
290aa219208SBarry Smith   Logically Collective on DMDA
29147c6ae99SBarry Smith 
29247c6ae99SBarry Smith   Input Parameter:
293aa219208SBarry Smith + da    - The DMDA
294aa219208SBarry Smith - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
29547c6ae99SBarry Smith 
29647c6ae99SBarry Smith   Level: intermediate
29747c6ae99SBarry Smith 
29847c6ae99SBarry Smith .keywords:  distributed array, stencil
29996e147daSBarry Smith .seealso: DMDACreate(), DMDestroy(), DMDA
30047c6ae99SBarry Smith @*/
3017087cfbeSBarry Smith PetscErrorCode  DMDASetStencilType(DM da, DMDAStencilType stype)
30247c6ae99SBarry Smith {
30347c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
30447c6ae99SBarry Smith 
30547c6ae99SBarry Smith   PetscFunctionBegin;
30647c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
30747c6ae99SBarry Smith   PetscValidLogicalCollectiveEnum(da,stype,2);
308cb630486SJed Brown   if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
30947c6ae99SBarry Smith   dd->stencil_type = stype;
31047c6ae99SBarry Smith   PetscFunctionReturn(0);
31147c6ae99SBarry Smith }
31247c6ae99SBarry Smith 
31347c6ae99SBarry Smith #undef __FUNCT__
314aa219208SBarry Smith #define __FUNCT__ "DMDASetStencilWidth"
31547c6ae99SBarry Smith /*@
316aa219208SBarry Smith   DMDASetStencilWidth - Sets the width of the communication stencil
31747c6ae99SBarry Smith 
318aa219208SBarry Smith   Logically Collective on DMDA
31947c6ae99SBarry Smith 
32047c6ae99SBarry Smith   Input Parameter:
321aa219208SBarry Smith + da    - The DMDA
32247c6ae99SBarry Smith - width - The stencil width
32347c6ae99SBarry Smith 
32447c6ae99SBarry Smith   Level: intermediate
32547c6ae99SBarry Smith 
32647c6ae99SBarry Smith .keywords:  distributed array, stencil
32796e147daSBarry Smith .seealso: DMDACreate(), DMDestroy(), DMDA
32847c6ae99SBarry Smith @*/
3297087cfbeSBarry Smith PetscErrorCode  DMDASetStencilWidth(DM da, PetscInt width)
33047c6ae99SBarry Smith {
33147c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
33247c6ae99SBarry Smith 
33347c6ae99SBarry Smith   PetscFunctionBegin;
33447c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
33547c6ae99SBarry Smith   PetscValidLogicalCollectiveInt(da,width,2);
336cb630486SJed Brown   if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
33747c6ae99SBarry Smith   dd->s = width;
33847c6ae99SBarry Smith   PetscFunctionReturn(0);
33947c6ae99SBarry Smith }
34047c6ae99SBarry Smith 
34147c6ae99SBarry Smith #undef __FUNCT__
342aa219208SBarry Smith #define __FUNCT__ "DMDACheckOwnershipRanges_Private"
343aa219208SBarry Smith static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
34447c6ae99SBarry Smith {
34547c6ae99SBarry Smith   PetscInt i,sum;
34647c6ae99SBarry Smith 
34747c6ae99SBarry Smith   PetscFunctionBegin;
34847c6ae99SBarry Smith   if (M < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"Global dimension not set");
34947c6ae99SBarry Smith   for (i=sum=0; i<m; i++) sum += lx[i];
35047c6ae99SBarry Smith   if (sum != M) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"Ownership ranges sum to %D but global dimension is %D",sum,M);
35147c6ae99SBarry Smith   PetscFunctionReturn(0);
35247c6ae99SBarry Smith }
35347c6ae99SBarry Smith 
35447c6ae99SBarry Smith #undef __FUNCT__
355aa219208SBarry Smith #define __FUNCT__ "DMDASetOwnershipRanges"
35647c6ae99SBarry Smith /*@
357aa219208SBarry Smith   DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
35847c6ae99SBarry Smith 
359aa219208SBarry Smith   Logically Collective on DMDA
36047c6ae99SBarry Smith 
36147c6ae99SBarry Smith   Input Parameter:
362aa219208SBarry Smith + da - The DMDA
36347c6ae99SBarry Smith . lx - array containing number of nodes in the X direction on each process, or PETSC_NULL. If non-null, must be of length da->m
36447c6ae99SBarry Smith . ly - array containing number of nodes in the Y direction on each process, or PETSC_NULL. If non-null, must be of length da->n
36547c6ae99SBarry Smith - lz - array containing number of nodes in the Z direction on each process, or PETSC_NULL. If non-null, must be of length da->p.
36647c6ae99SBarry Smith 
36747c6ae99SBarry Smith   Level: intermediate
36847c6ae99SBarry Smith 
369e3f3e4b6SBarry Smith   Note: these numbers are NOT multiplied by the number of dof per node.
370e3f3e4b6SBarry Smith 
37147c6ae99SBarry Smith .keywords:  distributed array
37296e147daSBarry Smith .seealso: DMDACreate(), DMDestroy(), DMDA
37347c6ae99SBarry Smith @*/
3747087cfbeSBarry Smith PetscErrorCode  DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
37547c6ae99SBarry Smith {
37647c6ae99SBarry Smith   PetscErrorCode ierr;
37747c6ae99SBarry Smith   DM_DA          *dd = (DM_DA*)da->data;
37847c6ae99SBarry Smith 
37947c6ae99SBarry Smith   PetscFunctionBegin;
38047c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
381cb630486SJed Brown   if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
38247c6ae99SBarry Smith   if (lx) {
38347c6ae99SBarry Smith     if (dd->m < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
384aa219208SBarry Smith     ierr = DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);CHKERRQ(ierr);
38547c6ae99SBarry Smith     if (!dd->lx) {
38647c6ae99SBarry Smith       ierr = PetscMalloc(dd->m*sizeof(PetscInt), &dd->lx);CHKERRQ(ierr);
38747c6ae99SBarry Smith     }
38847c6ae99SBarry Smith     ierr = PetscMemcpy(dd->lx, lx, dd->m*sizeof(PetscInt));CHKERRQ(ierr);
38947c6ae99SBarry Smith   }
39047c6ae99SBarry Smith   if (ly) {
39147c6ae99SBarry Smith     if (dd->n < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
392aa219208SBarry Smith     ierr = DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);CHKERRQ(ierr);
39347c6ae99SBarry Smith     if (!dd->ly) {
39447c6ae99SBarry Smith       ierr = PetscMalloc(dd->n*sizeof(PetscInt), &dd->ly);CHKERRQ(ierr);
39547c6ae99SBarry Smith     }
39647c6ae99SBarry Smith     ierr = PetscMemcpy(dd->ly, ly, dd->n*sizeof(PetscInt));CHKERRQ(ierr);
39747c6ae99SBarry Smith   }
39847c6ae99SBarry Smith   if (lz) {
39947c6ae99SBarry Smith     if (dd->p < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
400aa219208SBarry Smith     ierr = DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);CHKERRQ(ierr);
40147c6ae99SBarry Smith     if (!dd->lz) {
40247c6ae99SBarry Smith       ierr = PetscMalloc(dd->p*sizeof(PetscInt), &dd->lz);CHKERRQ(ierr);
40347c6ae99SBarry Smith     }
40447c6ae99SBarry Smith     ierr = PetscMemcpy(dd->lz, lz, dd->p*sizeof(PetscInt));CHKERRQ(ierr);
40547c6ae99SBarry Smith   }
40647c6ae99SBarry Smith   PetscFunctionReturn(0);
40747c6ae99SBarry Smith }
40847c6ae99SBarry Smith 
40947c6ae99SBarry Smith #undef __FUNCT__
410aa219208SBarry Smith #define __FUNCT__ "DMDASetInterpolationType"
41147c6ae99SBarry Smith /*@
412aa219208SBarry Smith        DMDASetInterpolationType - Sets the type of interpolation that will be
413e727c939SJed Brown           returned by DMCreateInterpolation()
41447c6ae99SBarry Smith 
415aa219208SBarry Smith    Logically Collective on DMDA
41647c6ae99SBarry Smith 
41747c6ae99SBarry Smith    Input Parameter:
41847c6ae99SBarry Smith +  da - initial distributed array
419aa219208SBarry Smith .  ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms
42047c6ae99SBarry Smith 
42147c6ae99SBarry Smith    Level: intermediate
42247c6ae99SBarry Smith 
423e727c939SJed Brown    Notes: you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()
42447c6ae99SBarry Smith 
42547c6ae99SBarry Smith .keywords:  distributed array, interpolation
42647c6ae99SBarry Smith 
42796e147daSBarry Smith .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
42847c6ae99SBarry Smith @*/
4297087cfbeSBarry Smith PetscErrorCode  DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
43047c6ae99SBarry Smith {
43147c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
43247c6ae99SBarry Smith 
43347c6ae99SBarry Smith   PetscFunctionBegin;
43447c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
43547c6ae99SBarry Smith   PetscValidLogicalCollectiveEnum(da,ctype,2);
43647c6ae99SBarry Smith   dd->interptype = ctype;
43747c6ae99SBarry Smith   PetscFunctionReturn(0);
43847c6ae99SBarry Smith }
43947c6ae99SBarry Smith 
4402dde6fd4SLisandro Dalcin #undef __FUNCT__
4412dde6fd4SLisandro Dalcin #define __FUNCT__ "DMDAGetInterpolationType"
4422dde6fd4SLisandro Dalcin /*@
4432dde6fd4SLisandro Dalcin        DMDAGetInterpolationType - Gets the type of interpolation that will be
444e727c939SJed Brown           used by DMCreateInterpolation()
4452dde6fd4SLisandro Dalcin 
4462dde6fd4SLisandro Dalcin    Not Collective
4472dde6fd4SLisandro Dalcin 
4482dde6fd4SLisandro Dalcin    Input Parameter:
4492dde6fd4SLisandro Dalcin .  da - distributed array
4502dde6fd4SLisandro Dalcin 
4512dde6fd4SLisandro Dalcin    Output Parameter:
4522dde6fd4SLisandro Dalcin .  ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)
4532dde6fd4SLisandro Dalcin 
4542dde6fd4SLisandro Dalcin    Level: intermediate
4552dde6fd4SLisandro Dalcin 
4562dde6fd4SLisandro Dalcin .keywords:  distributed array, interpolation
4572dde6fd4SLisandro Dalcin 
458e727c939SJed Brown .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
4592dde6fd4SLisandro Dalcin @*/
4602dde6fd4SLisandro Dalcin PetscErrorCode  DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
4612dde6fd4SLisandro Dalcin {
4622dde6fd4SLisandro Dalcin   DM_DA *dd = (DM_DA*)da->data;
4632dde6fd4SLisandro Dalcin 
4642dde6fd4SLisandro Dalcin   PetscFunctionBegin;
4652dde6fd4SLisandro Dalcin   PetscValidHeaderSpecific(da,DM_CLASSID,1);
4662dde6fd4SLisandro Dalcin   PetscValidPointer(ctype,2);
4672dde6fd4SLisandro Dalcin   *ctype = dd->interptype;
4682dde6fd4SLisandro Dalcin   PetscFunctionReturn(0);
4692dde6fd4SLisandro Dalcin }
47047c6ae99SBarry Smith 
47147c6ae99SBarry Smith #undef __FUNCT__
472aa219208SBarry Smith #define __FUNCT__ "DMDAGetNeighbors"
47347c6ae99SBarry Smith /*@C
474aa219208SBarry Smith       DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
47547c6ae99SBarry Smith         processes neighbors.
47647c6ae99SBarry Smith 
47747c6ae99SBarry Smith     Not Collective
47847c6ae99SBarry Smith 
47947c6ae99SBarry Smith    Input Parameter:
480aa219208SBarry Smith .     da - the DMDA object
48147c6ae99SBarry Smith 
48247c6ae99SBarry Smith    Output Parameters:
48347c6ae99SBarry Smith .     ranks - the neighbors ranks, stored with the x index increasing most rapidly.
48447c6ae99SBarry Smith               this process itself is in the list
48547c6ae99SBarry Smith 
48647c6ae99SBarry Smith    Notes: In 2d the array is of length 9, in 3d of length 27
48747c6ae99SBarry Smith           Not supported in 1d
488aa219208SBarry Smith           Do not free the array, it is freed when the DMDA is destroyed.
48947c6ae99SBarry Smith 
49047c6ae99SBarry Smith    Fortran Notes: In fortran you must pass in an array of the appropriate length.
49147c6ae99SBarry Smith 
49247c6ae99SBarry Smith    Level: intermediate
49347c6ae99SBarry Smith 
49447c6ae99SBarry Smith @*/
4957087cfbeSBarry Smith PetscErrorCode  DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
49647c6ae99SBarry Smith {
49747c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
4985fd66863SKarl Rupp 
49947c6ae99SBarry Smith   PetscFunctionBegin;
50047c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
50147c6ae99SBarry Smith   *ranks = dd->neighbors;
50247c6ae99SBarry Smith   PetscFunctionReturn(0);
50347c6ae99SBarry Smith }
50447c6ae99SBarry Smith 
50547c6ae99SBarry Smith #undef __FUNCT__
506aa219208SBarry Smith #define __FUNCT__ "DMDAGetOwnershipRanges"
50747c6ae99SBarry Smith /*@C
508aa219208SBarry Smith       DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process
50947c6ae99SBarry Smith 
51047c6ae99SBarry Smith     Not Collective
51147c6ae99SBarry Smith 
51247c6ae99SBarry Smith    Input Parameter:
513aa219208SBarry Smith .     da - the DMDA object
51447c6ae99SBarry Smith 
51547c6ae99SBarry Smith    Output Parameter:
51647c6ae99SBarry Smith +     lx - ownership along x direction (optional)
51747c6ae99SBarry Smith .     ly - ownership along y direction (optional)
51847c6ae99SBarry Smith -     lz - ownership along z direction (optional)
51947c6ae99SBarry Smith 
52047c6ae99SBarry Smith    Level: intermediate
52147c6ae99SBarry Smith 
522aa219208SBarry Smith     Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()
52347c6ae99SBarry Smith 
52447c6ae99SBarry Smith     In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
525aa219208SBarry Smith     eighth arguments from DMDAGetInfo()
52647c6ae99SBarry Smith 
52747c6ae99SBarry Smith      In C you should not free these arrays, nor change the values in them. They will only have valid values while the
528aa219208SBarry Smith     DMDA they came from still exists (has not been destroyed).
52947c6ae99SBarry Smith 
530e3f3e4b6SBarry Smith     These numbers are NOT multiplied by the number of dof per node.
531e3f3e4b6SBarry Smith 
532aa219208SBarry Smith .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
53347c6ae99SBarry Smith @*/
5347087cfbeSBarry Smith PetscErrorCode  DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
53547c6ae99SBarry Smith {
53647c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
53747c6ae99SBarry Smith 
53847c6ae99SBarry Smith   PetscFunctionBegin;
53947c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
54047c6ae99SBarry Smith   if (lx) *lx = dd->lx;
54147c6ae99SBarry Smith   if (ly) *ly = dd->ly;
54247c6ae99SBarry Smith   if (lz) *lz = dd->lz;
54347c6ae99SBarry Smith   PetscFunctionReturn(0);
54447c6ae99SBarry Smith }
54547c6ae99SBarry Smith 
54647c6ae99SBarry Smith #undef __FUNCT__
547aa219208SBarry Smith #define __FUNCT__ "DMDASetRefinementFactor"
54847c6ae99SBarry Smith /*@
549aa219208SBarry Smith      DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined
55047c6ae99SBarry Smith 
551aa219208SBarry Smith     Logically Collective on DMDA
55247c6ae99SBarry Smith 
55347c6ae99SBarry Smith   Input Parameters:
554aa219208SBarry Smith +    da - the DMDA object
55547c6ae99SBarry Smith .    refine_x - ratio of fine grid to coarse in x direction (2 by default)
55647c6ae99SBarry Smith .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
55747c6ae99SBarry Smith -    refine_z - ratio of fine grid to coarse in z direction (2 by default)
55847c6ae99SBarry Smith 
55947c6ae99SBarry Smith   Options Database:
56047c6ae99SBarry Smith +  -da_refine_x - refinement ratio in x direction
56147c6ae99SBarry Smith .  -da_refine_y - refinement ratio in y direction
56247c6ae99SBarry Smith -  -da_refine_z - refinement ratio in z direction
56347c6ae99SBarry Smith 
56447c6ae99SBarry Smith   Level: intermediate
56547c6ae99SBarry Smith 
56647c6ae99SBarry Smith     Notes: Pass PETSC_IGNORE to leave a value unchanged
56747c6ae99SBarry Smith 
568aa219208SBarry Smith .seealso: DMRefine(), DMDAGetRefinementFactor()
56947c6ae99SBarry Smith @*/
5707087cfbeSBarry Smith PetscErrorCode  DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
57147c6ae99SBarry Smith {
57247c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
57347c6ae99SBarry Smith 
57447c6ae99SBarry Smith   PetscFunctionBegin;
57547c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
57647c6ae99SBarry Smith   PetscValidLogicalCollectiveInt(da,refine_x,2);
57747c6ae99SBarry Smith   PetscValidLogicalCollectiveInt(da,refine_y,3);
57847c6ae99SBarry Smith   PetscValidLogicalCollectiveInt(da,refine_z,4);
57947c6ae99SBarry Smith 
58047c6ae99SBarry Smith   if (refine_x > 0) dd->refine_x = refine_x;
58147c6ae99SBarry Smith   if (refine_y > 0) dd->refine_y = refine_y;
58247c6ae99SBarry Smith   if (refine_z > 0) dd->refine_z = refine_z;
58347c6ae99SBarry Smith   PetscFunctionReturn(0);
58447c6ae99SBarry Smith }
58547c6ae99SBarry Smith 
58647c6ae99SBarry Smith #undef __FUNCT__
587aa219208SBarry Smith #define __FUNCT__ "DMDAGetRefinementFactor"
58847c6ae99SBarry Smith /*@C
589aa219208SBarry Smith      DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined
59047c6ae99SBarry Smith 
59147c6ae99SBarry Smith     Not Collective
59247c6ae99SBarry Smith 
59347c6ae99SBarry Smith   Input Parameter:
594aa219208SBarry Smith .    da - the DMDA object
59547c6ae99SBarry Smith 
59647c6ae99SBarry Smith   Output Parameters:
59747c6ae99SBarry Smith +    refine_x - ratio of fine grid to coarse in x direction (2 by default)
59847c6ae99SBarry Smith .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
59947c6ae99SBarry Smith -    refine_z - ratio of fine grid to coarse in z direction (2 by default)
60047c6ae99SBarry Smith 
60147c6ae99SBarry Smith   Level: intermediate
60247c6ae99SBarry Smith 
60347c6ae99SBarry Smith     Notes: Pass PETSC_NULL for values you do not need
60447c6ae99SBarry Smith 
605aa219208SBarry Smith .seealso: DMRefine(), DMDASetRefinementFactor()
60647c6ae99SBarry Smith @*/
6077087cfbeSBarry Smith PetscErrorCode  DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
60847c6ae99SBarry Smith {
60947c6ae99SBarry Smith   DM_DA *dd = (DM_DA*)da->data;
61047c6ae99SBarry Smith 
61147c6ae99SBarry Smith   PetscFunctionBegin;
61247c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
61347c6ae99SBarry Smith   if (refine_x) *refine_x = dd->refine_x;
61447c6ae99SBarry Smith   if (refine_y) *refine_y = dd->refine_y;
61547c6ae99SBarry Smith   if (refine_z) *refine_z = dd->refine_z;
61647c6ae99SBarry Smith   PetscFunctionReturn(0);
61747c6ae99SBarry Smith }
61847c6ae99SBarry Smith 
61947c6ae99SBarry Smith #undef __FUNCT__
620aa219208SBarry Smith #define __FUNCT__ "DMDASetGetMatrix"
62147c6ae99SBarry Smith /*@C
622aa219208SBarry Smith      DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.
62347c6ae99SBarry Smith 
624aa219208SBarry Smith     Logically Collective on DMDA
62547c6ae99SBarry Smith 
62647c6ae99SBarry Smith   Input Parameters:
627aa219208SBarry Smith +    da - the DMDA object
628aa219208SBarry Smith -    f - the function that allocates the matrix for that specific DMDA
62947c6ae99SBarry Smith 
63047c6ae99SBarry Smith   Level: developer
63147c6ae99SBarry Smith 
632aa219208SBarry Smith    Notes: See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
63347c6ae99SBarry Smith        the diagonal and off-diagonal blocks of the matrix
63447c6ae99SBarry Smith 
635950540a4SJed Brown .seealso: DMCreateMatrix(), DMDASetBlockFills()
63647c6ae99SBarry Smith @*/
63719fd82e9SBarry Smith PetscErrorCode  DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, MatType,Mat*))
63847c6ae99SBarry Smith {
63947c6ae99SBarry Smith   PetscFunctionBegin;
64047c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
64125296bd5SBarry Smith   da->ops->creatematrix = f;
64247c6ae99SBarry Smith   PetscFunctionReturn(0);
64347c6ae99SBarry Smith }
64447c6ae99SBarry Smith 
64547c6ae99SBarry Smith #undef __FUNCT__
64694013140SBarry Smith #define __FUNCT__ "DMDARefineOwnershipRanges"
64747c6ae99SBarry Smith /*
64847c6ae99SBarry Smith   Creates "balanced" ownership ranges after refinement, constrained by the need for the
64947c6ae99SBarry Smith   fine grid boundaries to fall within one stencil width of the coarse partition.
65047c6ae99SBarry Smith 
65147c6ae99SBarry Smith   Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
65247c6ae99SBarry Smith */
65394013140SBarry Smith static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
65447c6ae99SBarry Smith {
65547c6ae99SBarry Smith   PetscInt       i,totalc = 0,remaining,startc = 0,startf = 0;
65647c6ae99SBarry Smith   PetscErrorCode ierr;
65747c6ae99SBarry Smith 
65847c6ae99SBarry Smith   PetscFunctionBegin;
65947c6ae99SBarry Smith   if (ratio < 1) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
66047c6ae99SBarry Smith   if (ratio == 1) {
66147c6ae99SBarry Smith     ierr = PetscMemcpy(lf,lc,m*sizeof(lc[0]));CHKERRQ(ierr);
66247c6ae99SBarry Smith     PetscFunctionReturn(0);
66347c6ae99SBarry Smith   }
66447c6ae99SBarry Smith   for (i=0; i<m; i++) totalc += lc[i];
66547c6ae99SBarry Smith   remaining = (!periodic) + ratio * (totalc - (!periodic));
66647c6ae99SBarry Smith   for (i=0; i<m; i++) {
66747c6ae99SBarry Smith     PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
66847c6ae99SBarry Smith     if (i == m-1) lf[i] = want;
66947c6ae99SBarry Smith     else {
6707aca7175SJed Brown       const PetscInt nextc = startc + lc[i];
6717aca7175SJed Brown       /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
6727aca7175SJed Brown        * coarse stencil width of the first coarse node in the next subdomain. */
6737aca7175SJed Brown       while ((startf+want)/ratio < nextc - stencil_width) want++;
6747aca7175SJed Brown       /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
6757aca7175SJed Brown        * coarse stencil width of the last coarse node in the current subdomain. */
6767aca7175SJed Brown       while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
6777aca7175SJed Brown       /* Make sure all constraints are satisfied */
67830729d88SBarry Smith       if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
67930729d88SBarry Smith           || ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Could not find a compatible refined ownership range");
68047c6ae99SBarry Smith     }
68147c6ae99SBarry Smith     lf[i]      = want;
68247c6ae99SBarry Smith     startc    += lc[i];
68347c6ae99SBarry Smith     startf    += lf[i];
68447c6ae99SBarry Smith     remaining -= lf[i];
68547c6ae99SBarry Smith   }
68647c6ae99SBarry Smith   PetscFunctionReturn(0);
68747c6ae99SBarry Smith }
68847c6ae99SBarry Smith 
68947c6ae99SBarry Smith #undef __FUNCT__
6902be375d4SJed Brown #define __FUNCT__ "DMDACoarsenOwnershipRanges"
6912be375d4SJed Brown /*
6922be375d4SJed Brown   Creates "balanced" ownership ranges after coarsening, constrained by the need for the
6932be375d4SJed Brown   fine grid boundaries to fall within one stencil width of the coarse partition.
6942be375d4SJed Brown 
6952be375d4SJed Brown   Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
6962be375d4SJed Brown */
6972be375d4SJed Brown static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
6982be375d4SJed Brown {
6992be375d4SJed Brown   PetscInt       i,totalf,remaining,startc,startf;
7002be375d4SJed Brown   PetscErrorCode ierr;
7012be375d4SJed Brown 
7022be375d4SJed Brown   PetscFunctionBegin;
7032be375d4SJed Brown   if (ratio < 1) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
7042be375d4SJed Brown   if (ratio == 1) {
7052be375d4SJed Brown     ierr = PetscMemcpy(lc,lf,m*sizeof(lf[0]));CHKERRQ(ierr);
7062be375d4SJed Brown     PetscFunctionReturn(0);
7072be375d4SJed Brown   }
7082be375d4SJed Brown   for (i=0,totalf=0; i<m; i++) totalf += lf[i];
7092be375d4SJed Brown   remaining = (!periodic) + (totalf - (!periodic)) / ratio;
7102be375d4SJed Brown   for (i=0,startc=0,startf=0; i<m; i++) {
7112be375d4SJed Brown     PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
7122be375d4SJed Brown     if (i == m-1) lc[i] = want;
7132be375d4SJed Brown     else {
7142be375d4SJed Brown       const PetscInt nextf = startf+lf[i];
7152be375d4SJed Brown       /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
7162be375d4SJed Brown        * node is within one stencil width. */
7172be375d4SJed Brown       while (nextf/ratio < startc+want-stencil_width) want--;
7182be375d4SJed Brown       /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
7192be375d4SJed Brown        * fine node is within one stencil width. */
7202be375d4SJed Brown       while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
7212be375d4SJed Brown       if (want < 0 || want > remaining
72230729d88SBarry Smith           || (nextf/ratio < startc+want-stencil_width) || ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Could not find a compatible coarsened ownership range");
7232be375d4SJed Brown     }
7242be375d4SJed Brown     lc[i]      = want;
7252be375d4SJed Brown     startc    += lc[i];
7262be375d4SJed Brown     startf    += lf[i];
7272be375d4SJed Brown     remaining -= lc[i];
7282be375d4SJed Brown   }
7292be375d4SJed Brown   PetscFunctionReturn(0);
7302be375d4SJed Brown }
7312be375d4SJed Brown 
7322be375d4SJed Brown #undef __FUNCT__
73394013140SBarry Smith #define __FUNCT__ "DMRefine_DA"
7347087cfbeSBarry Smith PetscErrorCode  DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
73547c6ae99SBarry Smith {
73647c6ae99SBarry Smith   PetscErrorCode ierr;
737f3141302SJed Brown   PetscInt       M,N,P,i;
7389a42bb27SBarry Smith   DM             da2;
73947c6ae99SBarry Smith   DM_DA          *dd = (DM_DA*)da->data,*dd2;
74047c6ae99SBarry Smith 
74147c6ae99SBarry Smith   PetscFunctionBegin;
74247c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
74347c6ae99SBarry Smith   PetscValidPointer(daref,3);
74447c6ae99SBarry Smith 
7451321219cSEthan Coon   if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
74647c6ae99SBarry Smith     M = dd->refine_x*dd->M;
74747c6ae99SBarry Smith   } else {
74847c6ae99SBarry Smith     M = 1 + dd->refine_x*(dd->M - 1);
74947c6ae99SBarry Smith   }
7501321219cSEthan Coon   if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
7511860e6e9SBarry Smith     if (dd->dim > 1) {
75247c6ae99SBarry Smith       N = dd->refine_y*dd->N;
75347c6ae99SBarry Smith     } else {
7541860e6e9SBarry Smith       N = 1;
7551860e6e9SBarry Smith     }
7561860e6e9SBarry Smith   } else {
75747c6ae99SBarry Smith     N = 1 + dd->refine_y*(dd->N - 1);
75847c6ae99SBarry Smith   }
7591321219cSEthan Coon   if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
7601860e6e9SBarry Smith     if (dd->dim > 2) {
76147c6ae99SBarry Smith       P = dd->refine_z*dd->P;
76247c6ae99SBarry Smith     } else {
7631860e6e9SBarry Smith       P = 1;
7641860e6e9SBarry Smith     }
7651860e6e9SBarry Smith   } else {
76647c6ae99SBarry Smith     P = 1 + dd->refine_z*(dd->P - 1);
76747c6ae99SBarry Smith   }
768397b6216SJed Brown   ierr = DMDACreate(((PetscObject)da)->comm,&da2);CHKERRQ(ierr);
769397b6216SJed Brown   ierr = DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);CHKERRQ(ierr);
770397b6216SJed Brown   ierr = DMDASetDim(da2,dd->dim);CHKERRQ(ierr);
771397b6216SJed Brown   ierr = DMDASetSizes(da2,M,N,P);CHKERRQ(ierr);
772397b6216SJed Brown   ierr = DMDASetNumProcs(da2,dd->m,dd->n,dd->p);CHKERRQ(ierr);
773397b6216SJed Brown   ierr = DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);CHKERRQ(ierr);
774397b6216SJed Brown   ierr = DMDASetDof(da2,dd->w);CHKERRQ(ierr);
775397b6216SJed Brown   ierr = DMDASetStencilType(da2,dd->stencil_type);CHKERRQ(ierr);
776397b6216SJed Brown   ierr = DMDASetStencilWidth(da2,dd->s);CHKERRQ(ierr);
77747c6ae99SBarry Smith   if (dd->dim == 3) {
77847c6ae99SBarry Smith     PetscInt *lx,*ly,*lz;
77947c6ae99SBarry Smith     ierr = PetscMalloc3(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly,dd->p,PetscInt,&lz);CHKERRQ(ierr);
7801321219cSEthan Coon     ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
7811321219cSEthan Coon     ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);CHKERRQ(ierr);
7821321219cSEthan Coon     ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);CHKERRQ(ierr);
783397b6216SJed Brown     ierr = DMDASetOwnershipRanges(da2,lx,ly,lz);CHKERRQ(ierr);
78447c6ae99SBarry Smith     ierr = PetscFree3(lx,ly,lz);CHKERRQ(ierr);
78547c6ae99SBarry Smith   } else if (dd->dim == 2) {
78647c6ae99SBarry Smith     PetscInt *lx,*ly;
78747c6ae99SBarry Smith     ierr = PetscMalloc2(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly);CHKERRQ(ierr);
7881321219cSEthan Coon     ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
7891321219cSEthan Coon     ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);CHKERRQ(ierr);
790397b6216SJed Brown     ierr = DMDASetOwnershipRanges(da2,lx,ly,PETSC_NULL);CHKERRQ(ierr);
79147c6ae99SBarry Smith     ierr = PetscFree2(lx,ly);CHKERRQ(ierr);
79247c6ae99SBarry Smith   } else if (dd->dim == 1) {
79347c6ae99SBarry Smith     PetscInt *lx;
79447c6ae99SBarry Smith     ierr = PetscMalloc(dd->m*sizeof(PetscInt),&lx);CHKERRQ(ierr);
7951321219cSEthan Coon     ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
796397b6216SJed Brown     ierr = DMDASetOwnershipRanges(da2,lx,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
79747c6ae99SBarry Smith     ierr = PetscFree(lx);CHKERRQ(ierr);
79847c6ae99SBarry Smith   }
79947c6ae99SBarry Smith   dd2 = (DM_DA*)da2->data;
80047c6ae99SBarry Smith 
80147c6ae99SBarry Smith   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
80225296bd5SBarry Smith   da2->ops->creatematrix = da->ops->creatematrix;
80325296bd5SBarry Smith   /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
80447c6ae99SBarry Smith   da2->ops->getcoloring = da->ops->getcoloring;
80547c6ae99SBarry Smith   dd2->interptype       = dd->interptype;
80647c6ae99SBarry Smith 
80747c6ae99SBarry Smith   /* copy fill information if given */
80847c6ae99SBarry Smith   if (dd->dfill) {
80947c6ae99SBarry Smith     ierr = PetscMalloc((dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->dfill);CHKERRQ(ierr);
81047c6ae99SBarry Smith     ierr = PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));CHKERRQ(ierr);
81147c6ae99SBarry Smith   }
81247c6ae99SBarry Smith   if (dd->ofill) {
81347c6ae99SBarry Smith     ierr = PetscMalloc((dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->ofill);CHKERRQ(ierr);
81447c6ae99SBarry Smith     ierr = PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));CHKERRQ(ierr);
81547c6ae99SBarry Smith   }
81647c6ae99SBarry Smith   /* copy the refine information */
817397b6216SJed Brown   dd2->coarsen_x = dd2->refine_x = dd->refine_x;
818397b6216SJed Brown   dd2->coarsen_y = dd2->refine_y = dd->refine_y;
819397b6216SJed Brown   dd2->coarsen_z = dd2->refine_z = dd->refine_z;
82047c6ae99SBarry Smith 
82147c6ae99SBarry Smith   /* copy vector type information */
82247c6ae99SBarry Smith   ierr = PetscFree(da2->vectype);CHKERRQ(ierr);
82319fd82e9SBarry Smith   ierr = PetscStrallocpy(da->vectype,(char**)&da2->vectype);CHKERRQ(ierr);
824ddcf8b74SDave May 
825efd51863SBarry Smith   dd2->lf = dd->lf;
826efd51863SBarry Smith   dd2->lj = dd->lj;
827efd51863SBarry Smith 
8286e87535bSJed Brown   da2->leveldown = da->leveldown;
8296e87535bSJed Brown   da2->levelup   = da->levelup + 1;
830*8865f1eaSKarl Rupp 
8316e87535bSJed Brown   ierr = DMSetFromOptions(da2);CHKERRQ(ierr);
8326e87535bSJed Brown   ierr = DMSetUp(da2);CHKERRQ(ierr);
833ca266f36SBarry Smith   ierr = DMViewFromOptions(da2,"-dm_view");CHKERRQ(ierr);
8346e87535bSJed Brown 
835ddcf8b74SDave May   /* interpolate coordinates if they are set on the coarse grid */
8366636e97aSMatthew G Knepley   if (da->coordinates) {
837ddcf8b74SDave May     DM  cdaf,cdac;
838ddcf8b74SDave May     Vec coordsc,coordsf;
839ddcf8b74SDave May     Mat II;
840ddcf8b74SDave May 
8416636e97aSMatthew G Knepley     ierr = DMGetCoordinateDM(da,&cdac);CHKERRQ(ierr);
8426636e97aSMatthew G Knepley     ierr = DMGetCoordinates(da,&coordsc);CHKERRQ(ierr);
8436636e97aSMatthew G Knepley     ierr = DMGetCoordinateDM(da2,&cdaf);CHKERRQ(ierr);
844b61d3410SDave May     /* force creation of the coordinate vector */
845b61d3410SDave May     ierr = DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);
8466636e97aSMatthew G Knepley     ierr = DMGetCoordinates(da2,&coordsf);CHKERRQ(ierr);
847e727c939SJed Brown     ierr = DMCreateInterpolation(cdac,cdaf,&II,PETSC_NULL);CHKERRQ(ierr);
848ddcf8b74SDave May     ierr = MatInterpolate(II,coordsc,coordsf);CHKERRQ(ierr);
849fcfd50ebSBarry Smith     ierr = MatDestroy(&II);CHKERRQ(ierr);
850ddcf8b74SDave May   }
851397b6216SJed Brown 
852f3141302SJed Brown   for (i=0; i<da->bs; i++) {
853f3141302SJed Brown     const char *fieldname;
854f3141302SJed Brown     ierr = DMDAGetFieldName(da,i,&fieldname);CHKERRQ(ierr);
855f3141302SJed Brown     ierr = DMDASetFieldName(da2,i,fieldname);CHKERRQ(ierr);
856f3141302SJed Brown   }
857397b6216SJed Brown 
85847c6ae99SBarry Smith   *daref = da2;
85947c6ae99SBarry Smith   PetscFunctionReturn(0);
86047c6ae99SBarry Smith }
86147c6ae99SBarry Smith 
862397b6216SJed Brown 
86347c6ae99SBarry Smith #undef __FUNCT__
86494013140SBarry Smith #define __FUNCT__ "DMCoarsen_DA"
8657087cfbeSBarry Smith PetscErrorCode  DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
86647c6ae99SBarry Smith {
86747c6ae99SBarry Smith   PetscErrorCode ierr;
868397b6216SJed Brown   PetscInt       M,N,P,i;
8699a42bb27SBarry Smith   DM             da2;
87047c6ae99SBarry Smith   DM_DA          *dd = (DM_DA*)da->data,*dd2;
87147c6ae99SBarry Smith 
87247c6ae99SBarry Smith   PetscFunctionBegin;
87347c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
87447c6ae99SBarry Smith   PetscValidPointer(daref,3);
87547c6ae99SBarry Smith 
8761321219cSEthan Coon   if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
877397b6216SJed Brown     M = dd->M / dd->coarsen_x;
87847c6ae99SBarry Smith   } else {
879397b6216SJed Brown     M = 1 + (dd->M - 1) / dd->coarsen_x;
88047c6ae99SBarry Smith   }
8811321219cSEthan Coon   if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
8821860e6e9SBarry Smith     if (dd->dim > 1) {
883397b6216SJed Brown       N = dd->N / dd->coarsen_y;
88447c6ae99SBarry Smith     } else {
8851860e6e9SBarry Smith       N = 1;
8861860e6e9SBarry Smith     }
8871860e6e9SBarry Smith   } else {
888397b6216SJed Brown     N = 1 + (dd->N - 1) / dd->coarsen_y;
88947c6ae99SBarry Smith   }
8901321219cSEthan Coon   if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
8911860e6e9SBarry Smith     if (dd->dim > 2) {
892397b6216SJed Brown       P = dd->P / dd->coarsen_z;
89347c6ae99SBarry Smith     } else {
8941860e6e9SBarry Smith       P = 1;
8951860e6e9SBarry Smith     }
8961860e6e9SBarry Smith   } else {
897397b6216SJed Brown     P = 1 + (dd->P - 1) / dd->coarsen_z;
89847c6ae99SBarry Smith   }
899397b6216SJed Brown   ierr = DMDACreate(((PetscObject)da)->comm,&da2);CHKERRQ(ierr);
900397b6216SJed Brown   ierr = DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);CHKERRQ(ierr);
901397b6216SJed Brown   ierr = DMDASetDim(da2,dd->dim);CHKERRQ(ierr);
902397b6216SJed Brown   ierr = DMDASetSizes(da2,M,N,P);CHKERRQ(ierr);
903397b6216SJed Brown   ierr = DMDASetNumProcs(da2,dd->m,dd->n,dd->p);CHKERRQ(ierr);
904397b6216SJed Brown   ierr = DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);CHKERRQ(ierr);
905397b6216SJed Brown   ierr = DMDASetDof(da2,dd->w);CHKERRQ(ierr);
906397b6216SJed Brown   ierr = DMDASetStencilType(da2,dd->stencil_type);CHKERRQ(ierr);
907397b6216SJed Brown   ierr = DMDASetStencilWidth(da2,dd->s);CHKERRQ(ierr);
90847c6ae99SBarry Smith   if (dd->dim == 3) {
9092be375d4SJed Brown     PetscInt *lx,*ly,*lz;
9102be375d4SJed Brown     ierr = PetscMalloc3(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly,dd->p,PetscInt,&lz);CHKERRQ(ierr);
911397b6216SJed Brown     ierr = DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
912397b6216SJed Brown     ierr = DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);CHKERRQ(ierr);
913397b6216SJed Brown     ierr = DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);CHKERRQ(ierr);
914397b6216SJed Brown     ierr = DMDASetOwnershipRanges(da2,lx,ly,lz);CHKERRQ(ierr);
9152be375d4SJed Brown     ierr = PetscFree3(lx,ly,lz);CHKERRQ(ierr);
91647c6ae99SBarry Smith   } else if (dd->dim == 2) {
9172be375d4SJed Brown     PetscInt *lx,*ly;
9182be375d4SJed Brown     ierr = PetscMalloc2(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly);CHKERRQ(ierr);
919397b6216SJed Brown     ierr = DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
920397b6216SJed Brown     ierr = DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);CHKERRQ(ierr);
921397b6216SJed Brown     ierr = DMDASetOwnershipRanges(da2,lx,ly,PETSC_NULL);CHKERRQ(ierr);
9222be375d4SJed Brown     ierr = PetscFree2(lx,ly);CHKERRQ(ierr);
92347c6ae99SBarry Smith   } else if (dd->dim == 1) {
9242be375d4SJed Brown     PetscInt *lx;
9252be375d4SJed Brown     ierr = PetscMalloc(dd->m*sizeof(PetscInt),&lx);CHKERRQ(ierr);
926397b6216SJed Brown     ierr = DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
927397b6216SJed Brown     ierr = DMDASetOwnershipRanges(da2,lx,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
9282be375d4SJed Brown     ierr = PetscFree(lx);CHKERRQ(ierr);
92947c6ae99SBarry Smith   }
93047c6ae99SBarry Smith   dd2 = (DM_DA*)da2->data;
93147c6ae99SBarry Smith 
9324dcab191SBarry Smith   /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
93325296bd5SBarry Smith   /* da2->ops->createinterpolation = da->ops->createinterpolation; copying this one causes trouble for DMSetVI */
93425296bd5SBarry Smith   da2->ops->creatematrix = da->ops->creatematrix;
93547c6ae99SBarry Smith   da2->ops->getcoloring  = da->ops->getcoloring;
93647c6ae99SBarry Smith   dd2->interptype        = dd->interptype;
93747c6ae99SBarry Smith 
93847c6ae99SBarry Smith   /* copy fill information if given */
93947c6ae99SBarry Smith   if (dd->dfill) {
94047c6ae99SBarry Smith     ierr = PetscMalloc((dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->dfill);CHKERRQ(ierr);
94147c6ae99SBarry Smith     ierr = PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));CHKERRQ(ierr);
94247c6ae99SBarry Smith   }
94347c6ae99SBarry Smith   if (dd->ofill) {
94447c6ae99SBarry Smith     ierr = PetscMalloc((dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->ofill);CHKERRQ(ierr);
94547c6ae99SBarry Smith     ierr = PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));CHKERRQ(ierr);
94647c6ae99SBarry Smith   }
94747c6ae99SBarry Smith   /* copy the refine information */
948397b6216SJed Brown   dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
949397b6216SJed Brown   dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
950397b6216SJed Brown   dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
95147c6ae99SBarry Smith 
95247c6ae99SBarry Smith   /* copy vector type information */
95347c6ae99SBarry Smith   ierr = PetscFree(da2->vectype);CHKERRQ(ierr);
95419fd82e9SBarry Smith   ierr = PetscStrallocpy(da->vectype,(char**)&da2->vectype);CHKERRQ(ierr);
95547c6ae99SBarry Smith 
956644e2e5bSBarry Smith   dd2->lf = dd->lf;
957644e2e5bSBarry Smith   dd2->lj = dd->lj;
958644e2e5bSBarry Smith 
9596e87535bSJed Brown   da2->leveldown = da->leveldown + 1;
9606e87535bSJed Brown   da2->levelup   = da->levelup;
961*8865f1eaSKarl Rupp 
9626e87535bSJed Brown   ierr = DMSetFromOptions(da2);CHKERRQ(ierr);
9636e87535bSJed Brown   ierr = DMSetUp(da2);CHKERRQ(ierr);
964ca266f36SBarry Smith   ierr = DMViewFromOptions(da2,"-dm_view");CHKERRQ(ierr);
9656e87535bSJed Brown 
966ddcf8b74SDave May   /* inject coordinates if they are set on the fine grid */
9676636e97aSMatthew G Knepley   if (da->coordinates) {
968ddcf8b74SDave May     DM         cdaf,cdac;
969ddcf8b74SDave May     Vec        coordsc,coordsf;
970ddcf8b74SDave May     VecScatter inject;
971ddcf8b74SDave May 
9726636e97aSMatthew G Knepley     ierr = DMGetCoordinateDM(da,&cdaf);CHKERRQ(ierr);
9736636e97aSMatthew G Knepley     ierr = DMGetCoordinates(da,&coordsf);CHKERRQ(ierr);
9746636e97aSMatthew G Knepley     ierr = DMGetCoordinateDM(da2,&cdac);CHKERRQ(ierr);
975b61d3410SDave May     /* force creation of the coordinate vector */
976b61d3410SDave May     ierr = DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);
9776636e97aSMatthew G Knepley     ierr = DMGetCoordinates(da2,&coordsc);CHKERRQ(ierr);
978ddcf8b74SDave May 
979e727c939SJed Brown     ierr = DMCreateInjection(cdac,cdaf,&inject);CHKERRQ(ierr);
980ddcf8b74SDave May     ierr = VecScatterBegin(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
981ddcf8b74SDave May     ierr = VecScatterEnd(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
982fcfd50ebSBarry Smith     ierr = VecScatterDestroy(&inject);CHKERRQ(ierr);
983ddcf8b74SDave May   }
984f98405f7SJed Brown 
985f98405f7SJed Brown   for (i=0; i<da->bs; i++) {
986f98405f7SJed Brown     const char *fieldname;
987f98405f7SJed Brown     ierr = DMDAGetFieldName(da,i,&fieldname);CHKERRQ(ierr);
988f98405f7SJed Brown     ierr = DMDASetFieldName(da2,i,fieldname);CHKERRQ(ierr);
989f98405f7SJed Brown   }
990f98405f7SJed Brown 
99147c6ae99SBarry Smith   *daref = da2;
99247c6ae99SBarry Smith   PetscFunctionReturn(0);
99347c6ae99SBarry Smith }
99447c6ae99SBarry Smith 
99547c6ae99SBarry Smith #undef __FUNCT__
99694013140SBarry Smith #define __FUNCT__ "DMRefineHierarchy_DA"
9977087cfbeSBarry Smith PetscErrorCode  DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
99847c6ae99SBarry Smith {
99947c6ae99SBarry Smith   PetscErrorCode ierr;
100047c6ae99SBarry Smith   PetscInt       i,n,*refx,*refy,*refz;
100147c6ae99SBarry Smith 
100247c6ae99SBarry Smith   PetscFunctionBegin;
100347c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
100447c6ae99SBarry Smith   if (nlevels < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
100547c6ae99SBarry Smith   if (nlevels == 0) PetscFunctionReturn(0);
100647c6ae99SBarry Smith   PetscValidPointer(daf,3);
100747c6ae99SBarry Smith 
1008aa219208SBarry Smith   /* Get refinement factors, defaults taken from the coarse DMDA */
100947c6ae99SBarry Smith   ierr = PetscMalloc3(nlevels,PetscInt,&refx,nlevels,PetscInt,&refy,nlevels,PetscInt,&refz);CHKERRQ(ierr);
101047c6ae99SBarry Smith   for (i=0; i<nlevels; i++) {
1011aa219208SBarry Smith     ierr = DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);CHKERRQ(ierr);
101247c6ae99SBarry Smith   }
101347c6ae99SBarry Smith   n    = nlevels;
101447c6ae99SBarry Smith   ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,PETSC_NULL);CHKERRQ(ierr);
101547c6ae99SBarry Smith   n    = nlevels;
101647c6ae99SBarry Smith   ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,PETSC_NULL);CHKERRQ(ierr);
101747c6ae99SBarry Smith   n    = nlevels;
101847c6ae99SBarry Smith   ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,PETSC_NULL);CHKERRQ(ierr);
101947c6ae99SBarry Smith 
1020aa219208SBarry Smith   ierr = DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);CHKERRQ(ierr);
102194013140SBarry Smith   ierr = DMRefine(da,((PetscObject)da)->comm,&daf[0]);CHKERRQ(ierr);
102247c6ae99SBarry Smith   for (i=1; i<nlevels; i++) {
1023aa219208SBarry Smith     ierr = DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);CHKERRQ(ierr);
102494013140SBarry Smith     ierr = DMRefine(daf[i-1],((PetscObject)da)->comm,&daf[i]);CHKERRQ(ierr);
102547c6ae99SBarry Smith   }
102647c6ae99SBarry Smith   ierr = PetscFree3(refx,refy,refz);CHKERRQ(ierr);
102747c6ae99SBarry Smith   PetscFunctionReturn(0);
102847c6ae99SBarry Smith }
102947c6ae99SBarry Smith 
103047c6ae99SBarry Smith #undef __FUNCT__
103194013140SBarry Smith #define __FUNCT__ "DMCoarsenHierarchy_DA"
10327087cfbeSBarry Smith PetscErrorCode  DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
103347c6ae99SBarry Smith {
103447c6ae99SBarry Smith   PetscErrorCode ierr;
103547c6ae99SBarry Smith   PetscInt       i;
103647c6ae99SBarry Smith 
103747c6ae99SBarry Smith   PetscFunctionBegin;
103847c6ae99SBarry Smith   PetscValidHeaderSpecific(da,DM_CLASSID,1);
103947c6ae99SBarry Smith   if (nlevels < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
104047c6ae99SBarry Smith   if (nlevels == 0) PetscFunctionReturn(0);
104147c6ae99SBarry Smith   PetscValidPointer(dac,3);
104294013140SBarry Smith   ierr = DMCoarsen(da,((PetscObject)da)->comm,&dac[0]);CHKERRQ(ierr);
104347c6ae99SBarry Smith   for (i=1; i<nlevels; i++) {
104494013140SBarry Smith     ierr = DMCoarsen(dac[i-1],((PetscObject)da)->comm,&dac[i]);CHKERRQ(ierr);
104547c6ae99SBarry Smith   }
104647c6ae99SBarry Smith   PetscFunctionReturn(0);
104747c6ae99SBarry Smith }
1048