Actual source code: dainterp.c
petsc-3.9.4 2018-09-11
2: /*
3: Code for interpolating between grids represented by DMDAs
4: */
6: /*
7: For linear elements there are two branches of code to compute the interpolation. They should compute the same results but may not. The "new version" does
8: not work for periodic domains, the old does. Change NEWVERSION to 1 to compile in the new version. Eventually when we are sure the two produce identical results
9: we will remove/merge the new version. Based on current tests, these both produce the same results. We are leaving NEWVERSION for now in the code since some
10: consider it cleaner, but old version is turned on since it handles periodic case.
11: */
12: #define NEWVERSION 0
14: #include <petsc/private/dmdaimpl.h>
16: /*@
17: DMCreateInterpolationScale - Forms L = R*1/diag(R*1) - L.*v is like a coarse grid average of the
18: nearby fine grid points.
20: Input Parameters:
21: + dac - DM that defines a coarse mesh
22: . daf - DM that defines a fine mesh
23: - mat - the restriction (or interpolation operator) from fine to coarse
25: Output Parameter:
26: . scale - the scaled vector
28: Level: developer
30: .seealso: DMCreateInterpolation()
32: @*/
33: PetscErrorCode DMCreateInterpolationScale(DM dac,DM daf,Mat mat,Vec *scale)
34: {
36: Vec fine;
37: PetscScalar one = 1.0;
40: DMCreateGlobalVector(daf,&fine);
41: DMCreateGlobalVector(dac,scale);
42: VecSet(fine,one);
43: MatRestrict(mat,fine,*scale);
44: VecDestroy(&fine);
45: VecReciprocal(*scale);
46: return(0);
47: }
49: PetscErrorCode DMCreateInterpolation_DA_1D_Q1(DM dac,DM daf,Mat *A)
50: {
51: PetscErrorCode ierr;
52: PetscInt i,i_start,m_f,Mx;
53: const PetscInt *idx_f,*idx_c;
54: PetscInt m_ghost,m_ghost_c;
55: PetscInt row,col,i_start_ghost,mx,m_c,nc,ratio;
56: PetscInt i_c,i_start_c,i_start_ghost_c,cols[2],dof;
57: PetscScalar v[2],x;
58: Mat mat;
59: DMBoundaryType bx;
60: ISLocalToGlobalMapping ltog_f,ltog_c;
64: DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);
65: DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);
66: if (bx == DM_BOUNDARY_PERIODIC) {
67: ratio = mx/Mx;
68: if (ratio*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
69: } else {
70: ratio = (mx-1)/(Mx-1);
71: if (ratio*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
72: }
74: DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);
75: DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
76: DMGetLocalToGlobalMapping(daf,<og_f);
77: ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);
79: DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
80: DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
81: DMGetLocalToGlobalMapping(dac,<og_c);
82: ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);
84: /* create interpolation matrix */
85: MatCreate(PetscObjectComm((PetscObject)dac),&mat);
86: MatSetSizes(mat,m_f,m_c,mx,Mx);
87: MatSetType(mat,MATAIJ);
88: MatSeqAIJSetPreallocation(mat,2,NULL);
89: MatMPIAIJSetPreallocation(mat,2,NULL,1,NULL);
91: /* loop over local fine grid nodes setting interpolation for those*/
92: if (!NEWVERSION) {
94: for (i=i_start; i<i_start+m_f; i++) {
95: /* convert to local "natural" numbering and then to PETSc global numbering */
96: row = idx_f[i-i_start_ghost];
98: i_c = (i/ratio); /* coarse grid node to left of fine grid node */
99: if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
100: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
102: /*
103: Only include those interpolation points that are truly
104: nonzero. Note this is very important for final grid lines
105: in x direction; since they have no right neighbor
106: */
107: x = ((PetscReal)(i - i_c*ratio))/((PetscReal)ratio);
108: nc = 0;
109: /* one left and below; or we are right on it */
110: col = (i_c-i_start_ghost_c);
111: cols[nc] = idx_c[col];
112: v[nc++] = -x + 1.0;
113: /* one right? */
114: if (i_c*ratio != i) {
115: cols[nc] = idx_c[col+1];
116: v[nc++] = x;
117: }
118: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
119: }
121: } else {
122: PetscScalar *xi;
123: PetscInt li,nxi,n;
124: PetscScalar Ni[2];
126: /* compute local coordinate arrays */
127: nxi = ratio + 1;
128: PetscMalloc1(nxi,&xi);
129: for (li=0; li<nxi; li++) {
130: xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
131: }
133: for (i=i_start; i<i_start+m_f; i++) {
134: /* convert to local "natural" numbering and then to PETSc global numbering */
135: row = idx_f[(i-i_start_ghost)];
137: i_c = (i/ratio); /* coarse grid node to left of fine grid node */
138: if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
139: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
141: /* remainders */
142: li = i - ratio * (i/ratio);
143: if (i==mx-1) li = nxi-1;
145: /* corners */
146: col = (i_c-i_start_ghost_c);
147: cols[0] = idx_c[col];
148: Ni[0] = 1.0;
149: if ((li==0) || (li==nxi-1)) {
150: MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
151: continue;
152: }
154: /* edges + interior */
155: /* remainders */
156: if (i==mx-1) i_c--;
158: col = (i_c-i_start_ghost_c);
159: cols[0] = idx_c[col]; /* one left and below; or we are right on it */
160: cols[1] = idx_c[col+1];
162: Ni[0] = 0.5*(1.0-xi[li]);
163: Ni[1] = 0.5*(1.0+xi[li]);
164: for (n=0; n<2; n++) {
165: if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1;
166: }
167: MatSetValues(mat,1,&row,2,cols,Ni,INSERT_VALUES);
168: }
169: PetscFree(xi);
170: }
171: ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
172: ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
173: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
174: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
175: MatCreateMAIJ(mat,dof,A);
176: MatDestroy(&mat);
177: return(0);
178: }
180: PetscErrorCode DMCreateInterpolation_DA_1D_Q0(DM dac,DM daf,Mat *A)
181: {
182: PetscErrorCode ierr;
183: PetscInt i,i_start,m_f,Mx;
184: const PetscInt *idx_f,*idx_c;
185: ISLocalToGlobalMapping ltog_f,ltog_c;
186: PetscInt m_ghost,m_ghost_c;
187: PetscInt row,col,i_start_ghost,mx,m_c,nc,ratio;
188: PetscInt i_c,i_start_c,i_start_ghost_c,cols[2],dof;
189: PetscScalar v[2],x;
190: Mat mat;
191: DMBoundaryType bx;
194: DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);
195: DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);
196: if (bx == DM_BOUNDARY_PERIODIC) {
197: if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
198: ratio = mx/Mx;
199: if (ratio*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
200: } else {
201: if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx);
202: ratio = (mx-1)/(Mx-1);
203: if (ratio*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
204: }
206: DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);
207: DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
208: DMGetLocalToGlobalMapping(daf,<og_f);
209: ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);
211: DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
212: DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
213: DMGetLocalToGlobalMapping(dac,<og_c);
214: ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);
216: /* create interpolation matrix */
217: MatCreate(PetscObjectComm((PetscObject)dac),&mat);
218: MatSetSizes(mat,m_f,m_c,mx,Mx);
219: MatSetType(mat,MATAIJ);
220: MatSeqAIJSetPreallocation(mat,2,NULL);
221: MatMPIAIJSetPreallocation(mat,2,NULL,0,NULL);
223: /* loop over local fine grid nodes setting interpolation for those*/
224: for (i=i_start; i<i_start+m_f; i++) {
225: /* convert to local "natural" numbering and then to PETSc global numbering */
226: row = idx_f[(i-i_start_ghost)];
228: i_c = (i/ratio); /* coarse grid node to left of fine grid node */
230: /*
231: Only include those interpolation points that are truly
232: nonzero. Note this is very important for final grid lines
233: in x direction; since they have no right neighbor
234: */
235: x = ((PetscReal)(i - i_c*ratio))/((PetscReal)ratio);
236: nc = 0;
237: /* one left and below; or we are right on it */
238: col = (i_c-i_start_ghost_c);
239: cols[nc] = idx_c[col];
240: v[nc++] = -x + 1.0;
241: /* one right? */
242: if (i_c*ratio != i) {
243: cols[nc] = idx_c[col+1];
244: v[nc++] = x;
245: }
246: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
247: }
248: ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
249: ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
250: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
251: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
252: MatCreateMAIJ(mat,dof,A);
253: MatDestroy(&mat);
254: PetscLogFlops(5.0*m_f);
255: return(0);
256: }
258: PetscErrorCode DMCreateInterpolation_DA_2D_Q1(DM dac,DM daf,Mat *A)
259: {
260: PetscErrorCode ierr;
261: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof;
262: const PetscInt *idx_c,*idx_f;
263: ISLocalToGlobalMapping ltog_f,ltog_c;
264: PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c,*dnz,*onz;
265: PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
266: PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c,col_shift,col_scale;
267: PetscMPIInt size_c,size_f,rank_f;
268: PetscScalar v[4],x,y;
269: Mat mat;
270: DMBoundaryType bx,by;
273: DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
274: DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
275: if (bx == DM_BOUNDARY_PERIODIC) {
276: if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
277: ratioi = mx/Mx;
278: if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
279: } else {
280: if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx);
281: ratioi = (mx-1)/(Mx-1);
282: if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
283: }
284: if (by == DM_BOUNDARY_PERIODIC) {
285: if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
286: ratioj = my/My;
287: if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My);
288: } else {
289: if (My < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be greater than 1",My);
290: ratioj = (my-1)/(My-1);
291: if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
292: }
294: DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
295: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
296: DMGetLocalToGlobalMapping(daf,<og_f);
297: ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);
299: DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
300: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
301: DMGetLocalToGlobalMapping(dac,<og_c);
302: ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);
304: /*
305: Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
306: The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
307: processors). It's effective length is hence 4 times its normal length, this is
308: why the col_scale is multiplied by the interpolation matrix column sizes.
309: sol_shift allows each set of 1/4 processors do its own interpolation using ITS
310: copy of the coarse vector. A bit of a hack but you do better.
312: In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
313: */
314: MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);
315: MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);
316: MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);
317: col_scale = size_f/size_c;
318: col_shift = Mx*My*(rank_f/size_c);
320: MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f,col_scale*m_c*n_c,dnz,onz);
321: for (j=j_start; j<j_start+n_f; j++) {
322: for (i=i_start; i<i_start+m_f; i++) {
323: /* convert to local "natural" numbering and then to PETSc global numbering */
324: row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
326: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
327: j_c = (j/ratioj); /* coarse grid node below fine grid node */
329: if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
330: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
331: if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
332: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
334: /*
335: Only include those interpolation points that are truly
336: nonzero. Note this is very important for final grid lines
337: in x and y directions; since they have no right/top neighbors
338: */
339: nc = 0;
340: /* one left and below; or we are right on it */
341: col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
342: cols[nc++] = col_shift + idx_c[col];
343: /* one right and below */
344: if (i_c*ratioi != i) cols[nc++] = col_shift + idx_c[col+1];
345: /* one left and above */
346: if (j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+m_ghost_c];
347: /* one right and above */
348: if (i_c*ratioi != i && j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+(m_ghost_c+1)];
349: MatPreallocateSet(row,nc,cols,dnz,onz);
350: }
351: }
352: MatCreate(PetscObjectComm((PetscObject)daf),&mat);
353: MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);
354: MatSetType(mat,MATAIJ);
355: MatSeqAIJSetPreallocation(mat,0,dnz);
356: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
357: MatPreallocateFinalize(dnz,onz);
359: /* loop over local fine grid nodes setting interpolation for those*/
360: if (!NEWVERSION) {
362: for (j=j_start; j<j_start+n_f; j++) {
363: for (i=i_start; i<i_start+m_f; i++) {
364: /* convert to local "natural" numbering and then to PETSc global numbering */
365: row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
367: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
368: j_c = (j/ratioj); /* coarse grid node below fine grid node */
370: /*
371: Only include those interpolation points that are truly
372: nonzero. Note this is very important for final grid lines
373: in x and y directions; since they have no right/top neighbors
374: */
375: x = ((PetscReal)(i - i_c*ratioi))/((PetscReal)ratioi);
376: y = ((PetscReal)(j - j_c*ratioj))/((PetscReal)ratioj);
378: nc = 0;
379: /* one left and below; or we are right on it */
380: col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
381: cols[nc] = col_shift + idx_c[col];
382: v[nc++] = x*y - x - y + 1.0;
383: /* one right and below */
384: if (i_c*ratioi != i) {
385: cols[nc] = col_shift + idx_c[col+1];
386: v[nc++] = -x*y + x;
387: }
388: /* one left and above */
389: if (j_c*ratioj != j) {
390: cols[nc] = col_shift + idx_c[col+m_ghost_c];
391: v[nc++] = -x*y + y;
392: }
393: /* one right and above */
394: if (j_c*ratioj != j && i_c*ratioi != i) {
395: cols[nc] = col_shift + idx_c[col+(m_ghost_c+1)];
396: v[nc++] = x*y;
397: }
398: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
399: }
400: }
402: } else {
403: PetscScalar Ni[4];
404: PetscScalar *xi,*eta;
405: PetscInt li,nxi,lj,neta;
407: /* compute local coordinate arrays */
408: nxi = ratioi + 1;
409: neta = ratioj + 1;
410: PetscMalloc1(nxi,&xi);
411: PetscMalloc1(neta,&eta);
412: for (li=0; li<nxi; li++) {
413: xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
414: }
415: for (lj=0; lj<neta; lj++) {
416: eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
417: }
419: /* loop over local fine grid nodes setting interpolation for those*/
420: for (j=j_start; j<j_start+n_f; j++) {
421: for (i=i_start; i<i_start+m_f; i++) {
422: /* convert to local "natural" numbering and then to PETSc global numbering */
423: row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
425: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
426: j_c = (j/ratioj); /* coarse grid node below fine grid node */
428: /* remainders */
429: li = i - ratioi * (i/ratioi);
430: if (i==mx-1) li = nxi-1;
431: lj = j - ratioj * (j/ratioj);
432: if (j==my-1) lj = neta-1;
434: /* corners */
435: col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
436: cols[0] = col_shift + idx_c[col]; /* left, below */
437: Ni[0] = 1.0;
438: if ((li==0) || (li==nxi-1)) {
439: if ((lj==0) || (lj==neta-1)) {
440: MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
441: continue;
442: }
443: }
445: /* edges + interior */
446: /* remainders */
447: if (i==mx-1) i_c--;
448: if (j==my-1) j_c--;
450: col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
451: cols[0] = col_shift + idx_c[col]; /* left, below */
452: cols[1] = col_shift + idx_c[col+1]; /* right, below */
453: cols[2] = col_shift + idx_c[col+m_ghost_c]; /* left, above */
454: cols[3] = col_shift + idx_c[col+(m_ghost_c+1)]; /* right, above */
456: Ni[0] = 0.25*(1.0-xi[li])*(1.0-eta[lj]);
457: Ni[1] = 0.25*(1.0+xi[li])*(1.0-eta[lj]);
458: Ni[2] = 0.25*(1.0-xi[li])*(1.0+eta[lj]);
459: Ni[3] = 0.25*(1.0+xi[li])*(1.0+eta[lj]);
461: nc = 0;
462: if (PetscAbsScalar(Ni[0])<1.0e-32) cols[0]=-1;
463: if (PetscAbsScalar(Ni[1])<1.0e-32) cols[1]=-1;
464: if (PetscAbsScalar(Ni[2])<1.0e-32) cols[2]=-1;
465: if (PetscAbsScalar(Ni[3])<1.0e-32) cols[3]=-1;
467: MatSetValues(mat,1,&row,4,cols,Ni,INSERT_VALUES);
468: }
469: }
470: PetscFree(xi);
471: PetscFree(eta);
472: }
473: ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
474: ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
475: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
476: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
477: MatCreateMAIJ(mat,dof,A);
478: MatDestroy(&mat);
479: return(0);
480: }
482: /*
483: Contributed by Andrei Draganescu <aidraga@sandia.gov>
484: */
485: PetscErrorCode DMCreateInterpolation_DA_2D_Q0(DM dac,DM daf,Mat *A)
486: {
487: PetscErrorCode ierr;
488: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof;
489: const PetscInt *idx_c,*idx_f;
490: ISLocalToGlobalMapping ltog_f,ltog_c;
491: PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c,*dnz,*onz;
492: PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
493: PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c,col_shift,col_scale;
494: PetscMPIInt size_c,size_f,rank_f;
495: PetscScalar v[4];
496: Mat mat;
497: DMBoundaryType bx,by;
500: DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
501: DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
502: if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
503: if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
504: ratioi = mx/Mx;
505: ratioj = my/My;
506: if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x");
507: if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y");
508: if (ratioi != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 2");
509: if (ratioj != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 2");
511: DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
512: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
513: DMGetLocalToGlobalMapping(daf,<og_f);
514: ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);
516: DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
517: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
518: DMGetLocalToGlobalMapping(dac,<og_c);
519: ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);
521: /*
522: Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
523: The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
524: processors). It's effective length is hence 4 times its normal length, this is
525: why the col_scale is multiplied by the interpolation matrix column sizes.
526: sol_shift allows each set of 1/4 processors do its own interpolation using ITS
527: copy of the coarse vector. A bit of a hack but you do better.
529: In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
530: */
531: MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);
532: MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);
533: MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);
534: col_scale = size_f/size_c;
535: col_shift = Mx*My*(rank_f/size_c);
537: MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f,col_scale*m_c*n_c,dnz,onz);
538: for (j=j_start; j<j_start+n_f; j++) {
539: for (i=i_start; i<i_start+m_f; i++) {
540: /* convert to local "natural" numbering and then to PETSc global numbering */
541: row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
543: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
544: j_c = (j/ratioj); /* coarse grid node below fine grid node */
546: if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
547: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
548: if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
549: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
551: /*
552: Only include those interpolation points that are truly
553: nonzero. Note this is very important for final grid lines
554: in x and y directions; since they have no right/top neighbors
555: */
556: nc = 0;
557: /* one left and below; or we are right on it */
558: col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
559: cols[nc++] = col_shift + idx_c[col];
560: MatPreallocateSet(row,nc,cols,dnz,onz);
561: }
562: }
563: MatCreate(PetscObjectComm((PetscObject)daf),&mat);
564: MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);
565: MatSetType(mat,MATAIJ);
566: MatSeqAIJSetPreallocation(mat,0,dnz);
567: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
568: MatPreallocateFinalize(dnz,onz);
570: /* loop over local fine grid nodes setting interpolation for those*/
571: for (j=j_start; j<j_start+n_f; j++) {
572: for (i=i_start; i<i_start+m_f; i++) {
573: /* convert to local "natural" numbering and then to PETSc global numbering */
574: row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
576: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
577: j_c = (j/ratioj); /* coarse grid node below fine grid node */
578: nc = 0;
579: /* one left and below; or we are right on it */
580: col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
581: cols[nc] = col_shift + idx_c[col];
582: v[nc++] = 1.0;
584: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
585: }
586: }
587: ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
588: ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
589: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
590: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
591: MatCreateMAIJ(mat,dof,A);
592: MatDestroy(&mat);
593: PetscLogFlops(13.0*m_f*n_f);
594: return(0);
595: }
597: /*
598: Contributed by Jianming Yang <jianming-yang@uiowa.edu>
599: */
600: PetscErrorCode DMCreateInterpolation_DA_3D_Q0(DM dac,DM daf,Mat *A)
601: {
602: PetscErrorCode ierr;
603: PetscInt i,j,l,i_start,j_start,l_start,m_f,n_f,p_f,Mx,My,Mz,dof;
604: const PetscInt *idx_c,*idx_f;
605: ISLocalToGlobalMapping ltog_f,ltog_c;
606: PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c,nc,*dnz,*onz;
607: PetscInt row,col,i_start_ghost,j_start_ghost,l_start_ghost,cols[8],mx,m_c,my,n_c,mz,p_c,ratioi,ratioj,ratiol;
608: PetscInt i_c,j_c,l_c,i_start_c,j_start_c,l_start_c,i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,col_shift,col_scale;
609: PetscMPIInt size_c,size_f,rank_f;
610: PetscScalar v[8];
611: Mat mat;
612: DMBoundaryType bx,by,bz;
615: DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
616: if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
617: if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
618: if (!Mz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be positive",Mz);
619: DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
620: ratioi = mx/Mx;
621: ratioj = my/My;
622: ratiol = mz/Mz;
623: if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x");
624: if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y");
625: if (ratiol*Mz != mz) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in z");
626: if (ratioi != 2 && ratioi != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 1 or 2");
627: if (ratioj != 2 && ratioj != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 1 or 2");
628: if (ratiol != 2 && ratiol != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in z must be 1 or 2");
630: DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
631: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
632: DMGetLocalToGlobalMapping(daf,<og_f);
633: ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);
635: DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
636: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
637: DMGetLocalToGlobalMapping(dac,<og_c);
638: ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);
640: /*
641: Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
642: The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
643: processors). It's effective length is hence 4 times its normal length, this is
644: why the col_scale is multiplied by the interpolation matrix column sizes.
645: sol_shift allows each set of 1/4 processors do its own interpolation using ITS
646: copy of the coarse vector. A bit of a hack but you do better.
648: In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
649: */
650: MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);
651: MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);
652: MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);
653: col_scale = size_f/size_c;
654: col_shift = Mx*My*Mz*(rank_f/size_c);
656: MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f*p_f,col_scale*m_c*n_c*p_c,dnz,onz);
657: for (l=l_start; l<l_start+p_f; l++) {
658: for (j=j_start; j<j_start+n_f; j++) {
659: for (i=i_start; i<i_start+m_f; i++) {
660: /* convert to local "natural" numbering and then to PETSc global numbering */
661: row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
663: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
664: j_c = (j/ratioj); /* coarse grid node below fine grid node */
665: l_c = (l/ratiol);
667: if (l_c < l_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
668: l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
669: if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
670: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
671: if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
672: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
674: /*
675: Only include those interpolation points that are truly
676: nonzero. Note this is very important for final grid lines
677: in x and y directions; since they have no right/top neighbors
678: */
679: nc = 0;
680: /* one left and below; or we are right on it */
681: col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
682: cols[nc++] = col_shift + idx_c[col];
683: MatPreallocateSet(row,nc,cols,dnz,onz);
684: }
685: }
686: }
687: MatCreate(PetscObjectComm((PetscObject)daf),&mat);
688: MatSetSizes(mat,m_f*n_f*p_f,col_scale*m_c*n_c*p_c,mx*my*mz,col_scale*Mx*My*Mz);
689: MatSetType(mat,MATAIJ);
690: MatSeqAIJSetPreallocation(mat,0,dnz);
691: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
692: MatPreallocateFinalize(dnz,onz);
694: /* loop over local fine grid nodes setting interpolation for those*/
695: for (l=l_start; l<l_start+p_f; l++) {
696: for (j=j_start; j<j_start+n_f; j++) {
697: for (i=i_start; i<i_start+m_f; i++) {
698: /* convert to local "natural" numbering and then to PETSc global numbering */
699: row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
701: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
702: j_c = (j/ratioj); /* coarse grid node below fine grid node */
703: l_c = (l/ratiol);
704: nc = 0;
705: /* one left and below; or we are right on it */
706: col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
707: cols[nc] = col_shift + idx_c[col];
708: v[nc++] = 1.0;
710: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
711: }
712: }
713: }
714: ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
715: ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
716: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
717: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
718: MatCreateMAIJ(mat,dof,A);
719: MatDestroy(&mat);
720: PetscLogFlops(13.0*m_f*n_f*p_f);
721: return(0);
722: }
724: PetscErrorCode DMCreateInterpolation_DA_3D_Q1(DM dac,DM daf,Mat *A)
725: {
726: PetscErrorCode ierr;
727: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof,l;
728: const PetscInt *idx_c,*idx_f;
729: ISLocalToGlobalMapping ltog_f,ltog_c;
730: PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c,Mz,mz;
731: PetscInt row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok;
732: PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
733: PetscInt l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c;
734: PetscInt l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz;
735: PetscScalar v[8],x,y,z;
736: Mat mat;
737: DMBoundaryType bx,by,bz;
740: DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
741: DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
742: if (mx == Mx) {
743: ratioi = 1;
744: } else if (bx == DM_BOUNDARY_PERIODIC) {
745: if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
746: ratioi = mx/Mx;
747: if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
748: } else {
749: if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx);
750: ratioi = (mx-1)/(Mx-1);
751: if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
752: }
753: if (my == My) {
754: ratioj = 1;
755: } else if (by == DM_BOUNDARY_PERIODIC) {
756: if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
757: ratioj = my/My;
758: if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My);
759: } else {
760: if (My < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be greater than 1",My);
761: ratioj = (my-1)/(My-1);
762: if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
763: }
764: if (mz == Mz) {
765: ratiok = 1;
766: } else if (bz == DM_BOUNDARY_PERIODIC) {
767: if (!Mz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be positive",Mz);
768: ratiok = mz/Mz;
769: if (ratiok*Mz != mz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mz/Mz must be integer: mz %D Mz %D",mz,Mz);
770: } else {
771: if (Mz < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be greater than 1",Mz);
772: ratiok = (mz-1)/(Mz-1);
773: if (ratiok*(Mz-1) != mz-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %D Mz %D",mz,Mz);
774: }
776: DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
777: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
778: DMGetLocalToGlobalMapping(daf,<og_f);
779: ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);
781: DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
782: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
783: DMGetLocalToGlobalMapping(dac,<og_c);
784: ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);
786: /* create interpolation matrix, determining exact preallocation */
787: MatPreallocateInitialize(PetscObjectComm((PetscObject)dac),m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);
788: /* loop over local fine grid nodes counting interpolating points */
789: for (l=l_start; l<l_start+p_f; l++) {
790: for (j=j_start; j<j_start+n_f; j++) {
791: for (i=i_start; i<i_start+m_f; i++) {
792: /* convert to local "natural" numbering and then to PETSc global numbering */
793: row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
794: i_c = (i/ratioi);
795: j_c = (j/ratioj);
796: l_c = (l/ratiok);
797: if (l_c < l_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
798: l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
799: if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
800: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
801: if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
802: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
804: /*
805: Only include those interpolation points that are truly
806: nonzero. Note this is very important for final grid lines
807: in x and y directions; since they have no right/top neighbors
808: */
809: nc = 0;
810: col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
811: cols[nc++] = idx_c[col];
812: if (i_c*ratioi != i) {
813: cols[nc++] = idx_c[col+1];
814: }
815: if (j_c*ratioj != j) {
816: cols[nc++] = idx_c[col+m_ghost_c];
817: }
818: if (l_c*ratiok != l) {
819: cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c];
820: }
821: if (j_c*ratioj != j && i_c*ratioi != i) {
822: cols[nc++] = idx_c[col+(m_ghost_c+1)];
823: }
824: if (j_c*ratioj != j && l_c*ratiok != l) {
825: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)];
826: }
827: if (i_c*ratioi != i && l_c*ratiok != l) {
828: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)];
829: }
830: if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
831: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)];
832: }
833: MatPreallocateSet(row,nc,cols,dnz,onz);
834: }
835: }
836: }
837: MatCreate(PetscObjectComm((PetscObject)dac),&mat);
838: MatSetSizes(mat,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz);
839: MatSetType(mat,MATAIJ);
840: MatSeqAIJSetPreallocation(mat,0,dnz);
841: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
842: MatPreallocateFinalize(dnz,onz);
844: /* loop over local fine grid nodes setting interpolation for those*/
845: if (!NEWVERSION) {
847: for (l=l_start; l<l_start+p_f; l++) {
848: for (j=j_start; j<j_start+n_f; j++) {
849: for (i=i_start; i<i_start+m_f; i++) {
850: /* convert to local "natural" numbering and then to PETSc global numbering */
851: row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
853: i_c = (i/ratioi);
854: j_c = (j/ratioj);
855: l_c = (l/ratiok);
857: /*
858: Only include those interpolation points that are truly
859: nonzero. Note this is very important for final grid lines
860: in x and y directions; since they have no right/top neighbors
861: */
862: x = ((PetscReal)(i - i_c*ratioi))/((PetscReal)ratioi);
863: y = ((PetscReal)(j - j_c*ratioj))/((PetscReal)ratioj);
864: z = ((PetscReal)(l - l_c*ratiok))/((PetscReal)ratiok);
866: nc = 0;
867: /* one left and below; or we are right on it */
868: col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));
870: cols[nc] = idx_c[col];
871: v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
873: if (i_c*ratioi != i) {
874: cols[nc] = idx_c[col+1];
875: v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
876: }
878: if (j_c*ratioj != j) {
879: cols[nc] = idx_c[col+m_ghost_c];
880: v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
881: }
883: if (l_c*ratiok != l) {
884: cols[nc] = idx_c[col+m_ghost_c*n_ghost_c];
885: v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
886: }
888: if (j_c*ratioj != j && i_c*ratioi != i) {
889: cols[nc] = idx_c[col+(m_ghost_c+1)];
890: v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
891: }
893: if (j_c*ratioj != j && l_c*ratiok != l) {
894: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)];
895: v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
896: }
898: if (i_c*ratioi != i && l_c*ratiok != l) {
899: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)];
900: v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
901: }
903: if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
904: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)];
905: v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
906: }
907: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
908: }
909: }
910: }
912: } else {
913: PetscScalar *xi,*eta,*zeta;
914: PetscInt li,nxi,lj,neta,lk,nzeta,n;
915: PetscScalar Ni[8];
917: /* compute local coordinate arrays */
918: nxi = ratioi + 1;
919: neta = ratioj + 1;
920: nzeta = ratiok + 1;
921: PetscMalloc1(nxi,&xi);
922: PetscMalloc1(neta,&eta);
923: PetscMalloc1(nzeta,&zeta);
924: for (li=0; li<nxi; li++) xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
925: for (lj=0; lj<neta; lj++) eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
926: for (lk=0; lk<nzeta; lk++) zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));
928: for (l=l_start; l<l_start+p_f; l++) {
929: for (j=j_start; j<j_start+n_f; j++) {
930: for (i=i_start; i<i_start+m_f; i++) {
931: /* convert to local "natural" numbering and then to PETSc global numbering */
932: row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
934: i_c = (i/ratioi);
935: j_c = (j/ratioj);
936: l_c = (l/ratiok);
938: /* remainders */
939: li = i - ratioi * (i/ratioi);
940: if (i==mx-1) li = nxi-1;
941: lj = j - ratioj * (j/ratioj);
942: if (j==my-1) lj = neta-1;
943: lk = l - ratiok * (l/ratiok);
944: if (l==mz-1) lk = nzeta-1;
946: /* corners */
947: col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));
948: cols[0] = idx_c[col];
949: Ni[0] = 1.0;
950: if ((li==0) || (li==nxi-1)) {
951: if ((lj==0) || (lj==neta-1)) {
952: if ((lk==0) || (lk==nzeta-1)) {
953: MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
954: continue;
955: }
956: }
957: }
959: /* edges + interior */
960: /* remainders */
961: if (i==mx-1) i_c--;
962: if (j==my-1) j_c--;
963: if (l==mz-1) l_c--;
965: col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
966: cols[0] = idx_c[col]; /* one left and below; or we are right on it */
967: cols[1] = idx_c[col+1]; /* one right and below */
968: cols[2] = idx_c[col+m_ghost_c]; /* one left and above */
969: cols[3] = idx_c[col+(m_ghost_c+1)]; /* one right and above */
971: cols[4] = idx_c[col+m_ghost_c*n_ghost_c]; /* one left and below and front; or we are right on it */
972: cols[5] = idx_c[col+(m_ghost_c*n_ghost_c+1)]; /* one right and below, and front */
973: cols[6] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)]; /* one left and above and front*/
974: cols[7] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)]; /* one right and above and front */
976: Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
977: Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
978: Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
979: Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
981: Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
982: Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
983: Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
984: Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
986: for (n=0; n<8; n++) {
987: if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1;
988: }
989: MatSetValues(mat,1,&row,8,cols,Ni,INSERT_VALUES);
991: }
992: }
993: }
994: PetscFree(xi);
995: PetscFree(eta);
996: PetscFree(zeta);
997: }
998: ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
999: ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
1000: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1001: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
1003: MatCreateMAIJ(mat,dof,A);
1004: MatDestroy(&mat);
1005: return(0);
1006: }
1008: PetscErrorCode DMCreateInterpolation_DA(DM dac,DM daf,Mat *A,Vec *scale)
1009: {
1010: PetscErrorCode ierr;
1011: PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1012: DMBoundaryType bxc,byc,bzc,bxf,byf,bzf;
1013: DMDAStencilType stc,stf;
1014: DM_DA *ddc = (DM_DA*)dac->data;
1022: DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1023: DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1024: if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1025: if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1026: if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1027: if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1028: if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1029: if (Mc < 2 && Mf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
1030: if (dimc > 1 && Nc < 2 && Nf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction");
1031: if (dimc > 2 && Pc < 2 && Pf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction");
1033: if (ddc->interptype == DMDA_Q1) {
1034: if (dimc == 1) {
1035: DMCreateInterpolation_DA_1D_Q1(dac,daf,A);
1036: } else if (dimc == 2) {
1037: DMCreateInterpolation_DA_2D_Q1(dac,daf,A);
1038: } else if (dimc == 3) {
1039: DMCreateInterpolation_DA_3D_Q1(dac,daf,A);
1040: } else SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype);
1041: } else if (ddc->interptype == DMDA_Q0) {
1042: if (dimc == 1) {
1043: DMCreateInterpolation_DA_1D_Q0(dac,daf,A);
1044: } else if (dimc == 2) {
1045: DMCreateInterpolation_DA_2D_Q0(dac,daf,A);
1046: } else if (dimc == 3) {
1047: DMCreateInterpolation_DA_3D_Q0(dac,daf,A);
1048: } else SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype);
1049: }
1050: if (scale) {
1051: DMCreateInterpolationScale((DM)dac,(DM)daf,*A,scale);
1052: }
1053: return(0);
1054: }
1056: PetscErrorCode DMCreateInjection_DA_1D(DM dac,DM daf,VecScatter *inject)
1057: {
1058: PetscErrorCode ierr;
1059: PetscInt i,i_start,m_f,Mx,dof;
1060: const PetscInt *idx_f;
1061: ISLocalToGlobalMapping ltog_f;
1062: PetscInt m_ghost,m_ghost_c;
1063: PetscInt row,i_start_ghost,mx,m_c,nc,ratioi;
1064: PetscInt i_start_c,i_start_ghost_c;
1065: PetscInt *cols;
1066: DMBoundaryType bx;
1067: Vec vecf,vecc;
1068: IS isf;
1071: DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);
1072: DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);
1073: if (bx == DM_BOUNDARY_PERIODIC) {
1074: ratioi = mx/Mx;
1075: if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
1076: } else {
1077: ratioi = (mx-1)/(Mx-1);
1078: if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
1079: }
1081: DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);
1082: DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
1083: DMGetLocalToGlobalMapping(daf,<og_f);
1084: ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);
1086: DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
1087: DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
1090: /* loop over local fine grid nodes setting interpolation for those*/
1091: nc = 0;
1092: PetscMalloc1(m_f,&cols);
1095: for (i=i_start_c; i<i_start_c+m_c; i++) {
1096: PetscInt i_f = i*ratioi;
1098: if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\ni_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1100: row = idx_f[(i_f-i_start_ghost)];
1101: cols[nc++] = row;
1102: }
1104: ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
1105: ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);
1106: DMGetGlobalVector(dac,&vecc);
1107: DMGetGlobalVector(daf,&vecf);
1108: VecScatterCreate(vecf,isf,vecc,NULL,inject);
1109: DMRestoreGlobalVector(dac,&vecc);
1110: DMRestoreGlobalVector(daf,&vecf);
1111: ISDestroy(&isf);
1112: return(0);
1113: }
1115: PetscErrorCode DMCreateInjection_DA_2D(DM dac,DM daf,VecScatter *inject)
1116: {
1117: PetscErrorCode ierr;
1118: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof;
1119: const PetscInt *idx_c,*idx_f;
1120: ISLocalToGlobalMapping ltog_f,ltog_c;
1121: PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c;
1122: PetscInt row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj;
1123: PetscInt i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
1124: PetscInt *cols;
1125: DMBoundaryType bx,by;
1126: Vec vecf,vecc;
1127: IS isf;
1130: DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
1131: DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
1132: if (bx == DM_BOUNDARY_PERIODIC) {
1133: ratioi = mx/Mx;
1134: if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
1135: } else {
1136: ratioi = (mx-1)/(Mx-1);
1137: if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
1138: }
1139: if (by == DM_BOUNDARY_PERIODIC) {
1140: ratioj = my/My;
1141: if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My);
1142: } else {
1143: ratioj = (my-1)/(My-1);
1144: if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
1145: }
1147: DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
1148: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
1149: DMGetLocalToGlobalMapping(daf,<og_f);
1150: ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);
1152: DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
1153: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
1154: DMGetLocalToGlobalMapping(dac,<og_c);
1155: ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);
1157: /* loop over local fine grid nodes setting interpolation for those*/
1158: nc = 0;
1159: PetscMalloc1(n_f*m_f,&cols);
1160: for (j=j_start_c; j<j_start_c+n_c; j++) {
1161: for (i=i_start_c; i<i_start_c+m_c; i++) {
1162: PetscInt i_f = i*ratioi,j_f = j*ratioj;
1163: if (j_f < j_start_ghost || j_f >= j_start_ghost+n_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
1164: j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost);
1165: if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
1166: i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1167: row = idx_f[(m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))];
1168: cols[nc++] = row;
1169: }
1170: }
1171: ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
1172: ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
1174: ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);
1175: DMGetGlobalVector(dac,&vecc);
1176: DMGetGlobalVector(daf,&vecf);
1177: VecScatterCreate(vecf,isf,vecc,NULL,inject);
1178: DMRestoreGlobalVector(dac,&vecc);
1179: DMRestoreGlobalVector(daf,&vecf);
1180: ISDestroy(&isf);
1181: return(0);
1182: }
1184: PetscErrorCode DMCreateInjection_DA_3D(DM dac,DM daf,VecScatter *inject)
1185: {
1186: PetscErrorCode ierr;
1187: PetscInt i,j,k,i_start,j_start,k_start,m_f,n_f,p_f,Mx,My,Mz;
1188: PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c;
1189: PetscInt i_start_ghost,j_start_ghost,k_start_ghost;
1190: PetscInt mx,my,mz,ratioi,ratioj,ratiok;
1191: PetscInt i_start_c,j_start_c,k_start_c;
1192: PetscInt m_c,n_c,p_c;
1193: PetscInt i_start_ghost_c,j_start_ghost_c,k_start_ghost_c;
1194: PetscInt row,nc,dof;
1195: const PetscInt *idx_c,*idx_f;
1196: ISLocalToGlobalMapping ltog_f,ltog_c;
1197: PetscInt *cols;
1198: DMBoundaryType bx,by,bz;
1199: Vec vecf,vecc;
1200: IS isf;
1203: DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
1204: DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
1206: if (bx == DM_BOUNDARY_PERIODIC) {
1207: ratioi = mx/Mx;
1208: if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx);
1209: } else {
1210: ratioi = (mx-1)/(Mx-1);
1211: if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
1212: }
1213: if (by == DM_BOUNDARY_PERIODIC) {
1214: ratioj = my/My;
1215: if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My);
1216: } else {
1217: ratioj = (my-1)/(My-1);
1218: if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
1219: }
1220: if (bz == DM_BOUNDARY_PERIODIC) {
1221: ratiok = mz/Mz;
1222: if (ratiok*Mz != mz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mz/Mz must be integer: mz %D My %D",mz,Mz);
1223: } else {
1224: ratiok = (mz-1)/(Mz-1);
1225: if (ratiok*(Mz-1) != mz-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %D Mz %D",mz,Mz);
1226: }
1228: DMDAGetCorners(daf,&i_start,&j_start,&k_start,&m_f,&n_f,&p_f);
1229: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&k_start_ghost,&m_ghost,&n_ghost,&p_ghost);
1230: DMGetLocalToGlobalMapping(daf,<og_f);
1231: ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);
1233: DMDAGetCorners(dac,&i_start_c,&j_start_c,&k_start_c,&m_c,&n_c,&p_c);
1234: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&k_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
1235: DMGetLocalToGlobalMapping(dac,<og_c);
1236: ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);
1239: /* loop over local fine grid nodes setting interpolation for those*/
1240: nc = 0;
1241: PetscMalloc1(n_f*m_f*p_f,&cols);
1242: for (k=k_start_c; k<k_start_c+p_c; k++) {
1243: for (j=j_start_c; j<j_start_c+n_c; j++) {
1244: for (i=i_start_c; i<i_start_c+m_c; i++) {
1245: PetscInt i_f = i*ratioi,j_f = j*ratioj,k_f = k*ratiok;
1246: if (k_f < k_start_ghost || k_f >= k_start_ghost+p_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA "
1247: "k_c %D k_f %D fine ghost range [%D,%D]",k,k_f,k_start_ghost,k_start_ghost+p_ghost);
1248: if (j_f < j_start_ghost || j_f >= j_start_ghost+n_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA "
1249: "j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost);
1250: if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA "
1251: "i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1252: row = idx_f[(m_ghost*n_ghost*(k_f-k_start_ghost) + m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))];
1253: cols[nc++] = row;
1254: }
1255: }
1256: }
1257: ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
1258: ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
1260: ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);
1261: DMGetGlobalVector(dac,&vecc);
1262: DMGetGlobalVector(daf,&vecf);
1263: VecScatterCreate(vecf,isf,vecc,NULL,inject);
1264: DMRestoreGlobalVector(dac,&vecc);
1265: DMRestoreGlobalVector(daf,&vecf);
1266: ISDestroy(&isf);
1267: return(0);
1268: }
1270: PetscErrorCode DMCreateInjection_DA(DM dac,DM daf,Mat *mat)
1271: {
1272: PetscErrorCode ierr;
1273: PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1274: DMBoundaryType bxc,byc,bzc,bxf,byf,bzf;
1275: DMDAStencilType stc,stf;
1276: VecScatter inject = NULL;
1283: DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1284: DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1285: if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1286: if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1287: if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1288: if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1289: if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1290: if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
1291: if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction");
1292: if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction");
1294: if (dimc == 1) {
1295: DMCreateInjection_DA_1D(dac,daf,&inject);
1296: } else if (dimc == 2) {
1297: DMCreateInjection_DA_2D(dac,daf,&inject);
1298: } else if (dimc == 3) {
1299: DMCreateInjection_DA_3D(dac,daf,&inject);
1300: }
1301: MatCreateScatter(PetscObjectComm((PetscObject)inject), inject, mat);
1302: VecScatterDestroy(&inject);
1303: return(0);
1304: }
1306: PetscErrorCode DMCreateAggregates_DA(DM dac,DM daf,Mat *rest)
1307: {
1308: PetscErrorCode ierr;
1309: PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc;
1310: PetscInt dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1311: DMBoundaryType bxc,byc,bzc,bxf,byf,bzf;
1312: DMDAStencilType stc,stf;
1313: PetscInt i,j,l;
1314: PetscInt i_start,j_start,l_start, m_f,n_f,p_f;
1315: PetscInt i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost;
1316: const PetscInt *idx_f;
1317: PetscInt i_c,j_c,l_c;
1318: PetscInt i_start_c,j_start_c,l_start_c, m_c,n_c,p_c;
1319: PetscInt i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c;
1320: const PetscInt *idx_c;
1321: PetscInt d;
1322: PetscInt a;
1323: PetscInt max_agg_size;
1324: PetscInt *fine_nodes;
1325: PetscScalar *one_vec;
1326: PetscInt fn_idx;
1327: ISLocalToGlobalMapping ltogmf,ltogmc;
1334: DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1335: DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1336: if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1337: if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1338: if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1339: if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1340: if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1342: if (Mf < Mc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Mc %D, Mf %D", Mc, Mf);
1343: if (Nf < Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Nc %D, Nf %D", Nc, Nf);
1344: if (Pf < Pc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Pc %D, Pf %D", Pc, Pf);
1346: if (Pc < 0) Pc = 1;
1347: if (Pf < 0) Pf = 1;
1348: if (Nc < 0) Nc = 1;
1349: if (Nf < 0) Nf = 1;
1351: DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
1352: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
1354: DMGetLocalToGlobalMapping(daf,<ogmf);
1355: ISLocalToGlobalMappingGetIndices(ltogmf,&idx_f);
1357: DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
1358: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
1360: DMGetLocalToGlobalMapping(dac,<ogmc);
1361: ISLocalToGlobalMappingGetIndices(ltogmc,&idx_c);
1363: /*
1364: Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios
1365: for dimension 1 and 2 respectively.
1366: Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1)
1367: and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate.
1368: Each specific dof on the fine grid is mapped to one dof on the coarse grid.
1369: */
1371: max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1);
1373: /* create the matrix that will contain the restriction operator */
1374: MatCreateAIJ(PetscObjectComm((PetscObject)daf), m_c*n_c*p_c*dofc, m_f*n_f*p_f*doff, Mc*Nc*Pc*dofc, Mf*Nf*Pf*doff,
1375: max_agg_size, NULL, max_agg_size, NULL, rest);
1377: /* store nodes in the fine grid here */
1378: PetscMalloc2(max_agg_size, &one_vec,max_agg_size, &fine_nodes);
1379: for (i=0; i<max_agg_size; i++) one_vec[i] = 1.0;
1381: /* loop over all coarse nodes */
1382: for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) {
1383: for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) {
1384: for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) {
1385: for (d=0; d<dofc; d++) {
1386: /* convert to local "natural" numbering and then to PETSc global numbering */
1387: a = idx_c[dofc*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c))] + d;
1389: fn_idx = 0;
1390: /* Corresponding fine points are all points (i_f, j_f, l_f) such that
1391: i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc
1392: (same for other dimensions)
1393: */
1394: for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) {
1395: for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) {
1396: for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) {
1397: fine_nodes[fn_idx] = idx_f[doff*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))] + d;
1398: fn_idx++;
1399: }
1400: }
1401: }
1402: /* add all these points to one aggregate */
1403: MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);
1404: }
1405: }
1406: }
1407: }
1408: ISLocalToGlobalMappingRestoreIndices(ltogmf,&idx_f);
1409: ISLocalToGlobalMappingRestoreIndices(ltogmc,&idx_c);
1410: PetscFree2(one_vec,fine_nodes);
1411: MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);
1412: MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);
1413: return(0);
1414: }