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> /*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;
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: }
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,NULL,&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: DMDAGetGlobalIndices(dac,NULL,&idx_c);
83: /* create interpolation matrix */
84: MatCreate(PetscObjectComm((PetscObject)dac),&mat);
85: MatSetSizes(mat,m_f,m_c,mx,Mx);
86: MatSetType(mat,MATAIJ);
87: MatSeqAIJSetPreallocation(mat,2,NULL);
88: MatMPIAIJSetPreallocation(mat,2,NULL,1,NULL);
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;
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);
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: }
120: } else {
121: PetscScalar *xi;
122: PetscInt li,nxi,n;
123: PetscScalar Ni[2];
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;
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;
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: }
153: /* edges + interior */
154: /* remainders */
155: if (i==mx-1) i_c--;
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;
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,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,NULL,&idx_c);
209: /* create interpolation matrix */
210: MatCreate(PetscObjectComm((PetscObject)dac),&mat);
211: MatSetSizes(mat,m_f,m_c,mx,Mx);
212: MatSetType(mat,MATAIJ);
213: MatSeqAIJSetPreallocation(mat,2,NULL);
214: MatMPIAIJSetPreallocation(mat,2,NULL,0,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;
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: }
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,NULL,&idx_f);
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,NULL,&idx_c);
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.
298: In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
299: */
300: MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);
301: MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);
302: MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);
303: col_scale = size_f/size_c;
304: col_shift = Mx*My*(rank_f/size_c);
306: MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),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;
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 */
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);
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) cols[nc++] = col_shift + idx_c[col+dof]/dof;
331: /* one left and above */
332: if (j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+m_ghost_c*dof]/dof;
333: /* one right and above */
334: if (i_c*ratioi != i && j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof;
335: MatPreallocateSet(row,nc,cols,dnz,onz);
336: }
337: }
338: MatCreate(PetscObjectComm((PetscObject)daf),&mat);
339: MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);
340: MatSetType(mat,MATAIJ);
341: MatSeqAIJSetPreallocation(mat,0,dnz);
342: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
343: MatPreallocateFinalize(dnz,onz);
345: /* loop over local fine grid nodes setting interpolation for those*/
346: if (!NEWVERSION) {
348: for (j=j_start; j<j_start+n_f; j++) {
349: for (i=i_start; i<i_start+m_f; i++) {
350: /* convert to local "natural" numbering and then to PETSc global numbering */
351: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
353: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
354: j_c = (j/ratioj); /* coarse grid node below fine grid node */
356: /*
357: Only include those interpolation points that are truly
358: nonzero. Note this is very important for final grid lines
359: in x and y directions; since they have no right/top neighbors
360: */
361: x = ((double)(i - i_c*ratioi))/((double)ratioi);
362: y = ((double)(j - j_c*ratioj))/((double)ratioj);
364: nc = 0;
365: /* one left and below; or we are right on it */
366: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
367: cols[nc] = col_shift + idx_c[col]/dof;
368: v[nc++] = x*y - x - y + 1.0;
369: /* one right and below */
370: if (i_c*ratioi != i) {
371: cols[nc] = col_shift + idx_c[col+dof]/dof;
372: v[nc++] = -x*y + x;
373: }
374: /* one left and above */
375: if (j_c*ratioj != j) {
376: cols[nc] = col_shift + idx_c[col+m_ghost_c*dof]/dof;
377: v[nc++] = -x*y + y;
378: }
379: /* one right and above */
380: if (j_c*ratioj != j && i_c*ratioi != i) {
381: cols[nc] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof;
382: v[nc++] = x*y;
383: }
384: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
385: }
386: }
388: } else {
389: PetscScalar Ni[4];
390: PetscScalar *xi,*eta;
391: PetscInt li,nxi,lj,neta;
393: /* compute local coordinate arrays */
394: nxi = ratioi + 1;
395: neta = ratioj + 1;
396: PetscMalloc(sizeof(PetscScalar)*nxi,&xi);
397: PetscMalloc(sizeof(PetscScalar)*neta,&eta);
398: for (li=0; li<nxi; li++) {
399: xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
400: }
401: for (lj=0; lj<neta; lj++) {
402: eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
403: }
405: /* loop over local fine grid nodes setting interpolation for those*/
406: for (j=j_start; j<j_start+n_f; j++) {
407: for (i=i_start; i<i_start+m_f; i++) {
408: /* convert to local "natural" numbering and then to PETSc global numbering */
409: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
411: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
412: j_c = (j/ratioj); /* coarse grid node below fine grid node */
414: /* remainders */
415: li = i - ratioi * (i/ratioi);
416: if (i==mx-1) li = nxi-1;
417: lj = j - ratioj * (j/ratioj);
418: if (j==my-1) lj = neta-1;
420: /* corners */
421: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
422: cols[0] = col_shift + idx_c[col]/dof; /* left, below */
423: Ni[0] = 1.0;
424: if ((li==0) || (li==nxi-1)) {
425: if ((lj==0) || (lj==neta-1)) {
426: MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
427: continue;
428: }
429: }
431: /* edges + interior */
432: /* remainders */
433: if (i==mx-1) i_c--;
434: if (j==my-1) j_c--;
436: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
437: cols[0] = col_shift + idx_c[col]/dof; /* left, below */
438: cols[1] = col_shift + idx_c[col+dof]/dof; /* right, below */
439: cols[2] = col_shift + idx_c[col+m_ghost_c*dof]/dof; /* left, above */
440: cols[3] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof; /* right, above */
442: Ni[0] = 0.25*(1.0-xi[li])*(1.0-eta[lj]);
443: Ni[1] = 0.25*(1.0+xi[li])*(1.0-eta[lj]);
444: Ni[2] = 0.25*(1.0-xi[li])*(1.0+eta[lj]);
445: Ni[3] = 0.25*(1.0+xi[li])*(1.0+eta[lj]);
447: nc = 0;
448: if (PetscAbsScalar(Ni[0])<1.0e-32) cols[0]=-1;
449: if (PetscAbsScalar(Ni[1])<1.0e-32) cols[1]=-1;
450: if (PetscAbsScalar(Ni[2])<1.0e-32) cols[2]=-1;
451: if (PetscAbsScalar(Ni[3])<1.0e-32) cols[3]=-1;
453: MatSetValues(mat,1,&row,4,cols,Ni,INSERT_VALUES);
454: }
455: }
456: PetscFree(xi);
457: PetscFree(eta);
458: }
459: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
460: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
461: MatCreateMAIJ(mat,dof,A);
462: MatDestroy(&mat);
463: return(0);
464: }
466: /*
467: Contributed by Andrei Draganescu <aidraga@sandia.gov>
468: */
471: PetscErrorCode DMCreateInterpolation_DA_2D_Q0(DM dac,DM daf,Mat *A)472: {
473: PetscErrorCode ierr;
474: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof;
475: PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,*dnz,*onz;
476: PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
477: 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;
478: PetscMPIInt size_c,size_f,rank_f;
479: PetscScalar v[4];
480: Mat mat;
481: DMDABoundaryType bx,by;
484: DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
485: DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
486: if (bx == DMDA_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in x");
487: if (by == DMDA_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in y");
488: ratioi = mx/Mx;
489: ratioj = my/My;
490: if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x");
491: if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y");
492: if (ratioi != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 2");
493: if (ratioj != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 2");
495: DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
496: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
497: DMDAGetGlobalIndices(daf,NULL,&idx_f);
499: DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
500: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
501: DMDAGetGlobalIndices(dac,NULL,&idx_c);
503: /*
504: Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
505: The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
506: processors). It's effective length is hence 4 times its normal length, this is
507: why the col_scale is multiplied by the interpolation matrix column sizes.
508: sol_shift allows each set of 1/4 processors do its own interpolation using ITS
509: copy of the coarse vector. A bit of a hack but you do better.
511: In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
512: */
513: MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);
514: MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);
515: MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);
516: col_scale = size_f/size_c;
517: col_shift = Mx*My*(rank_f/size_c);
519: MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f,col_scale*m_c*n_c,dnz,onz);
520: for (j=j_start; j<j_start+n_f; j++) {
521: for (i=i_start; i<i_start+m_f; i++) {
522: /* convert to local "natural" numbering and then to PETSc global numbering */
523: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
525: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
526: j_c = (j/ratioj); /* coarse grid node below fine grid node */
528: 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\
529: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
530: 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\
531: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
533: /*
534: Only include those interpolation points that are truly
535: nonzero. Note this is very important for final grid lines
536: in x and y directions; since they have no right/top neighbors
537: */
538: nc = 0;
539: /* one left and below; or we are right on it */
540: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
541: cols[nc++] = col_shift + idx_c[col]/dof;
542: MatPreallocateSet(row,nc,cols,dnz,onz);
543: }
544: }
545: MatCreate(PetscObjectComm((PetscObject)daf),&mat);
546: MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);
547: MatSetType(mat,MATAIJ);
548: MatSeqAIJSetPreallocation(mat,0,dnz);
549: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
550: MatPreallocateFinalize(dnz,onz);
552: /* loop over local fine grid nodes setting interpolation for those*/
553: for (j=j_start; j<j_start+n_f; j++) {
554: for (i=i_start; i<i_start+m_f; i++) {
555: /* convert to local "natural" numbering and then to PETSc global numbering */
556: row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
558: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
559: j_c = (j/ratioj); /* coarse grid node below fine grid node */
560: nc = 0;
561: /* one left and below; or we are right on it */
562: col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
563: cols[nc] = col_shift + idx_c[col]/dof;
564: v[nc++] = 1.0;
566: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
567: }
568: }
569: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
570: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
571: MatCreateMAIJ(mat,dof,A);
572: MatDestroy(&mat);
573: PetscLogFlops(13.0*m_f*n_f);
574: return(0);
575: }
577: /*
578: Contributed by Jianming Yang <jianming-yang@uiowa.edu>
579: */
582: PetscErrorCode DMCreateInterpolation_DA_3D_Q0(DM dac,DM daf,Mat *A)583: {
584: PetscErrorCode ierr;
585: PetscInt i,j,l,i_start,j_start,l_start,m_f,n_f,p_f,Mx,My,Mz,*idx_f,dof;
586: PetscInt m_ghost,n_ghost,p_ghost,*idx_c,m_ghost_c,n_ghost_c,p_ghost_c,nc,*dnz,*onz;
587: 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;
588: 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;
589: PetscMPIInt size_c,size_f,rank_f;
590: PetscScalar v[8];
591: Mat mat;
592: DMDABoundaryType bx,by,bz;
595: DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
596: DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
597: if (bx == DMDA_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in x");
598: if (by == DMDA_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in y");
599: if (bz == DMDA_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in z");
600: ratioi = mx/Mx;
601: ratioj = my/My;
602: ratiol = mz/Mz;
603: if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x");
604: if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y");
605: if (ratiol*Mz != mz) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in z");
606: if (ratioi != 2 && ratioi != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 1 or 2");
607: if (ratioj != 2 && ratioj != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 1 or 2");
608: if (ratiol != 2 && ratiol != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in z must be 1 or 2");
610: DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
611: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
612: DMDAGetGlobalIndices(daf,NULL,&idx_f);
614: DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
615: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
616: DMDAGetGlobalIndices(dac,NULL,&idx_c);
617: /*
618: Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
619: The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
620: processors). It's effective length is hence 4 times its normal length, this is
621: why the col_scale is multiplied by the interpolation matrix column sizes.
622: sol_shift allows each set of 1/4 processors do its own interpolation using ITS
623: copy of the coarse vector. A bit of a hack but you do better.
625: In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
626: */
627: MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);
628: MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);
629: MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);
630: col_scale = size_f/size_c;
631: col_shift = Mx*My*Mz*(rank_f/size_c);
633: MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f*p_f,col_scale*m_c*n_c*p_c,dnz,onz);
634: for (l=l_start; l<l_start+p_f; l++) {
635: for (j=j_start; j<j_start+n_f; j++) {
636: for (i=i_start; i<i_start+m_f; i++) {
637: /* convert to local "natural" numbering and then to PETSc global numbering */
638: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
640: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
641: j_c = (j/ratioj); /* coarse grid node below fine grid node */
642: l_c = (l/ratiol);
644: 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\
645: l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
646: 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\
647: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
648: 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\
649: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
651: /*
652: Only include those interpolation points that are truly
653: nonzero. Note this is very important for final grid lines
654: in x and y directions; since they have no right/top neighbors
655: */
656: nc = 0;
657: /* one left and below; or we are right on it */
658: 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));
659: cols[nc++] = col_shift + idx_c[col]/dof;
660: MatPreallocateSet(row,nc,cols,dnz,onz);
661: }
662: }
663: }
664: MatCreate(PetscObjectComm((PetscObject)daf),&mat);
665: MatSetSizes(mat,m_f*n_f*p_f,col_scale*m_c*n_c*p_c,mx*my*mz,col_scale*Mx*My*Mz);
666: MatSetType(mat,MATAIJ);
667: MatSeqAIJSetPreallocation(mat,0,dnz);
668: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
669: MatPreallocateFinalize(dnz,onz);
671: /* loop over local fine grid nodes setting interpolation for those*/
672: for (l=l_start; l<l_start+p_f; l++) {
673: for (j=j_start; j<j_start+n_f; j++) {
674: for (i=i_start; i<i_start+m_f; i++) {
675: /* convert to local "natural" numbering and then to PETSc global numbering */
676: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
678: i_c = (i/ratioi); /* coarse grid node to left of fine grid node */
679: j_c = (j/ratioj); /* coarse grid node below fine grid node */
680: l_c = (l/ratiol);
681: nc = 0;
682: /* one left and below; or we are right on it */
683: 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));
684: cols[nc] = col_shift + idx_c[col]/dof;
685: v[nc++] = 1.0;
687: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
688: }
689: }
690: }
691: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
692: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
693: MatCreateMAIJ(mat,dof,A);
694: MatDestroy(&mat);
695: PetscLogFlops(13.0*m_f*n_f*p_f);
696: return(0);
697: }
701: PetscErrorCode DMCreateInterpolation_DA_3D_Q1(DM dac,DM daf,Mat *A)702: {
703: PetscErrorCode ierr;
704: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof,l;
705: PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,Mz,mz;
706: PetscInt row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok;
707: PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
708: PetscInt l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c;
709: PetscInt l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz;
710: PetscScalar v[8],x,y,z;
711: Mat mat;
712: DMDABoundaryType bx,by,bz;
715: DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
716: DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
717: if (mx == Mx) {
718: ratioi = 1;
719: } else if (bx == DMDA_BOUNDARY_PERIODIC) {
720: ratioi = mx/Mx;
721: 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);
722: } else {
723: ratioi = (mx-1)/(Mx-1);
724: 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);
725: }
726: if (my == My) {
727: ratioj = 1;
728: } else if (by == DMDA_BOUNDARY_PERIODIC) {
729: ratioj = my/My;
730: 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);
731: } else {
732: ratioj = (my-1)/(My-1);
733: 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);
734: }
735: if (mz == Mz) {
736: ratiok = 1;
737: } else if (bz == DMDA_BOUNDARY_PERIODIC) {
738: ratiok = mz/Mz;
739: 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);
740: } else {
741: ratiok = (mz-1)/(Mz-1);
742: 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);
743: }
745: DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
746: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
747: DMDAGetGlobalIndices(daf,NULL,&idx_f);
749: DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
750: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
751: DMDAGetGlobalIndices(dac,NULL,&idx_c);
753: /* create interpolation matrix, determining exact preallocation */
754: MatPreallocateInitialize(PetscObjectComm((PetscObject)dac),m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);
755: /* loop over local fine grid nodes counting interpolating points */
756: for (l=l_start; l<l_start+p_f; l++) {
757: for (j=j_start; j<j_start+n_f; j++) {
758: for (i=i_start; i<i_start+m_f; i++) {
759: /* convert to local "natural" numbering and then to PETSc global numbering */
760: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
761: i_c = (i/ratioi);
762: j_c = (j/ratioj);
763: l_c = (l/ratiok);
764: 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\
765: l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
766: 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\
767: j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
768: 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\
769: i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);
771: /*
772: Only include those interpolation points that are truly
773: nonzero. Note this is very important for final grid lines
774: in x and y directions; since they have no right/top neighbors
775: */
776: nc = 0;
777: 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));
778: cols[nc++] = idx_c[col]/dof;
779: if (i_c*ratioi != i) {
780: cols[nc++] = idx_c[col+dof]/dof;
781: }
782: if (j_c*ratioj != j) {
783: cols[nc++] = idx_c[col+m_ghost_c*dof]/dof;
784: }
785: if (l_c*ratiok != l) {
786: cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof;
787: }
788: if (j_c*ratioj != j && i_c*ratioi != i) {
789: cols[nc++] = idx_c[col+(m_ghost_c+1)*dof]/dof;
790: }
791: if (j_c*ratioj != j && l_c*ratiok != l) {
792: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;
793: }
794: if (i_c*ratioi != i && l_c*ratiok != l) {
795: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof;
796: }
797: if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
798: cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof;
799: }
800: MatPreallocateSet(row,nc,cols,dnz,onz);
801: }
802: }
803: }
804: MatCreate(PetscObjectComm((PetscObject)dac),&mat);
805: MatSetSizes(mat,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz);
806: MatSetType(mat,MATAIJ);
807: MatSeqAIJSetPreallocation(mat,0,dnz);
808: MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
809: MatPreallocateFinalize(dnz,onz);
811: /* loop over local fine grid nodes setting interpolation for those*/
812: if (!NEWVERSION) {
814: for (l=l_start; l<l_start+p_f; l++) {
815: for (j=j_start; j<j_start+n_f; j++) {
816: for (i=i_start; i<i_start+m_f; i++) {
817: /* convert to local "natural" numbering and then to PETSc global numbering */
818: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
820: i_c = (i/ratioi);
821: j_c = (j/ratioj);
822: l_c = (l/ratiok);
824: /*
825: Only include those interpolation points that are truly
826: nonzero. Note this is very important for final grid lines
827: in x and y directions; since they have no right/top neighbors
828: */
829: x = ((double)(i - i_c*ratioi))/((double)ratioi);
830: y = ((double)(j - j_c*ratioj))/((double)ratioj);
831: z = ((double)(l - l_c*ratiok))/((double)ratiok);
833: nc = 0;
834: /* one left and below; or we are right on it */
835: 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));
837: cols[nc] = idx_c[col]/dof;
838: v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
840: if (i_c*ratioi != i) {
841: cols[nc] = idx_c[col+dof]/dof;
842: v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
843: }
845: if (j_c*ratioj != j) {
846: cols[nc] = idx_c[col+m_ghost_c*dof]/dof;
847: v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
848: }
850: if (l_c*ratiok != l) {
851: cols[nc] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof;
852: v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
853: }
855: if (j_c*ratioj != j && i_c*ratioi != i) {
856: cols[nc] = idx_c[col+(m_ghost_c+1)*dof]/dof;
857: v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
858: }
860: if (j_c*ratioj != j && l_c*ratiok != l) {
861: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;
862: v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
863: }
865: if (i_c*ratioi != i && l_c*ratiok != l) {
866: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof;
867: v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
868: }
870: if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
871: cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof;
872: v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
873: }
874: MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
875: }
876: }
877: }
879: } else {
880: PetscScalar *xi,*eta,*zeta;
881: PetscInt li,nxi,lj,neta,lk,nzeta,n;
882: PetscScalar Ni[8];
884: /* compute local coordinate arrays */
885: nxi = ratioi + 1;
886: neta = ratioj + 1;
887: nzeta = ratiok + 1;
888: PetscMalloc(sizeof(PetscScalar)*nxi,&xi);
889: PetscMalloc(sizeof(PetscScalar)*neta,&eta);
890: PetscMalloc(sizeof(PetscScalar)*nzeta,&zeta);
891: for (li=0; li<nxi; li++) xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
892: for (lj=0; lj<neta; lj++) eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
893: for (lk=0; lk<nzeta; lk++) zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));
895: for (l=l_start; l<l_start+p_f; l++) {
896: for (j=j_start; j<j_start+n_f; j++) {
897: for (i=i_start; i<i_start+m_f; i++) {
898: /* convert to local "natural" numbering and then to PETSc global numbering */
899: row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
901: i_c = (i/ratioi);
902: j_c = (j/ratioj);
903: l_c = (l/ratiok);
905: /* remainders */
906: li = i - ratioi * (i/ratioi);
907: if (i==mx-1) li = nxi-1;
908: lj = j - ratioj * (j/ratioj);
909: if (j==my-1) lj = neta-1;
910: lk = l - ratiok * (l/ratiok);
911: if (l==mz-1) lk = nzeta-1;
913: /* corners */
914: 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));
915: cols[0] = idx_c[col]/dof;
916: Ni[0] = 1.0;
917: if ((li==0) || (li==nxi-1)) {
918: if ((lj==0) || (lj==neta-1)) {
919: if ((lk==0) || (lk==nzeta-1)) {
920: MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
921: continue;
922: }
923: }
924: }
926: /* edges + interior */
927: /* remainders */
928: if (i==mx-1) i_c--;
929: if (j==my-1) j_c--;
930: if (l==mz-1) l_c--;
932: 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));
933: cols[0] = idx_c[col]/dof; /* one left and below; or we are right on it */
934: cols[1] = idx_c[col+dof]/dof; /* one right and below */
935: cols[2] = idx_c[col+m_ghost_c*dof]/dof; /* one left and above */
936: cols[3] = idx_c[col+(m_ghost_c+1)*dof]/dof; /* one right and above */
938: 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 */
939: cols[5] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof; /* one right and below, and front */
940: cols[6] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof; /* one left and above and front*/
941: cols[7] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof; /* one right and above and front */
943: Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
944: Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
945: Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
946: Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
948: Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
949: Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
950: Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
951: Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
953: for (n=0; n<8; n++) {
954: if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1;
955: }
956: MatSetValues(mat,1,&row,8,cols,Ni,INSERT_VALUES);
958: }
959: }
960: }
961: PetscFree(xi);
962: PetscFree(eta);
963: PetscFree(zeta);
964: }
966: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
967: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
969: MatCreateMAIJ(mat,dof,A);
970: MatDestroy(&mat);
971: return(0);
972: }
976: PetscErrorCode DMCreateInterpolation_DA(DM dac,DM daf,Mat *A,Vec *scale)977: {
978: PetscErrorCode ierr;
979: PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
980: DMDABoundaryType bxc,byc,bzc,bxf,byf,bzf;
981: DMDAStencilType stc,stf;
982: DM_DA *ddc = (DM_DA*)dac->data;
990: DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
991: DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
992: if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
993: if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
994: if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
995: if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
996: if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
997: if (Mc < 2 && Mf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
998: 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");
999: 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");
1001: if (ddc->interptype == DMDA_Q1) {
1002: if (dimc == 1) {
1003: DMCreateInterpolation_DA_1D_Q1(dac,daf,A);
1004: } else if (dimc == 2) {
1005: DMCreateInterpolation_DA_2D_Q1(dac,daf,A);
1006: } else if (dimc == 3) {
1007: DMCreateInterpolation_DA_3D_Q1(dac,daf,A);
1008: } else SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype);
1009: } else if (ddc->interptype == DMDA_Q0) {
1010: if (dimc == 1) {
1011: DMCreateInterpolation_DA_1D_Q0(dac,daf,A);
1012: } else if (dimc == 2) {
1013: DMCreateInterpolation_DA_2D_Q0(dac,daf,A);
1014: } else if (dimc == 3) {
1015: DMCreateInterpolation_DA_3D_Q0(dac,daf,A);
1016: } else SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype);
1017: }
1018: if (scale) {
1019: DMCreateInterpolationScale((DM)dac,(DM)daf,*A,scale);
1020: }
1021: return(0);
1022: }
1026: PetscErrorCode DMCreateInjection_DA_1D(DM dac,DM daf,VecScatter *inject)1027: {
1028: PetscErrorCode ierr;
1029: PetscInt i,i_start,m_f,Mx,*idx_f,dof;
1030: PetscInt m_ghost,*idx_c,m_ghost_c;
1031: PetscInt row,i_start_ghost,mx,m_c,nc,ratioi;
1032: PetscInt i_start_c,i_start_ghost_c;
1033: PetscInt *cols;
1034: DMDABoundaryType bx;
1035: Vec vecf,vecc;
1036: IS isf;
1039: DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);
1040: DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);
1041: if (bx == DMDA_BOUNDARY_PERIODIC) {
1042: ratioi = mx/Mx;
1043: 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);
1044: } else {
1045: ratioi = (mx-1)/(Mx-1);
1046: 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);
1047: }
1049: DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);
1050: DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
1051: DMDAGetGlobalIndices(daf,NULL,&idx_f);
1053: DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
1054: DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
1055: DMDAGetGlobalIndices(dac,NULL,&idx_c);
1058: /* loop over local fine grid nodes setting interpolation for those*/
1059: nc = 0;
1060: PetscMalloc(m_f*sizeof(PetscInt),&cols);
1063: for (i=i_start_c; i<i_start_c+m_c; i++) {
1064: PetscInt i_f = i*ratioi;
1066: 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);
1068: row = idx_f[dof*(i_f-i_start_ghost)];
1069: cols[nc++] = row/dof;
1070: }
1073: ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);
1074: DMGetGlobalVector(dac,&vecc);
1075: DMGetGlobalVector(daf,&vecf);
1076: VecScatterCreate(vecf,isf,vecc,NULL,inject);
1077: DMRestoreGlobalVector(dac,&vecc);
1078: DMRestoreGlobalVector(daf,&vecf);
1079: ISDestroy(&isf);
1080: return(0);
1081: }
1085: PetscErrorCode DMCreateInjection_DA_2D(DM dac,DM daf,VecScatter *inject)1086: {
1087: PetscErrorCode ierr;
1088: PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof;
1089: PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c;
1090: PetscInt row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj;
1091: PetscInt i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
1092: PetscInt *cols;
1093: DMDABoundaryType bx,by;
1094: Vec vecf,vecc;
1095: IS isf;
1098: DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
1099: DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
1100: if (bx == DMDA_BOUNDARY_PERIODIC) {
1101: ratioi = mx/Mx;
1102: 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);
1103: } else {
1104: ratioi = (mx-1)/(Mx-1);
1105: 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);
1106: }
1107: if (by == DMDA_BOUNDARY_PERIODIC) {
1108: ratioj = my/My;
1109: 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);
1110: } else {
1111: ratioj = (my-1)/(My-1);
1112: 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);
1113: }
1115: DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
1116: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
1117: DMDAGetGlobalIndices(daf,NULL,&idx_f);
1119: DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
1120: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
1121: DMDAGetGlobalIndices(dac,NULL,&idx_c);
1124: /* loop over local fine grid nodes setting interpolation for those*/
1125: nc = 0;
1126: PetscMalloc(n_f*m_f*sizeof(PetscInt),&cols);
1127: for (j=j_start_c; j<j_start_c+n_c; j++) {
1128: for (i=i_start_c; i<i_start_c+m_c; i++) {
1129: PetscInt i_f = i*ratioi,j_f = j*ratioj;
1130: 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\
1131: j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost);
1132: 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\
1133: i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1134: row = idx_f[dof*(m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))];
1135: cols[nc++] = row/dof;
1136: }
1137: }
1139: ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);
1140: DMGetGlobalVector(dac,&vecc);
1141: DMGetGlobalVector(daf,&vecf);
1142: VecScatterCreate(vecf,isf,vecc,NULL,inject);
1143: DMRestoreGlobalVector(dac,&vecc);
1144: DMRestoreGlobalVector(daf,&vecf);
1145: ISDestroy(&isf);
1146: return(0);
1147: }
1151: PetscErrorCode DMCreateInjection_DA_3D(DM dac,DM daf,VecScatter *inject)1152: {
1153: PetscErrorCode ierr;
1154: PetscInt i,j,k,i_start,j_start,k_start,m_f,n_f,p_f,Mx,My,Mz;
1155: PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c;
1156: PetscInt i_start_ghost,j_start_ghost,k_start_ghost;
1157: PetscInt mx,my,mz,ratioi,ratioj,ratiok;
1158: PetscInt i_start_c,j_start_c,k_start_c;
1159: PetscInt m_c,n_c,p_c;
1160: PetscInt i_start_ghost_c,j_start_ghost_c,k_start_ghost_c;
1161: PetscInt row,nc,dof;
1162: PetscInt *idx_c,*idx_f;
1163: PetscInt *cols;
1164: DMDABoundaryType bx,by,bz;
1165: Vec vecf,vecc;
1166: IS isf;
1169: DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
1170: DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
1172: if (bx == DMDA_BOUNDARY_PERIODIC) {
1173: ratioi = mx/Mx;
1174: 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);
1175: } else {
1176: ratioi = (mx-1)/(Mx-1);
1177: 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);
1178: }
1179: if (by == DMDA_BOUNDARY_PERIODIC) {
1180: ratioj = my/My;
1181: 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);
1182: } else {
1183: ratioj = (my-1)/(My-1);
1184: 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);
1185: }
1186: if (bz == DMDA_BOUNDARY_PERIODIC) {
1187: ratiok = mz/Mz;
1188: 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);
1189: } else {
1190: ratiok = (mz-1)/(Mz-1);
1191: 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);
1192: }
1194: DMDAGetCorners(daf,&i_start,&j_start,&k_start,&m_f,&n_f,&p_f);
1195: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&k_start_ghost,&m_ghost,&n_ghost,&p_ghost);
1196: DMDAGetGlobalIndices(daf,NULL,&idx_f);
1198: DMDAGetCorners(dac,&i_start_c,&j_start_c,&k_start_c,&m_c,&n_c,&p_c);
1199: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&k_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
1200: DMDAGetGlobalIndices(dac,NULL,&idx_c);
1203: /* loop over local fine grid nodes setting interpolation for those*/
1204: nc = 0;
1205: PetscMalloc(n_f*m_f*p_f*sizeof(PetscInt),&cols);
1206: for (k=k_start_c; k<k_start_c+p_c; k++) {
1207: for (j=j_start_c; j<j_start_c+n_c; j++) {
1208: for (i=i_start_c; i<i_start_c+m_c; i++) {
1209: PetscInt i_f = i*ratioi,j_f = j*ratioj,k_f = k*ratiok;
1210: 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 "
1211: "k_c %D k_f %D fine ghost range [%D,%D]",k,k_f,k_start_ghost,k_start_ghost+p_ghost);
1212: 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 "
1213: "j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost);
1214: 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 "
1215: "i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1216: 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))];
1217: cols[nc++] = row/dof;
1218: }
1219: }
1220: }
1222: ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);
1223: DMGetGlobalVector(dac,&vecc);
1224: DMGetGlobalVector(daf,&vecf);
1225: VecScatterCreate(vecf,isf,vecc,NULL,inject);
1226: DMRestoreGlobalVector(dac,&vecc);
1227: DMRestoreGlobalVector(daf,&vecf);
1228: ISDestroy(&isf);
1229: return(0);
1230: }
1234: PetscErrorCode DMCreateInjection_DA(DM dac,DM daf,VecScatter *inject)1235: {
1236: PetscErrorCode ierr;
1237: PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1238: DMDABoundaryType bxc,byc,bzc,bxf,byf,bzf;
1239: DMDAStencilType stc,stf;
1246: DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1247: DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1248: if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1249: if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1250: if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1251: if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1252: if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1253: if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
1254: if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction");
1255: if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction");
1257: if (dimc == 1) {
1258: DMCreateInjection_DA_1D(dac,daf,inject);
1259: } else if (dimc == 2) {
1260: DMCreateInjection_DA_2D(dac,daf,inject);
1261: } else if (dimc == 3) {
1262: DMCreateInjection_DA_3D(dac,daf,inject);
1263: }
1264: return(0);
1265: }
1269: PetscErrorCode DMCreateAggregates_DA(DM dac,DM daf,Mat *rest)1270: {
1271: PetscErrorCode ierr;
1272: PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc;
1273: PetscInt dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1274: DMDABoundaryType bxc,byc,bzc,bxf,byf,bzf;
1275: DMDAStencilType stc,stf;
1276: PetscInt i,j,l;
1277: PetscInt i_start,j_start,l_start, m_f,n_f,p_f;
1278: PetscInt i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost;
1279: PetscInt *idx_f;
1280: PetscInt i_c,j_c,l_c;
1281: PetscInt i_start_c,j_start_c,l_start_c, m_c,n_c,p_c;
1282: PetscInt i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c;
1283: PetscInt *idx_c;
1284: PetscInt d;
1285: PetscInt a;
1286: PetscInt max_agg_size;
1287: PetscInt *fine_nodes;
1288: PetscScalar *one_vec;
1289: PetscInt fn_idx;
1296: DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1297: DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1298: if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1299: if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1300: if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1301: if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1302: if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1304: 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);
1305: 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);
1306: 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);
1308: if (Pc < 0) Pc = 1;
1309: if (Pf < 0) Pf = 1;
1310: if (Nc < 0) Nc = 1;
1311: if (Nf < 0) Nf = 1;
1313: DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
1314: DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
1315: DMDAGetGlobalIndices(daf,NULL,&idx_f);
1317: DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
1318: DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
1319: DMDAGetGlobalIndices(dac,NULL,&idx_c);
1321: /*
1322: Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios
1323: for dimension 1 and 2 respectively.
1324: Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1)
1325: and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate.
1326: Each specific dof on the fine grid is mapped to one dof on the coarse grid.
1327: */
1329: max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1);
1331: /* create the matrix that will contain the restriction operator */
1332: 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,
1333: max_agg_size, NULL, max_agg_size, NULL, rest);
1335: /* store nodes in the fine grid here */
1336: PetscMalloc2(max_agg_size,PetscScalar, &one_vec,max_agg_size,PetscInt, &fine_nodes);
1337: for (i=0; i<max_agg_size; i++) one_vec[i] = 1.0;
1339: /* loop over all coarse nodes */
1340: for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) {
1341: for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) {
1342: for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) {
1343: for (d=0; d<dofc; d++) {
1344: /* convert to local "natural" numbering and then to PETSc global numbering */
1345: 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;
1347: fn_idx = 0;
1348: /* Corresponding fine points are all points (i_f, j_f, l_f) such that
1349: i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc
1350: (same for other dimensions)
1351: */
1352: for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) {
1353: for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) {
1354: for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) {
1355: 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;
1356: fn_idx++;
1357: }
1358: }
1359: }
1360: /* add all these points to one aggregate */
1361: MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);
1362: }
1363: }
1364: }
1365: }
1366: PetscFree2(one_vec,fine_nodes);
1367: MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);
1368: MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);
1369: return(0);
1370: }