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/daimpl.h> /*I "petscdmda.h" I*/
18: /*@
19: DMCreateInterpolationScale - Forms L = R*1/diag(R*1) - L.*v is like a coarse grid average of the
20: nearby fine grid points.
22: Input Parameters:
23: + dac - DM that defines a coarse mesh
24: . daf - DM that defines a fine mesh
25: - mat - the restriction (or interpolation operator) from fine to coarse
27: Output Parameter:
28: . scale - the scaled vector
30: Level: developer
32: .seealso: DMCreateInterpolation()
34: @*/
35: PetscErrorCodeDMCreateInterpolationScale(DM dac,DM daf,Mat mat,Vec *scale) 36: {
38: Vec fine;
39: PetscScalar one = 1.0;
42: DMCreateGlobalVector(daf,&fine);
43: DMCreateGlobalVector(dac,scale);
44: VecSet(fine,one);
45: MatRestrict(mat,fine,*scale);
46: VecDestroy(&fine);
47: VecReciprocal(*scale);
48: return(0);
49: }
53: PetscErrorCode DMCreateInterpolation_DA_1D_Q1(DM dac,DM daf,Mat *A) 54: {
55: PetscErrorCode ierr;
56: PetscInt i,i_start,m_f,Mx,*idx_f;
57: PetscInt m_ghost,*idx_c,m_ghost_c;
58: PetscInt row,col,i_start_ghost,mx,m_c,nc,ratio;
59: PetscInt i_c,i_start_c,i_start_ghost_c,cols[2],dof;
60: PetscScalar v[2],x;
61: Mat mat;
62: DMDABoundaryType bx;
63: 65: DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);
66: DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);
67: if (bx == DMDA_BOUNDARY_PERIODIC){
68: ratio = mx/Mx;
69: 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);
70: } else {
71: ratio = (mx-1)/(Mx-1);
72: 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);
73: }
74: 75: DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);
76: DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
77: DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
78: 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: DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
82: 83: /* create interpolation matrix */
84: MatCreate(((PetscObject)dac)->comm,&mat);
85: MatSetSizes(mat,m_f,m_c,mx,Mx);
86: MatSetType(mat,MATAIJ);
87: MatSeqAIJSetPreallocation(mat,2,PETSC_NULL);
88: MatMPIAIJSetPreallocation(mat,2,PETSC_NULL,1,PETSC_NULL);
89: 90: /* loop over local fine grid nodes setting interpolation for those*/
91: if (!NEWVERSION) {
93: for (i=i_start; i<i_start+m_f; i++) {
94: /* convert to local "natural" numbering and then to PETSc global numbering */
95: row = idx_f[dof*(i-i_start_ghost)]/dof;
96: 97: i_c = (i/ratio); /* coarse grid node to left of fine grid node */
98: 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\
99: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
100: 101: /*
102: Only include those interpolation points that are truly
103: nonzero. Note this is very important for final grid lines
104: in x direction; since they have no right neighbor
105: */
106: x = ((double)(i - i_c*ratio))/((double)ratio);
107: nc = 0;
108: /* one left and below; or we are right on it */
109: col = dof*(i_c-i_start_ghost_c);
110: cols[nc] = idx_c[col]/dof;
111: v[nc++] = - x + 1.0;
112: /* one right? */
113: if (i_c*ratio != i) {
114: cols[nc] = idx_c[col+dof]/dof;
115: v[nc++] = x;
116: }
117: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
118: }
119: 120: } else {
121: PetscScalar *xi;
122: PetscInt li,nxi,n;
123: PetscScalar Ni[2];
124: 125: /* compute local coordinate arrays */
126: nxi = ratio + 1;
127: PetscMalloc(sizeof(PetscScalar)*nxi,&xi);
128: for (li=0; li<nxi; li++) {
129: xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
130: }
132: for (i=i_start; i<i_start+m_f; i++) {
133: /* convert to local "natural" numbering and then to PETSc global numbering */
134: row = idx_f[dof*(i-i_start_ghost)]/dof;
135: 136: i_c = (i/ratio); /* coarse grid node to left of fine grid node */
137: 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\
138: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
140: /* remainders */
141: li = i - ratio * (i/ratio);
142: if (i==mx-1){ li = nxi-1; }
143: 144: /* corners */
145: col = dof*(i_c-i_start_ghost_c);
146: cols[0] = idx_c[col]/dof;
147: Ni[0] = 1.0;
148: if ( (li==0) || (li==nxi-1) ) {
149: MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
150: continue;
151: }
152: 153: /* edges + interior */
154: /* remainders */
155: if (i==mx-1){ i_c--; }
156: 157: col = dof*(i_c-i_start_ghost_c);
158: cols[0] = idx_c[col]/dof; /* one left and below; or we are right on it */
159: cols[1] = idx_c[col+dof]/dof;
161: Ni[0] = 0.5*(1.0-xi[li]);
162: Ni[1] = 0.5*(1.0+xi[li]);
163: for (n=0; n<2; n++) {
164: if( PetscAbsScalar(Ni[n])<1.0e-32) { cols[n]=-1; }
165: }
166: MatSetValues(mat,1,&row,2,cols,Ni,INSERT_VALUES);
167: }
168: PetscFree(xi);
169: }
170: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
171: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
172: MatCreateMAIJ(mat,dof,A);
173: MatDestroy(&mat);
174: return(0);
175: }
179: PetscErrorCode DMCreateInterpolation_DA_1D_Q0(DM dac,DM daf,Mat *A)180: {
181: PetscErrorCode ierr;
182: PetscInt i,i_start,m_f,Mx,*idx_f;
183: PetscInt m_ghost,*idx_c,m_ghost_c;
184: PetscInt row,col,i_start_ghost,mx,m_c,nc,ratio;
185: PetscInt i_c,i_start_c,i_start_ghost_c,cols[2],dof;
186: PetscScalar v[2],x;
187: Mat mat;
188: DMDABoundaryType bx;
189: 191: DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);
192: DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);
193: if (bx == DMDA_BOUNDARY_PERIODIC){
194: ratio = mx/Mx;
195: 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);
196: } else {
197: ratio = (mx-1)/(Mx-1);
198: 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);
199: }
201: DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);
202: DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
203: DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
205: DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
206: DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
207: DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
209: /* create interpolation matrix */
210: MatCreate(((PetscObject)dac)->comm,&mat);
211: MatSetSizes(mat,m_f,m_c,mx,Mx);
212: MatSetType(mat,MATAIJ);
213: MatSeqAIJSetPreallocation(mat,2,PETSC_NULL);
214: MatMPIAIJSetPreallocation(mat,2,PETSC_NULL,0,PETSC_NULL);
216: /* loop over local fine grid nodes setting interpolation for those*/
217: for (i=i_start; i<i_start+m_f; i++) {
218: /* convert to local "natural" numbering and then to PETSc global numbering */
219: row = idx_f[dof*(i-i_start_ghost)]/dof;
221: i_c = (i/ratio); /* coarse grid node to left of fine grid node */
223: /*
224: Only include those interpolation points that are truly
225: nonzero. Note this is very important for final grid lines
226: in x direction; since they have no right neighbor
227: */
228: x = ((double)(i - i_c*ratio))/((double)ratio);
229: nc = 0;
230: /* one left and below; or we are right on it */
231: col = dof*(i_c-i_start_ghost_c);
232: cols[nc] = idx_c[col]/dof;
233: v[nc++] = - x + 1.0;
234: /* one right? */
235: if (i_c*ratio != i) {
236: cols[nc] = idx_c[col+dof]/dof;
237: v[nc++] = x;
238: }
239: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
240: }
241: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
242: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
243: MatCreateMAIJ(mat,dof,A);
244: MatDestroy(&mat);
245: PetscLogFlops(5.0*m_f);
246: return(0);
247: }
251: PetscErrorCode DMCreateInterpolation_DA_2D_Q1(DM dac,DM daf,Mat *A)252: {
253: PetscErrorCode ierr;
254: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof;
255: PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,*dnz,*onz;
256: PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
257: 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;
258: PetscMPIInt size_c,size_f,rank_f;
259: PetscScalar v[4],x,y;
260: Mat mat;
261: DMDABoundaryType bx,by;
262: 264: DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
265: DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
266: if (bx == DMDA_BOUNDARY_PERIODIC){
267: ratioi = mx/Mx;
268: 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);
269: } else {
270: ratioi = (mx-1)/(Mx-1);
271: 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);
272: }
273: if (by == DMDA_BOUNDARY_PERIODIC){
274: ratioj = my/My;
275: 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);
276: } else {
277: ratioj = (my-1)/(My-1);
278: 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);
279: }
280: 281: 282: DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
283: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
284: DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
285: 286: DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
287: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
288: DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
289: 290: /*
291: Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
292: The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
293: processors). It's effective length is hence 4 times its normal length, this is
294: why the col_scale is multiplied by the interpolation matrix column sizes.
295: sol_shift allows each set of 1/4 processors do its own interpolation using ITS
296: copy of the coarse vector. A bit of a hack but you do better.
297: 298: In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
299: */
300: MPI_Comm_size(((PetscObject)dac)->comm,&size_c);
301: MPI_Comm_size(((PetscObject)daf)->comm,&size_f);
302: MPI_Comm_rank(((PetscObject)daf)->comm,&rank_f);
303: col_scale = size_f/size_c;
304: col_shift = Mx*My*(rank_f/size_c);
305: 306: MatPreallocateInitialize(((PetscObject)daf)->comm,m_f*n_f,col_scale*m_c*n_c,dnz,onz);
307: for (j=j_start; j<j_start+n_f; j++) {
308: for (i=i_start; i<i_start+m_f; i++) {
309: /* convert to local "natural" numbering and then to PETSc global numbering */
310: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
311: 312: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
313: j_c = (j/ratioj); /* coarse grid node below fine grid node */
314: 315: 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\
316: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
317: 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\
318: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
319: 320: /*
321: Only include those interpolation points that are truly
322: nonzero. Note this is very important for final grid lines
323: in x and y directions; since they have no right/top neighbors
324: */
325: nc = 0;
326: /* one left and below; or we are right on it */
327: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
328: cols[nc++] = col_shift + idx_c[col]/dof;
329: /* one right and below */
330: if (i_c*ratioi != i) {
331: cols[nc++] = col_shift + idx_c[col+dof]/dof;
332: }
333: /* one left and above */
334: if (j_c*ratioj != j) {
335: cols[nc++] = col_shift + idx_c[col+m_ghost_c*dof]/dof;
336: }
337: /* one right and above */
338: if (i_c*ratioi != i && j_c*ratioj != j) {
339: cols[nc++] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof;
340: }
341: MatPreallocateSet(row,nc,cols,dnz,onz);
342: }
343: }
344: MatCreate(((PetscObject)daf)->comm,&mat);
345: MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);
346: MatSetType(mat,MATAIJ);
347: MatSeqAIJSetPreallocation(mat,0,dnz);
348: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
349: MatPreallocateFinalize(dnz,onz);
350: 351: /* loop over local fine grid nodes setting interpolation for those*/
352: if (!NEWVERSION) {
353: 354: for (j=j_start; j<j_start+n_f; j++) {
355: for (i=i_start; i<i_start+m_f; i++) {
356: /* convert to local "natural" numbering and then to PETSc global numbering */
357: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
358: 359: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
360: j_c = (j/ratioj); /* coarse grid node below fine grid node */
361: 362: /*
363: Only include those interpolation points that are truly
364: nonzero. Note this is very important for final grid lines
365: in x and y directions; since they have no right/top neighbors
366: */
367: x = ((double)(i - i_c*ratioi))/((double)ratioi);
368: y = ((double)(j - j_c*ratioj))/((double)ratioj);
369: 370: nc = 0;
371: /* one left and below; or we are right on it */
372: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
373: cols[nc] = col_shift + idx_c[col]/dof;
374: v[nc++] = x*y - x - y + 1.0;
375: /* one right and below */
376: if (i_c*ratioi != i) {
377: cols[nc] = col_shift + idx_c[col+dof]/dof;
378: v[nc++] = -x*y + x;
379: }
380: /* one left and above */
381: if (j_c*ratioj != j) {
382: cols[nc] = col_shift + idx_c[col+m_ghost_c*dof]/dof;
383: v[nc++] = -x*y + y;
384: }
385: /* one right and above */
386: if (j_c*ratioj != j && i_c*ratioi != i) {
387: cols[nc] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof;
388: v[nc++] = x*y;
389: }
390: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
391: }
392: }
393: 394: } else {
395: PetscScalar Ni[4];
396: PetscScalar *xi,*eta;
397: PetscInt li,nxi,lj,neta;
398: 399: /* compute local coordinate arrays */
400: nxi = ratioi + 1;
401: neta = ratioj + 1;
402: PetscMalloc(sizeof(PetscScalar)*nxi,&xi);
403: PetscMalloc(sizeof(PetscScalar)*neta,&eta);
404: for (li=0; li<nxi; li++) {
405: xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
406: }
407: for (lj=0; lj<neta; lj++) {
408: eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
409: }
411: /* loop over local fine grid nodes setting interpolation for those*/
412: for (j=j_start; j<j_start+n_f; j++) {
413: for (i=i_start; i<i_start+m_f; i++) {
414: /* convert to local "natural" numbering and then to PETSc global numbering */
415: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
416: 417: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
418: j_c = (j/ratioj); /* coarse grid node below fine grid node */
419: 420: /* remainders */
421: li = i - ratioi * (i/ratioi);
422: if (i==mx-1){ li = nxi-1; }
423: lj = j - ratioj * (j/ratioj);
424: if (j==my-1){ lj = neta-1; }
425: 426: /* corners */
427: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
428: cols[0] = col_shift + idx_c[col]/dof; /* left, below */
429: Ni[0] = 1.0;
430: if ( (li==0) || (li==nxi-1) ) {
431: if ( (lj==0) || (lj==neta-1) ) {
432: MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
433: continue;
434: }
435: }
436: 437: /* edges + interior */
438: /* remainders */
439: if (i==mx-1){ i_c--; }
440: if (j==my-1){ j_c--; }
442: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
443: cols[0] = col_shift + idx_c[col]/dof; /* left, below */
444: cols[1] = col_shift + idx_c[col+dof]/dof; /* right, below */
445: cols[2] = col_shift + idx_c[col+m_ghost_c*dof]/dof; /* left, above */
446: cols[3] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof; /* right, above */
448: Ni[0] = 0.25*(1.0-xi[li])*(1.0-eta[lj]);
449: Ni[1] = 0.25*(1.0+xi[li])*(1.0-eta[lj]);
450: Ni[2] = 0.25*(1.0-xi[li])*(1.0+eta[lj]);
451: Ni[3] = 0.25*(1.0+xi[li])*(1.0+eta[lj]);
453: nc = 0;
454: if( PetscAbsScalar(Ni[0])<1.0e-32) { cols[0]=-1; }
455: if( PetscAbsScalar(Ni[1])<1.0e-32) { cols[1]=-1; }
456: if( PetscAbsScalar(Ni[2])<1.0e-32) { cols[2]=-1; }
457: if( PetscAbsScalar(Ni[3])<1.0e-32) { cols[3]=-1; }
458: 459: MatSetValues(mat,1,&row,4,cols,Ni,INSERT_VALUES);
460: }
461: }
462: PetscFree(xi);
463: PetscFree(eta);
464: }
465: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
466: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
467: MatCreateMAIJ(mat,dof,A);
468: MatDestroy(&mat);
469: return(0);
470: }
472: /*
473: Contributed by Andrei Draganescu <aidraga@sandia.gov>
474: */
477: PetscErrorCode DMCreateInterpolation_DA_2D_Q0(DM dac,DM daf,Mat *A)478: {
479: PetscErrorCode ierr;
480: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof;
481: PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,*dnz,*onz;
482: PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
483: 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;
484: PetscMPIInt size_c,size_f,rank_f;
485: PetscScalar v[4];
486: Mat mat;
487: DMDABoundaryType bx,by;
490: DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
491: DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
492: if (bx == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in x");
493: if (by == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in y");
494: ratioi = mx/Mx;
495: ratioj = my/My;
496: if (ratioi*Mx != mx) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x");
497: if (ratioj*My != my) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y");
498: if (ratioi != 2) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 2");
499: if (ratioj != 2) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 2");
501: DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
502: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
503: DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
505: DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
506: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
507: DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
509: /*
510: Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
511: The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
512: processors). It's effective length is hence 4 times its normal length, this is
513: why the col_scale is multiplied by the interpolation matrix column sizes.
514: sol_shift allows each set of 1/4 processors do its own interpolation using ITS
515: copy of the coarse vector. A bit of a hack but you do better.
517: In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
518: */
519: MPI_Comm_size(((PetscObject)dac)->comm,&size_c);
520: MPI_Comm_size(((PetscObject)daf)->comm,&size_f);
521: MPI_Comm_rank(((PetscObject)daf)->comm,&rank_f);
522: col_scale = size_f/size_c;
523: col_shift = Mx*My*(rank_f/size_c);
525: MatPreallocateInitialize(((PetscObject)daf)->comm,m_f*n_f,col_scale*m_c*n_c,dnz,onz);
526: for (j=j_start; j<j_start+n_f; j++) {
527: for (i=i_start; i<i_start+m_f; i++) {
528: /* convert to local "natural" numbering and then to PETSc global numbering */
529: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
531: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
532: j_c = (j/ratioj); /* coarse grid node below fine grid node */
534: 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\
535: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
536: 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\
537: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
539: /*
540: Only include those interpolation points that are truly
541: nonzero. Note this is very important for final grid lines
542: in x and y directions; since they have no right/top neighbors
543: */
544: nc = 0;
545: /* one left and below; or we are right on it */
546: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
547: cols[nc++] = col_shift + idx_c[col]/dof;
548: MatPreallocateSet(row,nc,cols,dnz,onz);
549: }
550: }
551: MatCreate(((PetscObject)daf)->comm,&mat);
552: MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);
553: MatSetType(mat,MATAIJ);
554: MatSeqAIJSetPreallocation(mat,0,dnz);
555: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
556: MatPreallocateFinalize(dnz,onz);
558: /* loop over local fine grid nodes setting interpolation for those*/
559: for (j=j_start; j<j_start+n_f; j++) {
560: for (i=i_start; i<i_start+m_f; i++) {
561: /* convert to local "natural" numbering and then to PETSc global numbering */
562: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
564: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
565: j_c = (j/ratioj); /* coarse grid node below fine grid node */
566: nc = 0;
567: /* one left and below; or we are right on it */
568: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
569: cols[nc] = col_shift + idx_c[col]/dof;
570: v[nc++] = 1.0;
571: 572: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
573: }
574: }
575: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
576: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
577: MatCreateMAIJ(mat,dof,A);
578: MatDestroy(&mat);
579: PetscLogFlops(13.0*m_f*n_f);
580: return(0);
581: }
583: /*
584: Contributed by Jianming Yang <jianming-yang@uiowa.edu>
585: */
588: PetscErrorCode DMCreateInterpolation_DA_3D_Q0(DM dac,DM daf,Mat *A)589: {
590: PetscErrorCode ierr;
591: PetscInt i,j,l,i_start,j_start,l_start,m_f,n_f,p_f,Mx,My,Mz,*idx_f,dof;
592: PetscInt m_ghost,n_ghost,p_ghost,*idx_c,m_ghost_c,n_ghost_c,p_ghost_c,nc,*dnz,*onz;
593: 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;
594: 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;
595: PetscMPIInt size_c,size_f,rank_f;
596: PetscScalar v[8];
597: Mat mat;
598: DMDABoundaryType bx,by,bz;
599: 601: DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
602: DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
603: if (bx == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in x");
604: if (by == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in y");
605: if (bz == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in z");
606: ratioi = mx/Mx;
607: ratioj = my/My;
608: ratiol = mz/Mz;
609: if (ratioi*Mx != mx) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x");
610: if (ratioj*My != my) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y");
611: if (ratiol*Mz != mz) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in z");
612: if (ratioi != 2 && ratioi != 1) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 1 or 2");
613: if (ratioj != 2 && ratioj != 1) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 1 or 2");
614: if (ratiol != 2 && ratiol != 1) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in z must be 1 or 2");
616: DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
617: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
618: DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
620: DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
621: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
622: DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
623: /*
624: Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
625: The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
626: processors). It's effective length is hence 4 times its normal length, this is
627: why the col_scale is multiplied by the interpolation matrix column sizes.
628: sol_shift allows each set of 1/4 processors do its own interpolation using ITS
629: copy of the coarse vector. A bit of a hack but you do better.
631: In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
632: */
633: MPI_Comm_size(((PetscObject)dac)->comm,&size_c);
634: MPI_Comm_size(((PetscObject)daf)->comm,&size_f);
635: MPI_Comm_rank(((PetscObject)daf)->comm,&rank_f);
636: col_scale = size_f/size_c;
637: col_shift = Mx*My*Mz*(rank_f/size_c);
639: MatPreallocateInitialize(((PetscObject)daf)->comm,m_f*n_f*p_f,col_scale*m_c*n_c*p_c,dnz,onz);
640: for (l=l_start; l<l_start+p_f; l++) {
641: for (j=j_start; j<j_start+n_f; j++) {
642: for (i=i_start; i<i_start+m_f; i++) {
643: /* convert to local "natural" numbering and then to PETSc global numbering */
644: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
646: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
647: j_c = (j/ratioj); /* coarse grid node below fine grid node */
648: l_c = (l/ratiol);
650: 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\
651: l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
652: 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\
653: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
654: 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\
655: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
657: /*
658: Only include those interpolation points that are truly
659: nonzero. Note this is very important for final grid lines
660: in x and y directions; since they have no right/top neighbors
661: */
662: nc = 0;
663: /* one left and below; or we are right on it */
664: col = dof*(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));
665: cols[nc++] = col_shift + idx_c[col]/dof;
666: MatPreallocateSet(row,nc,cols,dnz,onz);
667: }
668: }
669: }
670: MatCreate(((PetscObject)daf)->comm,&mat);
671: MatSetSizes(mat,m_f*n_f*p_f,col_scale*m_c*n_c*p_c,mx*my*mz,col_scale*Mx*My*Mz);
672: MatSetType(mat,MATAIJ);
673: MatSeqAIJSetPreallocation(mat,0,dnz);
674: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
675: MatPreallocateFinalize(dnz,onz);
677: /* loop over local fine grid nodes setting interpolation for those*/
678: for (l=l_start; l<l_start+p_f; l++) {
679: for (j=j_start; j<j_start+n_f; j++) {
680: for (i=i_start; i<i_start+m_f; i++) {
681: /* convert to local "natural" numbering and then to PETSc global numbering */
682: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
683: 684: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
685: j_c = (j/ratioj); /* coarse grid node below fine grid node */
686: l_c = (l/ratiol);
687: nc = 0;
688: /* one left and below; or we are right on it */
689: col = dof*(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));
690: cols[nc] = col_shift + idx_c[col]/dof;
691: v[nc++] = 1.0;
692: 693: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
694: }
695: }
696: }
697: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
698: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
699: MatCreateMAIJ(mat,dof,A);
700: MatDestroy(&mat);
701: PetscLogFlops(13.0*m_f*n_f*p_f);
702: return(0);
703: }
707: PetscErrorCode DMCreateInterpolation_DA_3D_Q1(DM dac,DM daf,Mat *A)708: {
709: PetscErrorCode ierr;
710: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof,l;
711: PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,Mz,mz;
712: PetscInt row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok;
713: PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
714: PetscInt l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c;
715: PetscInt l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz;
716: PetscScalar v[8],x,y,z;
717: Mat mat;
718: DMDABoundaryType bx,by,bz;
719: 721: DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
722: DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
723: if (mx == Mx) {
724: ratioi = 1;
725: } else if (bx == DMDA_BOUNDARY_PERIODIC) {
726: ratioi = mx/Mx;
727: 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);
728: } else {
729: ratioi = (mx-1)/(Mx-1);
730: 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);
731: }
732: if (my == My) {
733: ratioj = 1;
734: } else if (by == DMDA_BOUNDARY_PERIODIC) {
735: ratioj = my/My;
736: 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);
737: } else {
738: ratioj = (my-1)/(My-1);
739: 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);
740: }
741: if (mz == Mz) {
742: ratiok = 1;
743: } else if (bz == DMDA_BOUNDARY_PERIODIC) {
744: ratiok = mz/Mz;
745: 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);
746: } else {
747: ratiok = (mz-1)/(Mz-1);
748: 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);
749: }
750: 751: DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
752: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
753: DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
754: 755: DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
756: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
757: DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
758: 759: /* create interpolation matrix, determining exact preallocation */
760: MatPreallocateInitialize(((PetscObject)dac)->comm,m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);
761: /* loop over local fine grid nodes counting interpolating points */
762: for (l=l_start; l<l_start+p_f; l++) {
763: for (j=j_start; j<j_start+n_f; j++) {
764: for (i=i_start; i<i_start+m_f; i++) {
765: /* convert to local "natural" numbering and then to PETSc global numbering */
766: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
767: i_c = (i/ratioi);
768: j_c = (j/ratioj);
769: l_c = (l/ratiok);
770: 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\
771: l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
772: 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\
773: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
774: 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\
775: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
776: 777: /*
778: Only include those interpolation points that are truly
779: nonzero. Note this is very important for final grid lines
780: in x and y directions; since they have no right/top neighbors
781: */
782: nc = 0;
783: col = dof*(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));
784: cols[nc++] = idx_c[col]/dof;
785: if (i_c*ratioi != i) {
786: cols[nc++] = idx_c[col+dof]/dof;
787: }
788: if (j_c*ratioj != j) {
789: cols[nc++] = idx_c[col+m_ghost_c*dof]/dof;
790: }
791: if (l_c*ratiok != l) {
792: cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof;
793: }
794: if (j_c*ratioj != j && i_c*ratioi != i) {
795: cols[nc++] = idx_c[col+(m_ghost_c+1)*dof]/dof;
796: }
797: if (j_c*ratioj != j && l_c*ratiok != l) {
798: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;
799: }
800: if (i_c*ratioi != i && l_c*ratiok != l) {
801: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof;
802: }
803: if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
804: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof;
805: }
806: MatPreallocateSet(row,nc,cols,dnz,onz);
807: }
808: }
809: }
810: MatCreate(((PetscObject)dac)->comm,&mat);
811: MatSetSizes(mat,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz);
812: MatSetType(mat,MATAIJ);
813: MatSeqAIJSetPreallocation(mat,0,dnz);
814: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
815: MatPreallocateFinalize(dnz,onz);
816: 817: /* loop over local fine grid nodes setting interpolation for those*/
818: if (!NEWVERSION) {
820: for (l=l_start; l<l_start+p_f; l++) {
821: for (j=j_start; j<j_start+n_f; j++) {
822: for (i=i_start; i<i_start+m_f; i++) {
823: /* convert to local "natural" numbering and then to PETSc global numbering */
824: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
825: 826: i_c = (i/ratioi);
827: j_c = (j/ratioj);
828: l_c = (l/ratiok);
829: 830: /*
831: Only include those interpolation points that are truly
832: nonzero. Note this is very important for final grid lines
833: in x and y directions; since they have no right/top neighbors
834: */
835: x = ((double)(i - i_c*ratioi))/((double)ratioi);
836: y = ((double)(j - j_c*ratioj))/((double)ratioj);
837: z = ((double)(l - l_c*ratiok))/((double)ratiok);
839: nc = 0;
840: /* one left and below; or we are right on it */
841: col = dof*(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));
842: 843: cols[nc] = idx_c[col]/dof;
844: v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
845: 846: if (i_c*ratioi != i) {
847: cols[nc] = idx_c[col+dof]/dof;
848: v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
849: }
850: 851: if (j_c*ratioj != j) {
852: cols[nc] = idx_c[col+m_ghost_c*dof]/dof;
853: v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
854: }
855: 856: if (l_c*ratiok != l) {
857: cols[nc] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof;
858: v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
859: }
860: 861: if (j_c*ratioj != j && i_c*ratioi != i) {
862: cols[nc] = idx_c[col+(m_ghost_c+1)*dof]/dof;
863: v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
864: }
865: 866: if (j_c*ratioj != j && l_c*ratiok != l) {
867: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;
868: v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
869: }
870: 871: if (i_c*ratioi != i && l_c*ratiok != l) {
872: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof;
873: v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
874: }
875: 876: if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
877: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof;
878: v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
879: }
880: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
881: }
882: }
883: }
884: 885: } else {
886: PetscScalar *xi,*eta,*zeta;
887: PetscInt li,nxi,lj,neta,lk,nzeta,n;
888: PetscScalar Ni[8];
889: 890: /* compute local coordinate arrays */
891: nxi = ratioi + 1;
892: neta = ratioj + 1;
893: nzeta = ratiok + 1;
894: PetscMalloc(sizeof(PetscScalar)*nxi,&xi);
895: PetscMalloc(sizeof(PetscScalar)*neta,&eta);
896: PetscMalloc(sizeof(PetscScalar)*nzeta,&zeta);
897: for (li=0; li<nxi; li++) {
898: xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
899: }
900: for (lj=0; lj<neta; lj++) {
901: eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
902: }
903: for (lk=0; lk<nzeta; lk++) {
904: zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));
905: }
906: 907: for (l=l_start; l<l_start+p_f; l++) {
908: for (j=j_start; j<j_start+n_f; j++) {
909: for (i=i_start; i<i_start+m_f; i++) {
910: /* convert to local "natural" numbering and then to PETSc global numbering */
911: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
912: 913: i_c = (i/ratioi);
914: j_c = (j/ratioj);
915: l_c = (l/ratiok);
917: /* remainders */
918: li = i - ratioi * (i/ratioi);
919: if (i==mx-1){ li = nxi-1; }
920: lj = j - ratioj * (j/ratioj);
921: if (j==my-1){ lj = neta-1; }
922: lk = l - ratiok * (l/ratiok);
923: if (l==mz-1){ lk = nzeta-1; }
924: 925: /* corners */
926: col = dof*(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));
927: cols[0] = idx_c[col]/dof;
928: Ni[0] = 1.0;
929: if ( (li==0) || (li==nxi-1) ) {
930: if ( (lj==0) || (lj==neta-1) ) {
931: if ( (lk==0) || (lk==nzeta-1) ) {
932: MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
933: continue;
934: }
935: }
936: }
937: 938: /* edges + interior */
939: /* remainders */
940: if (i==mx-1){ i_c--; }
941: if (j==my-1){ j_c--; }
942: if (l==mz-1){ l_c--; }
943: 944: col = dof*(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));
945: cols[0] = idx_c[col]/dof; /* one left and below; or we are right on it */
946: cols[1] = idx_c[col+dof]/dof; /* one right and below */
947: cols[2] = idx_c[col+m_ghost_c*dof]/dof; /* one left and above */
948: cols[3] = idx_c[col+(m_ghost_c+1)*dof]/dof; /* one right and above */
950: cols[4] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof; /* one left and below and front; or we are right on it */
951: cols[5] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof; /* one right and below, and front */
952: cols[6] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;/* one left and above and front*/
953: cols[7] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof; /* one right and above and front */
955: Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
956: Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
957: Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
958: Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
960: Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
961: Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
962: Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
963: Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
965: for (n=0; n<8; n++) {
966: if( PetscAbsScalar(Ni[n])<1.0e-32) { cols[n]=-1; }
967: }
968: MatSetValues(mat,1,&row,8,cols,Ni,INSERT_VALUES);
969: 970: }
971: }
972: }
973: PetscFree(xi);
974: PetscFree(eta);
975: PetscFree(zeta);
976: }
977: 978: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
979: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
980: 981: MatCreateMAIJ(mat,dof,A);
982: MatDestroy(&mat);
983: return(0);
984: }
988: PetscErrorCode DMCreateInterpolation_DA(DM dac,DM daf,Mat *A,Vec *scale)989: {
990: PetscErrorCode ierr;
991: PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
992: DMDABoundaryType bxc,byc,bzc,bxf,byf,bzf;
993: DMDAStencilType stc,stf;
994: DM_DA *ddc = (DM_DA*)dac->data;
1002: DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1003: DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1004: if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1005: if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1006: if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1007: if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1008: if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1009: if (Mc < 2 && Mf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
1010: 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");
1011: 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");
1013: if (ddc->interptype == DMDA_Q1){
1014: if (dimc == 1){
1015: DMCreateInterpolation_DA_1D_Q1(dac,daf,A);
1016: } else if (dimc == 2){
1017: DMCreateInterpolation_DA_2D_Q1(dac,daf,A);
1018: } else if (dimc == 3){
1019: DMCreateInterpolation_DA_3D_Q1(dac,daf,A);
1020: } else SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype);
1021: } else if (ddc->interptype == DMDA_Q0){
1022: if (dimc == 1){
1023: DMCreateInterpolation_DA_1D_Q0(dac,daf,A);
1024: } else if (dimc == 2){
1025: DMCreateInterpolation_DA_2D_Q0(dac,daf,A);
1026: } else if (dimc == 3){
1027: DMCreateInterpolation_DA_3D_Q0(dac,daf,A);
1028: } else SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype);
1029: }
1030: if (scale) {
1031: DMCreateInterpolationScale((DM)dac,(DM)daf,*A,scale);
1032: }
1033: return(0);
1034: }
1038: PetscErrorCode DMCreateInjection_DA_1D(DM dac,DM daf,VecScatter *inject)1039: {
1040: PetscErrorCode ierr;
1041: PetscInt i,i_start,m_f,Mx,*idx_f,dof;
1042: PetscInt m_ghost,*idx_c,m_ghost_c;
1043: PetscInt row,i_start_ghost,mx,m_c,nc,ratioi;
1044: PetscInt i_start_c,i_start_ghost_c;
1045: PetscInt *cols;
1046: DMDABoundaryType bx;
1047: Vec vecf,vecc;
1048: IS isf;
1049: 1051: DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);
1052: DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);
1053: if (bx == DMDA_BOUNDARY_PERIODIC) {
1054: ratioi = mx/Mx;
1055: 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);
1056: } else {
1057: ratioi = (mx-1)/(Mx-1);
1058: 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);
1059: }
1060: 1061: DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);
1062: DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
1063: DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
1064: 1065: DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
1066: DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
1067: DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
1068: 1069: 1070: /* loop over local fine grid nodes setting interpolation for those*/
1071: nc = 0;
1072: PetscMalloc(m_f*sizeof(PetscInt),&cols);
1073: 1074: 1075: for (i=i_start_c; i<i_start_c+m_c; i++) {
1076: PetscInt i_f = i*ratioi;
1078: 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\
1079: i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1080: row = idx_f[dof*(i_f-i_start_ghost)];
1081: cols[nc++] = row/dof;
1082: }
1083: 1085: ISCreateBlock(((PetscObject)daf)->comm,dof,nc,cols,PETSC_OWN_POINTER,&isf);
1086: DMGetGlobalVector(dac,&vecc);
1087: DMGetGlobalVector(daf,&vecf);
1088: VecScatterCreate(vecf,isf,vecc,PETSC_NULL,inject);
1089: DMRestoreGlobalVector(dac,&vecc);
1090: DMRestoreGlobalVector(daf,&vecf);
1091: ISDestroy(&isf);
1092: return(0);
1093: }
1097: PetscErrorCode DMCreateInjection_DA_2D(DM dac,DM daf,VecScatter *inject)1098: {
1099: PetscErrorCode ierr;
1100: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof;
1101: PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c;
1102: PetscInt row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj;
1103: PetscInt i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
1104: PetscInt *cols;
1105: DMDABoundaryType bx,by;
1106: Vec vecf,vecc;
1107: IS isf;
1110: DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
1111: DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
1112: if (bx == DMDA_BOUNDARY_PERIODIC) {
1113: ratioi = mx/Mx;
1114: 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);
1115: } else {
1116: ratioi = (mx-1)/(Mx-1);
1117: 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);
1118: }
1119: if (by == DMDA_BOUNDARY_PERIODIC) {
1120: ratioj = my/My;
1121: 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);
1122: } else {
1123: ratioj = (my-1)/(My-1);
1124: 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);
1125: }
1127: DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
1128: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
1129: DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
1131: DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
1132: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
1133: DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
1136: /* loop over local fine grid nodes setting interpolation for those*/
1137: nc = 0;
1138: PetscMalloc(n_f*m_f*sizeof(PetscInt),&cols);
1139: for (j=j_start_c; j<j_start_c+n_c; j++) {
1140: for (i=i_start_c; i<i_start_c+m_c; i++) {
1141: PetscInt i_f = i*ratioi,j_f = j*ratioj;
1142: 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\
1143: j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost);
1144: 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\
1145: i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1146: row = idx_f[dof*(m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))];
1147: cols[nc++] = row/dof;
1148: }
1149: }
1151: ISCreateBlock(((PetscObject)daf)->comm,dof,nc,cols,PETSC_OWN_POINTER,&isf);
1152: DMGetGlobalVector(dac,&vecc);
1153: DMGetGlobalVector(daf,&vecf);
1154: VecScatterCreate(vecf,isf,vecc,PETSC_NULL,inject);
1155: DMRestoreGlobalVector(dac,&vecc);
1156: DMRestoreGlobalVector(daf,&vecf);
1157: ISDestroy(&isf);
1158: return(0);
1159: }
1163: PetscErrorCode DMCreateInjection_DA_3D(DM dac,DM daf,VecScatter *inject)1164: {
1165: PetscErrorCode ierr;
1166: PetscInt i,j,k,i_start,j_start,k_start,m_f,n_f,p_f,Mx,My,Mz;
1167: PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c;
1168: PetscInt i_start_ghost,j_start_ghost,k_start_ghost;
1169: PetscInt mx,my,mz,ratioi,ratioj,ratiok;
1170: PetscInt i_start_c,j_start_c,k_start_c;
1171: PetscInt m_c,n_c,p_c;
1172: PetscInt i_start_ghost_c,j_start_ghost_c,k_start_ghost_c;
1173: PetscInt row,nc,dof;
1174: PetscInt *idx_c,*idx_f;
1175: PetscInt *cols;
1176: DMDABoundaryType bx,by,bz;
1177: Vec vecf,vecc;
1178: IS isf;
1181: DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
1182: DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
1184: if (bx == DMDA_BOUNDARY_PERIODIC){
1185: ratioi = mx/Mx;
1186: 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);
1187: } else {
1188: ratioi = (mx-1)/(Mx-1);
1189: 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);
1190: }
1191: if (by == DMDA_BOUNDARY_PERIODIC){
1192: ratioj = my/My;
1193: 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);
1194: } else {
1195: ratioj = (my-1)/(My-1);
1196: 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);
1197: }
1198: if (bz == DMDA_BOUNDARY_PERIODIC){
1199: ratiok = mz/Mz;
1200: 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);
1201: } else {
1202: ratiok = (mz-1)/(Mz-1);
1203: 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);
1204: }
1206: DMDAGetCorners(daf,&i_start,&j_start,&k_start,&m_f,&n_f,&p_f);
1207: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&k_start_ghost,&m_ghost,&n_ghost,&p_ghost);
1208: DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
1210: DMDAGetCorners(dac,&i_start_c,&j_start_c,&k_start_c,&m_c,&n_c,&p_c);
1211: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&k_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
1212: DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
1215: /* loop over local fine grid nodes setting interpolation for those*/
1216: nc = 0;
1217: PetscMalloc(n_f*m_f*p_f*sizeof(PetscInt),&cols);
1218: for (k=k_start_c; k<k_start_c+p_c; k++) {
1219: for (j=j_start_c; j<j_start_c+n_c; j++) {
1220: for (i=i_start_c; i<i_start_c+m_c; i++) {
1221: PetscInt i_f = i*ratioi,j_f = j*ratioj,k_f = k*ratiok;
1222: 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 "
1223: "k_c %D k_f %D fine ghost range [%D,%D]",k,k_f,k_start_ghost,k_start_ghost+p_ghost);
1224: 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 "
1225: "j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost);
1226: 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 "
1227: "i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1228: row = idx_f[dof*(m_ghost*n_ghost*(k_f-k_start_ghost) + m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))];
1229: cols[nc++] = row/dof;
1230: }
1231: }
1232: }
1234: ISCreateBlock(((PetscObject)daf)->comm,dof,nc,cols,PETSC_OWN_POINTER,&isf);
1235: DMGetGlobalVector(dac,&vecc);
1236: DMGetGlobalVector(daf,&vecf);
1237: VecScatterCreate(vecf,isf,vecc,PETSC_NULL,inject);
1238: DMRestoreGlobalVector(dac,&vecc);
1239: DMRestoreGlobalVector(daf,&vecf);
1240: ISDestroy(&isf);
1241: return(0);
1242: }
1246: PetscErrorCode DMCreateInjection_DA(DM dac,DM daf,VecScatter *inject)1247: {
1248: PetscErrorCode ierr;
1249: PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1250: DMDABoundaryType bxc,byc,bzc,bxf,byf,bzf;
1251: DMDAStencilType stc,stf;
1258: DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1259: DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1260: if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1261: if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1262: if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1263: if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1264: if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1265: if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
1266: if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction");
1267: if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction");
1269: if (dimc == 1){
1270: DMCreateInjection_DA_1D(dac,daf,inject);
1271: } else if (dimc == 2) {
1272: DMCreateInjection_DA_2D(dac,daf,inject);
1273: } else if (dimc == 3) {
1274: DMCreateInjection_DA_3D(dac,daf,inject);
1275: }
1276: return(0);
1277: }
1281: PetscErrorCode DMCreateAggregates_DA(DM dac,DM daf,Mat *rest)1282: {
1283: PetscErrorCode ierr;
1284: PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc;
1285: PetscInt dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1286: DMDABoundaryType bxc,byc,bzc,bxf,byf,bzf;
1287: DMDAStencilType stc,stf;
1288: PetscInt i,j,l;
1289: PetscInt i_start,j_start,l_start, m_f,n_f,p_f;
1290: PetscInt i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost;
1291: PetscInt *idx_f;
1292: PetscInt i_c,j_c,l_c;
1293: PetscInt i_start_c,j_start_c,l_start_c, m_c,n_c,p_c;
1294: PetscInt i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c;
1295: PetscInt *idx_c;
1296: PetscInt d;
1297: PetscInt a;
1298: PetscInt max_agg_size;
1299: PetscInt *fine_nodes;
1300: PetscScalar *one_vec;
1301: PetscInt fn_idx;
1308: DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1309: DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1310: if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1311: if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1312: if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1313: if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1314: if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1316: 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);
1317: 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);
1318: 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);
1320: if (Pc < 0) Pc = 1;
1321: if (Pf < 0) Pf = 1;
1322: if (Nc < 0) Nc = 1;
1323: if (Nf < 0) Nf = 1;
1325: DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
1326: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
1327: DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
1329: DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
1330: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
1331: DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
1333: /*
1334: Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios
1335: for dimension 1 and 2 respectively.
1336: Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1)
1337: and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate.
1338: Each specific dof on the fine grid is mapped to one dof on the coarse grid.
1339: */
1341: max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1);
1343: /* create the matrix that will contain the restriction operator */
1344: MatCreateAIJ( ((PetscObject)daf)->comm, m_c*n_c*p_c*dofc, m_f*n_f*p_f*doff, Mc*Nc*Pc*dofc, Mf*Nf*Pf*doff,
1345: max_agg_size, PETSC_NULL, max_agg_size, PETSC_NULL, rest);
1347: /* store nodes in the fine grid here */
1348: PetscMalloc2(max_agg_size,PetscScalar, &one_vec,max_agg_size,PetscInt, &fine_nodes);
1349: for(i=0; i<max_agg_size; i++) one_vec[i] = 1.0;
1350: 1351: /* loop over all coarse nodes */
1352: for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) {
1353: for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) {
1354: for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) {
1355: for(d=0; d<dofc; d++) {
1356: /* convert to local "natural" numbering and then to PETSc global numbering */
1357: 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;
1359: fn_idx = 0;
1360: /* Corresponding fine points are all points (i_f, j_f, l_f) such that
1361: i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc
1362: (same for other dimensions)
1363: */
1364: for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) {
1365: for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) {
1366: for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) {
1367: 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;
1368: fn_idx++;
1369: }
1370: }
1371: }
1372: /* add all these points to one aggregate */
1373: MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);
1374: }
1375: }
1376: }
1377: }
1378: PetscFree2(one_vec,fine_nodes);
1379: MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);
1380: MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);
1381: return(0);
1382: }