Actual source code: dainterp.c
1: /*$Id: dainterp.c,v 1.25 2001/08/07 03:04:39 balay Exp $*/
2:
3: /*
4: Code for interpolating between grids represented by DAs
5: */
7: #include src/dm/da/daimpl.h
8: #include petscmg.h
10: #undef __FUNCT__
12: int DMGetInterpolationScale(DM dac,DM daf,Mat mat,Vec *scale)
13: {
14: int ierr;
15: Vec fine;
16: PetscScalar one = 1.0;
19: DMCreateGlobalVector(daf,&fine);
20: DMCreateGlobalVector(dac,scale);
21: VecSet(&one,fine);
22: MatRestrict(mat,fine,*scale);
23: VecDestroy(fine);
24: VecReciprocal(*scale);
25: return(0);
26: }
28: #undef __FUNCT__
30: int DAGetInterpolation_1D_Q1(DA dac,DA daf,Mat *A)
31: {
32: int ierr,i,i_start,m_f,Mx,*idx_f;
33: int m_ghost,*idx_c,m_ghost_c;
34: int row,col,i_start_ghost,mx,m_c,nc,ratio;
35: int i_c,i_start_c,i_start_ghost_c,cols[2],dof;
36: PetscScalar v[2],x;
37: Mat mat;
38: DAPeriodicType pt;
41: DAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&pt,0);
42: DAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0);
43: if (pt == DA_XPERIODIC) {
44: ratio = mx/Mx;
45: if (ratio*Mx != mx) SETERRQ2(1,"Ratio between levels: mx/Mx must be integer: mx %d Mx %d",mx,Mx);
46: } else {
47: ratio = (mx-1)/(Mx-1);
48: if (ratio*(Mx-1) != mx-1) SETERRQ2(1,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %d Mx %d",mx,Mx);
49: }
51: DAGetCorners(daf,&i_start,0,0,&m_f,0,0);
52: DAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
53: DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
55: DAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
56: DAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
57: DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
59: /* create interpolation matrix */
60: MatCreateMPIAIJ(dac->comm,m_f,m_c,mx,Mx,2,0,0,0,&mat);
61: if (!DAXPeriodic(pt)){MatSetOption(mat,MAT_COLUMNS_SORTED);}
63: /* loop over local fine grid nodes setting interpolation for those*/
64: for (i=i_start; i<i_start+m_f; i++) {
65: /* convert to local "natural" numbering and then to PETSc global numbering */
66: row = idx_f[dof*(i-i_start_ghost)]/dof;
68: i_c = (i/ratio); /* coarse grid node to left of fine grid node */
69: if (i_c < i_start_ghost_c) SETERRQ3(1,"Processor's coarse DA must lie over fine DAn
70: i_start %d i_c %d i_start_ghost_c %d",i_start,i_c,i_start_ghost_c);
72: /*
73: Only include those interpolation points that are truly
74: nonzero. Note this is very important for final grid lines
75: in x direction; since they have no right neighbor
76: */
77: x = ((double)(i - i_c*ratio))/((double)ratio);
78: /* printf("i j %d %d %g %gn",i,j,x,y); */
79: nc = 0;
80: /* one left and below; or we are right on it */
81: col = dof*(i_c-i_start_ghost_c);
82: cols[nc] = idx_c[col]/dof;
83: v[nc++] = - x + 1.0;
84: /* one right? */
85: if (i_c*ratio != i) {
86: cols[nc] = idx_c[col+dof]/dof;
87: v[nc++] = x;
88: }
89: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
90: }
91: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
92: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
93: MatCreateMAIJ(mat,dof,A);
94: MatDestroy(mat);
95: PetscLogFlops(5*m_f);
96: return(0);
97: }
99: #undef __FUNCT__
101: int DAGetInterpolation_1D_Q0(DA dac,DA daf,Mat *A)
102: {
103: int ierr,i,i_start,m_f,Mx,*idx_f;
104: int m_ghost,*idx_c,m_ghost_c;
105: int row,col,i_start_ghost,mx,m_c,nc,ratio;
106: int i_c,i_start_c,i_start_ghost_c,cols[2],dof;
107: PetscScalar v[2],x;
108: Mat mat;
109: DAPeriodicType pt;
112: DAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&pt,0);
113: DAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0);
114: if (pt == DA_XPERIODIC) {
115: ratio = mx/Mx;
116: if (ratio*Mx != mx) SETERRQ2(1,"Ratio between levels: mx/Mx must be integer: mx %d Mx %d",mx,Mx);
117: } else {
118: ratio = (mx-1)/(Mx-1);
119: if (ratio*(Mx-1) != mx-1) SETERRQ2(1,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %d Mx %d",mx,Mx);
120: }
122: DAGetCorners(daf,&i_start,0,0,&m_f,0,0);
123: DAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
124: DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
126: DAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
127: DAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
128: DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
130: /* create interpolation matrix */
131: MatCreateMPIAIJ(dac->comm,m_f,m_c,mx,Mx,2,0,0,0,&mat);
132: if (!DAXPeriodic(pt)) {MatSetOption(mat,MAT_COLUMNS_SORTED);}
134: /* loop over local fine grid nodes setting interpolation for those*/
135: for (i=i_start; i<i_start+m_f; i++) {
136: /* convert to local "natural" numbering and then to PETSc global numbering */
137: row = idx_f[dof*(i-i_start_ghost)]/dof;
139: i_c = (i/ratio); /* coarse grid node to left of fine grid node */
141: /*
142: Only include those interpolation points that are truly
143: nonzero. Note this is very important for final grid lines
144: in x direction; since they have no right neighbor
145: */
146: x = ((double)(i - i_c*ratio))/((double)ratio);
147: /* printf("i j %d %d %g %gn",i,j,x,y); */
148: nc = 0;
149: /* one left and below; or we are right on it */
150: col = dof*(i_c-i_start_ghost_c);
151: cols[nc] = idx_c[col]/dof;
152: v[nc++] = - x + 1.0;
153: /* one right? */
154: if (i_c*ratio != i) {
155: cols[nc] = idx_c[col+dof]/dof;
156: v[nc++] = x;
157: }
158: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
159: }
160: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
161: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
162: MatCreateMAIJ(mat,dof,A);
163: MatDestroy(mat);
164: PetscLogFlops(5*m_f);
165: return(0);
166: }
169: /* dof degree of freedom per node, nonperiodic */
170: #undef __FUNCT__
172: int DAGetInterpolation_2D_Q1(DA dac,DA daf,Mat *A)
173: {
174: int ierr,i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof;
175: int m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,*dnz,*onz;
176: int row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
177: int i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
178: int size_c,size_f,rank_f,col_shift,col_scale;
179: PetscScalar v[4],x,y;
180: Mat mat;
181: DAPeriodicType pt;
185: DAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&pt,0);
186: DAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0);
187: if (DAXPeriodic(pt)){
188: ratioi = mx/Mx;
189: if (ratioi*Mx != mx) SETERRQ2(1,"Ratio between levels: mx/Mx must be integer: mx %d Mx %d",mx,Mx);
190: } else {
191: ratioi = (mx-1)/(Mx-1);
192: if (ratioi*(Mx-1) != mx-1) SETERRQ2(1,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %d Mx %d",mx,Mx);
193: }
194: if (DAYPeriodic(pt)){
195: ratioj = my/My;
196: if (ratioj*My != my) SETERRQ2(1,"Ratio between levels: my/My must be integer: my %d My %d",my,My);
197: } else {
198: ratioj = (my-1)/(My-1);
199: if (ratioj*(My-1) != my-1) SETERRQ2(1,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %d My %d",my,My);
200: }
203: DAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
204: DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
205: DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
207: DAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
208: DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
209: DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
211: /*
212: Used for handling a coarse DA that lives on 1/4 the processors of the fine DA.
213: The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
214: processors). It's effective length is hence 4 times its normal length, this is
215: why the col_scale is multiplied by the interpolation matrix column sizes.
216: sol_shift allows each set of 1/4 processors do its own interpolation using ITS
217: copy of the coarse vector. A bit of a hack but you do better.
219: In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
220: */
221: MPI_Comm_size(dac->comm,&size_c);
222: MPI_Comm_size(daf->comm,&size_f);
223: MPI_Comm_rank(daf->comm,&rank_f);
224: col_scale = size_f/size_c;
225: col_shift = Mx*My*(rank_f/size_c);
227: MatPreallocateInitialize(daf->comm,m_f*n_f,col_scale*m_c*n_c,dnz,onz);
228: for (j=j_start; j<j_start+n_f; j++) {
229: for (i=i_start; i<i_start+m_f; i++) {
230: /* convert to local "natural" numbering and then to PETSc global numbering */
231: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
233: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
234: j_c = (j/ratioj); /* coarse grid node below fine grid node */
236: if (j_c < j_start_ghost_c) SETERRQ3(1,"Processor's coarse DA must lie over fine DAn
237: j_start %d j_c %d j_start_ghost_c %d",j_start,j_c,j_start_ghost_c);
238: if (i_c < i_start_ghost_c) SETERRQ3(1,"Processor's coarse DA must lie over fine DAn
239: i_start %d i_c %d i_start_ghost_c %d",i_start,i_c,i_start_ghost_c);
241: /*
242: Only include those interpolation points that are truly
243: nonzero. Note this is very important for final grid lines
244: in x and y directions; since they have no right/top neighbors
245: */
246: nc = 0;
247: /* one left and below; or we are right on it */
248: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
249: cols[nc++] = col_shift + idx_c[col]/dof;
250: /* one right and below */
251: if (i_c*ratioi != i) {
252: cols[nc++] = col_shift + idx_c[col+dof]/dof;
253: }
254: /* one left and above */
255: if (j_c*ratioj != j) {
256: cols[nc++] = col_shift + idx_c[col+m_ghost_c*dof]/dof;
257: }
258: /* one right and above */
259: if (j_c*ratioi != j && i_c*ratioj != i) {
260: cols[nc++] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof;
261: }
262: MatPreallocateSet(row,nc,cols,dnz,onz);
263: }
264: }
265: MatCreateMPIAIJ(daf->comm,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My,0,dnz,0,onz,&mat);
266: MatPreallocateFinalize(dnz,onz);
267: if (!DAXPeriodic(pt) && !DAYPeriodic(pt)) {MatSetOption(mat,MAT_COLUMNS_SORTED);}
269: /* loop over local fine grid nodes setting interpolation for those*/
270: for (j=j_start; j<j_start+n_f; j++) {
271: for (i=i_start; i<i_start+m_f; i++) {
272: /* convert to local "natural" numbering and then to PETSc global numbering */
273: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
275: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
276: j_c = (j/ratioj); /* coarse grid node below fine grid node */
278: /*
279: Only include those interpolation points that are truly
280: nonzero. Note this is very important for final grid lines
281: in x and y directions; since they have no right/top neighbors
282: */
283: x = ((double)(i - i_c*ratioi))/((double)ratioi);
284: y = ((double)(j - j_c*ratioj))/((double)ratioj);
285: /* printf("i j %d %d %g %gn",i,j,x,y); */
286: nc = 0;
287: /* one left and below; or we are right on it */
288: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
289: cols[nc] = col_shift + idx_c[col]/dof;
290: v[nc++] = x*y - x - y + 1.0;
291: /* one right and below */
292: if (i_c*ratioi != i) {
293: cols[nc] = col_shift + idx_c[col+dof]/dof;
294: v[nc++] = -x*y + x;
295: }
296: /* one left and above */
297: if (j_c*ratioj != j) {
298: cols[nc] = col_shift + idx_c[col+m_ghost_c*dof]/dof;
299: v[nc++] = -x*y + y;
300: }
301: /* one right and above */
302: if (j_c*ratioj != j && i_c*ratioi != i) {
303: cols[nc] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof;
304: v[nc++] = x*y;
305: }
306: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
307: }
308: }
309: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
310: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
311: MatCreateMAIJ(mat,dof,A);
312: MatDestroy(mat);
313: PetscLogFlops(13*m_f*n_f);
314: return(0);
315: }
318: /* dof degree of freedom per node, nonperiodic */
319: #undef __FUNCT__
321: int DAGetInterpolation_3D_Q1(DA dac,DA daf,Mat *A)
322: {
323: int ierr,i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof,l;
324: int m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,Mz,mz;
325: int row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok;
326: int i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
327: int l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c;
328: int l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz;
329: PetscScalar v[8],x,y,z;
330: Mat mat;
331: DAPeriodicType pt;
334: DAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&pt,0);
335: DAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);
336: if (DAXPeriodic(pt)){
337: ratioi = mx/Mx;
338: if (ratioi*Mx != mx) SETERRQ2(1,"Ratio between levels: mx/Mx must be integer: mx %d Mx %d",mx,Mx);
339: } else {
340: ratioi = (mx-1)/(Mx-1);
341: if (ratioi*(Mx-1) != mx-1) SETERRQ2(1,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %d Mx %d",mx,Mx);
342: }
343: if (DAYPeriodic(pt)){
344: ratioj = my/My;
345: if (ratioj*My != my) SETERRQ2(1,"Ratio between levels: my/My must be integer: my %d My %d",my,My);
346: } else {
347: ratioj = (my-1)/(My-1);
348: if (ratioj*(My-1) != my-1) SETERRQ2(1,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %d My %d",my,My);
349: }
350: if (DAZPeriodic(pt)){
351: ratiok = mz/Mz;
352: if (ratiok*Mz != mz) SETERRQ2(1,"Ratio between levels: mz/Mz must be integer: mz %d Mz %d",mz,Mz);
353: } else {
354: ratiok = (mz-1)/(Mz-1);
355: if (ratiok*(Mz-1) != mz-1) SETERRQ2(1,"Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %d Mz %d",mz,Mz);
356: }
358: DAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
359: DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
360: DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);
362: DAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
363: DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
364: DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);
366: /* create interpolation matrix, determining exact preallocation */
367: MatPreallocateInitialize(dac->comm,m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);
368: /* loop over local fine grid nodes counting interpolating points */
369: for (l=l_start; l<l_start+p_f; l++) {
370: for (j=j_start; j<j_start+n_f; j++) {
371: for (i=i_start; i<i_start+m_f; i++) {
372: /* convert to local "natural" numbering and then to PETSc global numbering */
373: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
374: i_c = (i/ratioi);
375: j_c = (j/ratioj);
376: l_c = (l/ratiok);
377: if (l_c < l_start_ghost_c) SETERRQ3(1,"Processor's coarse DA must lie over fine DAn
378: l_start %d l_c %d l_start_ghost_c %d",l_start,l_c,l_start_ghost_c);
379: if (j_c < j_start_ghost_c) SETERRQ3(1,"Processor's coarse DA must lie over fine DAn
380: j_start %d j_c %d j_start_ghost_c %d",j_start,j_c,j_start_ghost_c);
381: if (i_c < i_start_ghost_c) SETERRQ3(1,"Processor's coarse DA must lie over fine DAn
382: i_start %d i_c %d i_start_ghost_c %d",i_start,i_c,i_start_ghost_c);
384: /*
385: Only include those interpolation points that are truly
386: nonzero. Note this is very important for final grid lines
387: in x and y directions; since they have no right/top neighbors
388: */
389: nc = 0;
390: 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));
391: cols[nc++] = idx_c[col]/dof;
392: if (i_c*ratioi != i) {
393: cols[nc++] = idx_c[col+dof]/dof;
394: }
395: if (j_c*ratioj != j) {
396: cols[nc++] = idx_c[col+m_ghost_c*dof]/dof;
397: }
398: if (l_c*ratiok != l) {
399: cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof;
400: }
401: if (j_c*ratioj != j && i_c*ratioi != i) {
402: cols[nc++] = idx_c[col+(m_ghost_c+1)*dof]/dof;
403: }
404: if (j_c*ratioj != j && l_c*ratiok != l) {
405: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;
406: }
407: if (i_c*ratioi != i && l_c*ratiok != l) {
408: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof;
409: }
410: if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
411: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof;
412: }
413: MatPreallocateSet(row,nc,cols,dnz,onz);
414: }
415: }
416: }
417: MatCreateMPIAIJ(dac->comm,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz,0,dnz,0,onz,&mat);
418: MatPreallocateFinalize(dnz,onz);
419: if (!DAXPeriodic(pt) && !DAYPeriodic(pt) && !DAZPeriodic(pt)) {MatSetOption(mat,MAT_COLUMNS_SORTED);}
421: /* loop over local fine grid nodes setting interpolation for those*/
422: for (l=l_start; l<l_start+p_f; l++) {
423: for (j=j_start; j<j_start+n_f; j++) {
424: for (i=i_start; i<i_start+m_f; i++) {
425: /* convert to local "natural" numbering and then to PETSc global numbering */
426: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
428: i_c = (i/ratioi);
429: j_c = (j/ratioj);
430: l_c = (l/ratiok);
432: /*
433: Only include those interpolation points that are truly
434: nonzero. Note this is very important for final grid lines
435: in x and y directions; since they have no right/top neighbors
436: */
437: x = ((double)(i - i_c*ratioi))/((double)ratioi);
438: y = ((double)(j - j_c*ratioj))/((double)ratioj);
439: z = ((double)(l - l_c*ratiok))/((double)ratiok);
440: /* printf("i j l %d %d %d %g %g %gn",i,j,l,x,y,z); */
441: nc = 0;
442: /* one left and below; or we are right on it */
443: 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));
445: cols[nc] = idx_c[col]/dof;
446: v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
448: if (i_c*ratioi != i) {
449: cols[nc] = idx_c[col+dof]/dof;
450: v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
451: }
453: if (j_c*ratioj != j) {
454: cols[nc] = idx_c[col+m_ghost_c*dof]/dof;
455: v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
456: }
458: if (l_c*ratiok != l) {
459: cols[nc] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof;
460: v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
461: }
463: if (j_c*ratioj != j && i_c*ratioi != i) {
464: cols[nc] = idx_c[col+(m_ghost_c+1)*dof]/dof;
465: v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
466: }
468: if (j_c*ratioj != j && l_c*ratiok != l) {
469: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;
470: v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
471: }
473: if (i_c*ratioi != i && l_c*ratiok != l) {
474: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof;
475: v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
476: }
478: if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
479: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof;
480: v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
481: }
482: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
483: }
484: }
485: }
486: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
487: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
489: MatCreateMAIJ(mat,dof,A);
490: MatDestroy(mat);
491: PetscLogFlops(13*m_f*n_f);
492: return(0);
493: }
495: #undef __FUNCT__
497: /*@C
498: DAGetInterpolation - Gets an interpolation matrix that maps between
499: grids associated with two DAs.
501: Collective on DA
503: Input Parameters:
504: + dac - the coarse grid DA
505: - daf - the fine grid DA
507: Output Parameters:
508: + A - the interpolation matrix
509: - scale - a scaling vector used to scale the coarse grid restricted vector before applying the
510: grid function or grid Jacobian to it.
512: Level: intermediate
514: .keywords: interpolation, restriction, multigrid
516: .seealso: DARefine()
517: @*/
518: int DAGetInterpolation(DA dac,DA daf,Mat *A,Vec *scale)
519: {
520: int ierr,dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
521: DAPeriodicType wrapc,wrapf;
522: DAStencilType stc,stf;
528: DAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&wrapc,&stc);
529: DAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&wrapf,&stf);
530: if (dimc != dimf) SETERRQ2(1,"Dimensions of DA do not match %d %d",dimc,dimf);
531: /* if (mc != mf) SETERRQ2(1,"Processor dimensions of DA in X %d %d do not match",mc,mf);
532: if (nc != nf) SETERRQ2(1,"Processor dimensions of DA in Y %d %d do not match",nc,nf);
533: if (pc != pf) SETERRQ2(1,"Processor dimensions of DA in Z %d %d do not match",pc,pf); */
534: if (dofc != doff) SETERRQ2(1,"DOF of DA do not match %d %d",dofc,doff);
535: if (sc != sf) SETERRQ2(1,"Stencil width of DA do not match %d %d",sc,sf);
536: if (wrapc != wrapf) SETERRQ(1,"Periodic type different in two DAs");
537: if (stc != stf) SETERRQ(1,"Stencil type different in two DAs");
538: if (Mc < 2) SETERRQ(1,"Coarse grid requires at least 2 points in x direction");
539: if (dimc > 1 && Nc < 2) SETERRQ(1,"Coarse grid requires at least 2 points in y direction");
540: if (dimc > 2 && Pc < 2) SETERRQ(1,"Coarse grid requires at least 2 points in z direction");
542: if (dac->interptype == DA_Q1){
543: if (dimc == 1){
544: DAGetInterpolation_1D_Q1(dac,daf,A);
545: } else if (dimc == 2){
546: DAGetInterpolation_2D_Q1(dac,daf,A);
547: } else if (dimc == 3){
548: DAGetInterpolation_3D_Q1(dac,daf,A);
549: } else {
550: SETERRQ2(1,"No support for this DA dimension %d for interpolation type %d",dimc,dac->interptype);
551: }
552: } else if (dac->interptype == DA_Q0){
553: if (dimc == 1){
554: DAGetInterpolation_1D_Q0(dac,daf,A);
555: } else {
556: SETERRQ2(1,"No support for this DA dimension %d for interpolation type %d",dimc,dac->interptype);
557: }
558: }
559: if (scale) {
560: DMGetInterpolationScale((DM)dac,(DM)daf,*A,scale);
561: }
562: return(0);
563: }