Actual source code: dainterp.c

petsc-3.10.5 2019-03-28
Report Typos and Errors

  2: /*
  3:   Code for interpolating between grids represented by DMDAs
  4: */

  6: /*
  7:       For linear elements there are two branches of code to compute the interpolation. They should compute the same results but may not. The "new version" does
  8:    not work for periodic domains, the old does. Change NEWVERSION to 1 to compile in the new version. Eventually when we are sure the two produce identical results
  9:    we will remove/merge the new version. Based on current tests, these both produce the same results. We are leaving NEWVERSION for now in the code since some
 10:    consider it cleaner, but old version is turned on since it handles periodic case.
 11: */
 12: #define NEWVERSION 0

 14:  #include <petsc/private/dmdaimpl.h>

 16: /*@
 17:     DMCreateInterpolationScale - Forms L = R*1/diag(R*1) - L.*v is like a coarse grid average of the
 18:       nearby fine grid points.

 20:   Input Parameters:
 21: +      dac - DM that defines a coarse mesh
 22: .      daf - DM that defines a fine mesh
 23: -      mat - the restriction (or interpolation operator) from fine to coarse

 25:   Output Parameter:
 26: .    scale - the scaled vector

 28:   Level: developer

 30: .seealso: DMCreateInterpolation()

 32: @*/
 33: PetscErrorCode  DMCreateInterpolationScale(DM dac,DM daf,Mat mat,Vec *scale)
 34: {
 36:   Vec            fine;
 37:   PetscScalar    one = 1.0;

 40:   DMCreateGlobalVector(daf,&fine);
 41:   DMCreateGlobalVector(dac,scale);
 42:   VecSet(fine,one);
 43:   MatRestrict(mat,fine,*scale);
 44:   VecDestroy(&fine);
 45:   VecReciprocal(*scale);
 46:   return(0);
 47: }

 49: PetscErrorCode DMCreateInterpolation_DA_1D_Q1(DM dac,DM daf,Mat *A)
 50: {
 51:   PetscErrorCode         ierr;
 52:   PetscInt               i,i_start,m_f,Mx;
 53:   const PetscInt         *idx_f,*idx_c;
 54:   PetscInt               m_ghost,m_ghost_c;
 55:   PetscInt               row,col,i_start_ghost,mx,m_c,nc,ratio;
 56:   PetscInt               i_c,i_start_c,i_start_ghost_c,cols[2],dof;
 57:   PetscScalar            v[2],x;
 58:   Mat                    mat;
 59:   DMBoundaryType         bx;
 60:   ISLocalToGlobalMapping ltog_f,ltog_c;


 64:   DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);
 65:   DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);
 66:   if (bx == DM_BOUNDARY_PERIODIC) {
 67:     ratio = mx/Mx;
 68:     if (ratio*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
 69:   } else {
 70:     ratio = (mx-1)/(Mx-1);
 71:     if (ratio*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
 72:   }

 74:   DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);
 75:   DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
 76:   DMGetLocalToGlobalMapping(daf,&ltog_f);
 77:   ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);

 79:   DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
 80:   DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
 81:   DMGetLocalToGlobalMapping(dac,&ltog_c);
 82:   ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);

 84:   /* create interpolation matrix */
 85:   MatCreate(PetscObjectComm((PetscObject)dac),&mat);
 86:   MatSetSizes(mat,m_f,m_c,mx,Mx);
 87:   MatSetType(mat,MATAIJ);
 88:   MatSeqAIJSetPreallocation(mat,2,NULL);
 89:   MatMPIAIJSetPreallocation(mat,2,NULL,1,NULL);

 91:   /* loop over local fine grid nodes setting interpolation for those*/
 92:   if (!NEWVERSION) {

 94:     for (i=i_start; i<i_start+m_f; i++) {
 95:       /* convert to local "natural" numbering and then to PETSc global numbering */
 96:       row = idx_f[i-i_start_ghost];

 98:       i_c = (i/ratio);    /* coarse grid node to left of fine grid node */
 99:       if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
100:                                           i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);

102:       /*
103:        Only include those interpolation points that are truly
104:        nonzero. Note this is very important for final grid lines
105:        in x direction; since they have no right neighbor
106:        */
107:       x  = ((PetscReal)(i - i_c*ratio))/((PetscReal)ratio);
108:       nc = 0;
109:       /* one left and below; or we are right on it */
110:       col      = (i_c-i_start_ghost_c);
111:       cols[nc] = idx_c[col];
112:       v[nc++]  = -x + 1.0;
113:       /* one right? */
114:       if (i_c*ratio != i) {
115:         cols[nc] = idx_c[col+1];
116:         v[nc++]  = x;
117:       }
118:       MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
119:     }

121:   } else {
122:     PetscScalar *xi;
123:     PetscInt    li,nxi,n;
124:     PetscScalar Ni[2];

126:     /* compute local coordinate arrays */
127:     nxi  = ratio + 1;
128:     PetscMalloc1(nxi,&xi);
129:     for (li=0; li<nxi; li++) {
130:       xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
131:     }

133:     for (i=i_start; i<i_start+m_f; i++) {
134:       /* convert to local "natural" numbering and then to PETSc global numbering */
135:       row = idx_f[(i-i_start_ghost)];

137:       i_c = (i/ratio);    /* coarse grid node to left of fine grid node */
138:       if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
139:                                           i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);

141:       /* remainders */
142:       li = i - ratio * (i/ratio);
143:       if (i==mx-1) li = nxi-1;

145:       /* corners */
146:       col     = (i_c-i_start_ghost_c);
147:       cols[0] = idx_c[col];
148:       Ni[0]   = 1.0;
149:       if ((li==0) || (li==nxi-1)) {
150:         MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
151:         continue;
152:       }

154:       /* edges + interior */
155:       /* remainders */
156:       if (i==mx-1) i_c--;

158:       col     = (i_c-i_start_ghost_c);
159:       cols[0] = idx_c[col]; /* one left and below; or we are right on it */
160:       cols[1] = idx_c[col+1];

162:       Ni[0] = 0.5*(1.0-xi[li]);
163:       Ni[1] = 0.5*(1.0+xi[li]);
164:       for (n=0; n<2; n++) {
165:         if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1;
166:       }
167:       MatSetValues(mat,1,&row,2,cols,Ni,INSERT_VALUES);
168:     }
169:     PetscFree(xi);
170:   }
171:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
172:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
173:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
174:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
175:   MatCreateMAIJ(mat,dof,A);
176:   MatDestroy(&mat);
177:   return(0);
178: }

180: PetscErrorCode DMCreateInterpolation_DA_1D_Q0(DM dac,DM daf,Mat *A)
181: {
182:   PetscErrorCode         ierr;
183:   PetscInt               i,i_start,m_f,Mx;
184:   const PetscInt         *idx_f,*idx_c;
185:   ISLocalToGlobalMapping ltog_f,ltog_c;
186:   PetscInt               m_ghost,m_ghost_c;
187:   PetscInt               row,col,i_start_ghost,mx,m_c,nc,ratio;
188:   PetscInt               i_c,i_start_c,i_start_ghost_c,cols[2],dof;
189:   PetscScalar            v[2],x;
190:   Mat                    mat;
191:   DMBoundaryType         bx;

194:   DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);
195:   DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);
196:   if (bx == DM_BOUNDARY_PERIODIC) {
197:     if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
198:     ratio = mx/Mx;
199:     if (ratio*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
200:   } else {
201:     if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx);
202:     ratio = (mx-1)/(Mx-1);
203:     if (ratio*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
204:   }

206:   DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);
207:   DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
208:   DMGetLocalToGlobalMapping(daf,&ltog_f);
209:   ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);

211:   DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
212:   DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
213:   DMGetLocalToGlobalMapping(dac,&ltog_c);
214:   ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);

216:   /* create interpolation matrix */
217:   MatCreate(PetscObjectComm((PetscObject)dac),&mat);
218:   MatSetSizes(mat,m_f,m_c,mx,Mx);
219:   MatSetType(mat,MATAIJ);
220:   MatSeqAIJSetPreallocation(mat,2,NULL);
221:   MatMPIAIJSetPreallocation(mat,2,NULL,0,NULL);

223:   /* loop over local fine grid nodes setting interpolation for those*/
224:   for (i=i_start; i<i_start+m_f; i++) {
225:     /* convert to local "natural" numbering and then to PETSc global numbering */
226:     row = idx_f[(i-i_start_ghost)];

228:     i_c = (i/ratio);    /* coarse grid node to left of fine grid node */

230:     /*
231:          Only include those interpolation points that are truly
232:          nonzero. Note this is very important for final grid lines
233:          in x direction; since they have no right neighbor
234:     */
235:     x  = ((PetscReal)(i - i_c*ratio))/((PetscReal)ratio);
236:     nc = 0;
237:     /* one left and below; or we are right on it */
238:     col      = (i_c-i_start_ghost_c);
239:     cols[nc] = idx_c[col];
240:     v[nc++]  = -x + 1.0;
241:     /* one right? */
242:     if (i_c*ratio != i) {
243:       cols[nc] = idx_c[col+1];
244:       v[nc++]  = x;
245:     }
246:     MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
247:   }
248:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
249:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
250:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
251:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
252:   MatCreateMAIJ(mat,dof,A);
253:   MatDestroy(&mat);
254:   PetscLogFlops(5.0*m_f);
255:   return(0);
256: }

258: PetscErrorCode DMCreateInterpolation_DA_2D_Q1(DM dac,DM daf,Mat *A)
259: {
260:   PetscErrorCode         ierr;
261:   PetscInt               i,j,i_start,j_start,m_f,n_f,Mx,My,dof;
262:   const PetscInt         *idx_c,*idx_f;
263:   ISLocalToGlobalMapping ltog_f,ltog_c;
264:   PetscInt               m_ghost,n_ghost,m_ghost_c,n_ghost_c,*dnz,*onz;
265:   PetscInt               row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
266:   PetscInt               i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c,col_shift,col_scale;
267:   PetscMPIInt            size_c,size_f,rank_f;
268:   PetscScalar            v[4],x,y;
269:   Mat                    mat;
270:   DMBoundaryType         bx,by;

273:   DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
274:   DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
275:   if (bx == DM_BOUNDARY_PERIODIC) {
276:     if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
277:     ratioi = mx/Mx;
278:     if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
279:   } else {
280:     if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx);
281:     ratioi = (mx-1)/(Mx-1);
282:     if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
283:   }
284:   if (by == DM_BOUNDARY_PERIODIC) {
285:     if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
286:     ratioj = my/My;
287:     if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My  must be integer: my %D My %D",my,My);
288:   } else {
289:     if (My < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be greater than 1",My);
290:     ratioj = (my-1)/(My-1);
291:     if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
292:   }

294:   DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
295:   DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
296:   DMGetLocalToGlobalMapping(daf,&ltog_f);
297:   ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);

299:   DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
300:   DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
301:   DMGetLocalToGlobalMapping(dac,&ltog_c);
302:   ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);

304:   /*
305:    Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
306:    The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
307:    processors). It's effective length is hence 4 times its normal length, this is
308:    why the col_scale is multiplied by the interpolation matrix column sizes.
309:    sol_shift allows each set of 1/4 processors do its own interpolation using ITS
310:    copy of the coarse vector. A bit of a hack but you do better.

312:    In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
313:    */
314:   MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);
315:   MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);
316:   MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);
317:   col_scale = size_f/size_c;
318:   col_shift = Mx*My*(rank_f/size_c);

320:   MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f,col_scale*m_c*n_c,dnz,onz);
321:   for (j=j_start; j<j_start+n_f; j++) {
322:     for (i=i_start; i<i_start+m_f; i++) {
323:       /* convert to local "natural" numbering and then to PETSc global numbering */
324:       row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];

326:       i_c = (i/ratioi);    /* coarse grid node to left of fine grid node */
327:       j_c = (j/ratioj);    /* coarse grid node below fine grid node */

329:       if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
330:                                           j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
331:       if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
332:                                           i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);

334:       /*
335:        Only include those interpolation points that are truly
336:        nonzero. Note this is very important for final grid lines
337:        in x and y directions; since they have no right/top neighbors
338:        */
339:       nc = 0;
340:       /* one left and below; or we are right on it */
341:       col        = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
342:       cols[nc++] = col_shift + idx_c[col];
343:       /* one right and below */
344:       if (i_c*ratioi != i) cols[nc++] = col_shift + idx_c[col+1];
345:       /* one left and above */
346:       if (j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+m_ghost_c];
347:       /* one right and above */
348:       if (i_c*ratioi != i && j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+(m_ghost_c+1)];
349:       MatPreallocateSet(row,nc,cols,dnz,onz);
350:     }
351:   }
352:   MatCreate(PetscObjectComm((PetscObject)daf),&mat);
353:   MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);
354:   MatSetType(mat,MATAIJ);
355:   MatSeqAIJSetPreallocation(mat,0,dnz);
356:   MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
357:   MatPreallocateFinalize(dnz,onz);

359:   /* loop over local fine grid nodes setting interpolation for those*/
360:   if (!NEWVERSION) {

362:     for (j=j_start; j<j_start+n_f; j++) {
363:       for (i=i_start; i<i_start+m_f; i++) {
364:         /* convert to local "natural" numbering and then to PETSc global numbering */
365:         row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];

367:         i_c = (i/ratioi);    /* coarse grid node to left of fine grid node */
368:         j_c = (j/ratioj);    /* coarse grid node below fine grid node */

370:         /*
371:          Only include those interpolation points that are truly
372:          nonzero. Note this is very important for final grid lines
373:          in x and y directions; since they have no right/top neighbors
374:          */
375:         x = ((PetscReal)(i - i_c*ratioi))/((PetscReal)ratioi);
376:         y = ((PetscReal)(j - j_c*ratioj))/((PetscReal)ratioj);

378:         nc = 0;
379:         /* one left and below; or we are right on it */
380:         col      = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
381:         cols[nc] = col_shift + idx_c[col];
382:         v[nc++]  = x*y - x - y + 1.0;
383:         /* one right and below */
384:         if (i_c*ratioi != i) {
385:           cols[nc] = col_shift + idx_c[col+1];
386:           v[nc++]  = -x*y + x;
387:         }
388:         /* one left and above */
389:         if (j_c*ratioj != j) {
390:           cols[nc] = col_shift + idx_c[col+m_ghost_c];
391:           v[nc++]  = -x*y + y;
392:         }
393:         /* one right and above */
394:         if (j_c*ratioj != j && i_c*ratioi != i) {
395:           cols[nc] = col_shift + idx_c[col+(m_ghost_c+1)];
396:           v[nc++]  = x*y;
397:         }
398:         MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
399:       }
400:     }

402:   } else {
403:     PetscScalar Ni[4];
404:     PetscScalar *xi,*eta;
405:     PetscInt    li,nxi,lj,neta;

407:     /* compute local coordinate arrays */
408:     nxi  = ratioi + 1;
409:     neta = ratioj + 1;
410:     PetscMalloc1(nxi,&xi);
411:     PetscMalloc1(neta,&eta);
412:     for (li=0; li<nxi; li++) {
413:       xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
414:     }
415:     for (lj=0; lj<neta; lj++) {
416:       eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
417:     }

419:     /* loop over local fine grid nodes setting interpolation for those*/
420:     for (j=j_start; j<j_start+n_f; j++) {
421:       for (i=i_start; i<i_start+m_f; i++) {
422:         /* convert to local "natural" numbering and then to PETSc global numbering */
423:         row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];

425:         i_c = (i/ratioi);    /* coarse grid node to left of fine grid node */
426:         j_c = (j/ratioj);    /* coarse grid node below fine grid node */

428:         /* remainders */
429:         li = i - ratioi * (i/ratioi);
430:         if (i==mx-1) li = nxi-1;
431:         lj = j - ratioj * (j/ratioj);
432:         if (j==my-1) lj = neta-1;

434:         /* corners */
435:         col     = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
436:         cols[0] = col_shift + idx_c[col]; /* left, below */
437:         Ni[0]   = 1.0;
438:         if ((li==0) || (li==nxi-1)) {
439:           if ((lj==0) || (lj==neta-1)) {
440:             MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
441:             continue;
442:           }
443:         }

445:         /* edges + interior */
446:         /* remainders */
447:         if (i==mx-1) i_c--;
448:         if (j==my-1) j_c--;

450:         col     = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
451:         cols[0] = col_shift + idx_c[col]; /* left, below */
452:         cols[1] = col_shift + idx_c[col+1]; /* right, below */
453:         cols[2] = col_shift + idx_c[col+m_ghost_c]; /* left, above */
454:         cols[3] = col_shift + idx_c[col+(m_ghost_c+1)]; /* right, above */

456:         Ni[0] = 0.25*(1.0-xi[li])*(1.0-eta[lj]);
457:         Ni[1] = 0.25*(1.0+xi[li])*(1.0-eta[lj]);
458:         Ni[2] = 0.25*(1.0-xi[li])*(1.0+eta[lj]);
459:         Ni[3] = 0.25*(1.0+xi[li])*(1.0+eta[lj]);

461:         nc = 0;
462:         if (PetscAbsScalar(Ni[0])<1.0e-32) cols[0]=-1;
463:         if (PetscAbsScalar(Ni[1])<1.0e-32) cols[1]=-1;
464:         if (PetscAbsScalar(Ni[2])<1.0e-32) cols[2]=-1;
465:         if (PetscAbsScalar(Ni[3])<1.0e-32) cols[3]=-1;

467:         MatSetValues(mat,1,&row,4,cols,Ni,INSERT_VALUES);
468:       }
469:     }
470:     PetscFree(xi);
471:     PetscFree(eta);
472:   }
473:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
474:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
475:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
476:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
477:   MatCreateMAIJ(mat,dof,A);
478:   MatDestroy(&mat);
479:   return(0);
480: }

482: /*
483:        Contributed by Andrei Draganescu <aidraga@sandia.gov>
484: */
485: PetscErrorCode DMCreateInterpolation_DA_2D_Q0(DM dac,DM daf,Mat *A)
486: {
487:   PetscErrorCode         ierr;
488:   PetscInt               i,j,i_start,j_start,m_f,n_f,Mx,My,dof;
489:   const PetscInt         *idx_c,*idx_f;
490:   ISLocalToGlobalMapping ltog_f,ltog_c;
491:   PetscInt               m_ghost,n_ghost,m_ghost_c,n_ghost_c,*dnz,*onz;
492:   PetscInt               row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
493:   PetscInt               i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c,col_shift,col_scale;
494:   PetscMPIInt            size_c,size_f,rank_f;
495:   PetscScalar            v[4];
496:   Mat                    mat;
497:   DMBoundaryType         bx,by;

500:   DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
501:   DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
502:   if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
503:   if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
504:   ratioi = mx/Mx;
505:   ratioj = my/My;
506:   if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x");
507:   if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y");
508:   if (ratioi != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 2");
509:   if (ratioj != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 2");

511:   DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
512:   DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
513:   DMGetLocalToGlobalMapping(daf,&ltog_f);
514:   ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);

516:   DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
517:   DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
518:   DMGetLocalToGlobalMapping(dac,&ltog_c);
519:   ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);

521:   /*
522:      Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
523:      The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
524:      processors). It's effective length is hence 4 times its normal length, this is
525:      why the col_scale is multiplied by the interpolation matrix column sizes.
526:      sol_shift allows each set of 1/4 processors do its own interpolation using ITS
527:      copy of the coarse vector. A bit of a hack but you do better.

529:      In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
530:   */
531:   MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);
532:   MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);
533:   MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);
534:   col_scale = size_f/size_c;
535:   col_shift = Mx*My*(rank_f/size_c);

537:   MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f,col_scale*m_c*n_c,dnz,onz);
538:   for (j=j_start; j<j_start+n_f; j++) {
539:     for (i=i_start; i<i_start+m_f; i++) {
540:       /* convert to local "natural" numbering and then to PETSc global numbering */
541:       row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];

543:       i_c = (i/ratioi);    /* coarse grid node to left of fine grid node */
544:       j_c = (j/ratioj);    /* coarse grid node below fine grid node */

546:       if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
547:     j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
548:       if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
549:     i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);

551:       /*
552:          Only include those interpolation points that are truly
553:          nonzero. Note this is very important for final grid lines
554:          in x and y directions; since they have no right/top neighbors
555:       */
556:       nc = 0;
557:       /* one left and below; or we are right on it */
558:       col        = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
559:       cols[nc++] = col_shift + idx_c[col];
560:       MatPreallocateSet(row,nc,cols,dnz,onz);
561:     }
562:   }
563:   MatCreate(PetscObjectComm((PetscObject)daf),&mat);
564:   MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);
565:   MatSetType(mat,MATAIJ);
566:   MatSeqAIJSetPreallocation(mat,0,dnz);
567:   MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
568:   MatPreallocateFinalize(dnz,onz);

570:   /* loop over local fine grid nodes setting interpolation for those*/
571:   for (j=j_start; j<j_start+n_f; j++) {
572:     for (i=i_start; i<i_start+m_f; i++) {
573:       /* convert to local "natural" numbering and then to PETSc global numbering */
574:       row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];

576:       i_c = (i/ratioi);    /* coarse grid node to left of fine grid node */
577:       j_c = (j/ratioj);    /* coarse grid node below fine grid node */
578:       nc  = 0;
579:       /* one left and below; or we are right on it */
580:       col      = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
581:       cols[nc] = col_shift + idx_c[col];
582:       v[nc++]  = 1.0;

584:       MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
585:     }
586:   }
587:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
588:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
589:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
590:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
591:   MatCreateMAIJ(mat,dof,A);
592:   MatDestroy(&mat);
593:   PetscLogFlops(13.0*m_f*n_f);
594:   return(0);
595: }

597: /*
598:        Contributed by Jianming Yang <jianming-yang@uiowa.edu>
599: */
600: PetscErrorCode DMCreateInterpolation_DA_3D_Q0(DM dac,DM daf,Mat *A)
601: {
602:   PetscErrorCode         ierr;
603:   PetscInt               i,j,l,i_start,j_start,l_start,m_f,n_f,p_f,Mx,My,Mz,dof;
604:   const PetscInt         *idx_c,*idx_f;
605:   ISLocalToGlobalMapping ltog_f,ltog_c;
606:   PetscInt               m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c,nc,*dnz,*onz;
607:   PetscInt               row,col,i_start_ghost,j_start_ghost,l_start_ghost,cols[8],mx,m_c,my,n_c,mz,p_c,ratioi,ratioj,ratiol;
608:   PetscInt               i_c,j_c,l_c,i_start_c,j_start_c,l_start_c,i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,col_shift,col_scale;
609:   PetscMPIInt            size_c,size_f,rank_f;
610:   PetscScalar            v[8];
611:   Mat                    mat;
612:   DMBoundaryType         bx,by,bz;

615:   DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
616:   if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
617:   if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
618:   if (!Mz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be positive",Mz);
619:   DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
620:   ratioi = mx/Mx;
621:   ratioj = my/My;
622:   ratiol = mz/Mz;
623:   if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x");
624:   if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y");
625:   if (ratiol*Mz != mz) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in z");
626:   if (ratioi != 2 && ratioi != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 1 or 2");
627:   if (ratioj != 2 && ratioj != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 1 or 2");
628:   if (ratiol != 2 && ratiol != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in z must be 1 or 2");

630:   DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
631:   DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
632:   DMGetLocalToGlobalMapping(daf,&ltog_f);
633:   ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);

635:   DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
636:   DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
637:   DMGetLocalToGlobalMapping(dac,&ltog_c);
638:   ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);

640:   /*
641:      Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA.
642:      The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
643:      processors). It's effective length is hence 4 times its normal length, this is
644:      why the col_scale is multiplied by the interpolation matrix column sizes.
645:      sol_shift allows each set of 1/4 processors do its own interpolation using ITS
646:      copy of the coarse vector. A bit of a hack but you do better.

648:      In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
649:   */
650:   MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);
651:   MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);
652:   MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);
653:   col_scale = size_f/size_c;
654:   col_shift = Mx*My*Mz*(rank_f/size_c);

656:   MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f*p_f,col_scale*m_c*n_c*p_c,dnz,onz);
657:   for (l=l_start; l<l_start+p_f; l++) {
658:     for (j=j_start; j<j_start+n_f; j++) {
659:       for (i=i_start; i<i_start+m_f; i++) {
660:         /* convert to local "natural" numbering and then to PETSc global numbering */
661:         row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];

663:         i_c = (i/ratioi);    /* coarse grid node to left of fine grid node */
664:         j_c = (j/ratioj);    /* coarse grid node below fine grid node */
665:         l_c = (l/ratiol);

667:         if (l_c < l_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
668:     l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
669:         if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
670:     j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
671:         if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
672:     i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);

674:         /*
675:            Only include those interpolation points that are truly
676:            nonzero. Note this is very important for final grid lines
677:            in x and y directions; since they have no right/top neighbors
678:         */
679:         nc = 0;
680:         /* one left and below; or we are right on it */
681:         col        = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
682:         cols[nc++] = col_shift + idx_c[col];
683:         MatPreallocateSet(row,nc,cols,dnz,onz);
684:       }
685:     }
686:   }
687:   MatCreate(PetscObjectComm((PetscObject)daf),&mat);
688:   MatSetSizes(mat,m_f*n_f*p_f,col_scale*m_c*n_c*p_c,mx*my*mz,col_scale*Mx*My*Mz);
689:   MatSetType(mat,MATAIJ);
690:   MatSeqAIJSetPreallocation(mat,0,dnz);
691:   MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
692:   MatPreallocateFinalize(dnz,onz);

694:   /* loop over local fine grid nodes setting interpolation for those*/
695:   for (l=l_start; l<l_start+p_f; l++) {
696:     for (j=j_start; j<j_start+n_f; j++) {
697:       for (i=i_start; i<i_start+m_f; i++) {
698:         /* convert to local "natural" numbering and then to PETSc global numbering */
699:         row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];

701:         i_c = (i/ratioi);    /* coarse grid node to left of fine grid node */
702:         j_c = (j/ratioj);    /* coarse grid node below fine grid node */
703:         l_c = (l/ratiol);
704:         nc  = 0;
705:         /* one left and below; or we are right on it */
706:         col      = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
707:         cols[nc] = col_shift + idx_c[col];
708:         v[nc++]  = 1.0;

710:         MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
711:       }
712:     }
713:   }
714:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
715:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
716:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
717:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
718:   MatCreateMAIJ(mat,dof,A);
719:   MatDestroy(&mat);
720:   PetscLogFlops(13.0*m_f*n_f*p_f);
721:   return(0);
722: }

724: PetscErrorCode DMCreateInterpolation_DA_3D_Q1(DM dac,DM daf,Mat *A)
725: {
726:   PetscErrorCode         ierr;
727:   PetscInt               i,j,i_start,j_start,m_f,n_f,Mx,My,dof,l;
728:   const PetscInt         *idx_c,*idx_f;
729:   ISLocalToGlobalMapping ltog_f,ltog_c;
730:   PetscInt               m_ghost,n_ghost,m_ghost_c,n_ghost_c,Mz,mz;
731:   PetscInt               row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok;
732:   PetscInt               i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
733:   PetscInt               l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c;
734:   PetscInt               l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz;
735:   PetscScalar            v[8],x,y,z;
736:   Mat                    mat;
737:   DMBoundaryType         bx,by,bz;

740:   DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
741:   DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
742:   if (mx == Mx) {
743:     ratioi = 1;
744:   } else if (bx == DM_BOUNDARY_PERIODIC) {
745:     if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx);
746:     ratioi = mx/Mx;
747:     if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
748:   } else {
749:     if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx);
750:     ratioi = (mx-1)/(Mx-1);
751:     if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
752:   }
753:   if (my == My) {
754:     ratioj = 1;
755:   } else if (by == DM_BOUNDARY_PERIODIC) {
756:     if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My);
757:     ratioj = my/My;
758:     if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My  must be integer: my %D My %D",my,My);
759:   } else {
760:     if (My < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be greater than 1",My);
761:     ratioj = (my-1)/(My-1);
762:     if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
763:   }
764:   if (mz == Mz) {
765:     ratiok = 1;
766:   } else if (bz == DM_BOUNDARY_PERIODIC) {
767:     if (!Mz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be positive",Mz);
768:     ratiok = mz/Mz;
769:     if (ratiok*Mz != mz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mz/Mz  must be integer: mz %D Mz %D",mz,Mz);
770:   } else {
771:     if (Mz < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be greater than 1",Mz);
772:     ratiok = (mz-1)/(Mz-1);
773:     if (ratiok*(Mz-1) != mz-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %D Mz %D",mz,Mz);
774:   }

776:   DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
777:   DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
778:   DMGetLocalToGlobalMapping(daf,&ltog_f);
779:   ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);

781:   DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
782:   DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
783:   DMGetLocalToGlobalMapping(dac,&ltog_c);
784:   ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);

786:   /* create interpolation matrix, determining exact preallocation */
787:   MatPreallocateInitialize(PetscObjectComm((PetscObject)dac),m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);
788:   /* loop over local fine grid nodes counting interpolating points */
789:   for (l=l_start; l<l_start+p_f; l++) {
790:     for (j=j_start; j<j_start+n_f; j++) {
791:       for (i=i_start; i<i_start+m_f; i++) {
792:         /* convert to local "natural" numbering and then to PETSc global numbering */
793:         row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
794:         i_c = (i/ratioi);
795:         j_c = (j/ratioj);
796:         l_c = (l/ratiok);
797:         if (l_c < l_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
798:                                             l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
799:         if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
800:                                             j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
801:         if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
802:                                             i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);

804:         /*
805:          Only include those interpolation points that are truly
806:          nonzero. Note this is very important for final grid lines
807:          in x and y directions; since they have no right/top neighbors
808:          */
809:         nc         = 0;
810:         col        = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
811:         cols[nc++] = idx_c[col];
812:         if (i_c*ratioi != i) {
813:           cols[nc++] = idx_c[col+1];
814:         }
815:         if (j_c*ratioj != j) {
816:           cols[nc++] = idx_c[col+m_ghost_c];
817:         }
818:         if (l_c*ratiok != l) {
819:           cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c];
820:         }
821:         if (j_c*ratioj != j && i_c*ratioi != i) {
822:           cols[nc++] = idx_c[col+(m_ghost_c+1)];
823:         }
824:         if (j_c*ratioj != j && l_c*ratiok != l) {
825:           cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)];
826:         }
827:         if (i_c*ratioi != i && l_c*ratiok != l) {
828:           cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)];
829:         }
830:         if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
831:           cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)];
832:         }
833:         MatPreallocateSet(row,nc,cols,dnz,onz);
834:       }
835:     }
836:   }
837:   MatCreate(PetscObjectComm((PetscObject)dac),&mat);
838:   MatSetSizes(mat,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz);
839:   MatSetType(mat,MATAIJ);
840:   MatSeqAIJSetPreallocation(mat,0,dnz);
841:   MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
842:   MatPreallocateFinalize(dnz,onz);

844:   /* loop over local fine grid nodes setting interpolation for those*/
845:   if (!NEWVERSION) {

847:     for (l=l_start; l<l_start+p_f; l++) {
848:       for (j=j_start; j<j_start+n_f; j++) {
849:         for (i=i_start; i<i_start+m_f; i++) {
850:           /* convert to local "natural" numbering and then to PETSc global numbering */
851:           row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];

853:           i_c = (i/ratioi);
854:           j_c = (j/ratioj);
855:           l_c = (l/ratiok);

857:           /*
858:            Only include those interpolation points that are truly
859:            nonzero. Note this is very important for final grid lines
860:            in x and y directions; since they have no right/top neighbors
861:            */
862:           x = ((PetscReal)(i - i_c*ratioi))/((PetscReal)ratioi);
863:           y = ((PetscReal)(j - j_c*ratioj))/((PetscReal)ratioj);
864:           z = ((PetscReal)(l - l_c*ratiok))/((PetscReal)ratiok);

866:           nc = 0;
867:           /* one left and below; or we are right on it */
868:           col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));

870:           cols[nc] = idx_c[col];
871:           v[nc++]  = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));

873:           if (i_c*ratioi != i) {
874:             cols[nc] = idx_c[col+1];
875:             v[nc++]  = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
876:           }

878:           if (j_c*ratioj != j) {
879:             cols[nc] = idx_c[col+m_ghost_c];
880:             v[nc++]  = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
881:           }

883:           if (l_c*ratiok != l) {
884:             cols[nc] = idx_c[col+m_ghost_c*n_ghost_c];
885:             v[nc++]  = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
886:           }

888:           if (j_c*ratioj != j && i_c*ratioi != i) {
889:             cols[nc] = idx_c[col+(m_ghost_c+1)];
890:             v[nc++]  = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
891:           }

893:           if (j_c*ratioj != j && l_c*ratiok != l) {
894:             cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)];
895:             v[nc++]  = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
896:           }

898:           if (i_c*ratioi != i && l_c*ratiok != l) {
899:             cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)];
900:             v[nc++]  = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
901:           }

903:           if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
904:             cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)];
905:             v[nc++]  = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
906:           }
907:           MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
908:         }
909:       }
910:     }

912:   } else {
913:     PetscScalar *xi,*eta,*zeta;
914:     PetscInt    li,nxi,lj,neta,lk,nzeta,n;
915:     PetscScalar Ni[8];

917:     /* compute local coordinate arrays */
918:     nxi   = ratioi + 1;
919:     neta  = ratioj + 1;
920:     nzeta = ratiok + 1;
921:     PetscMalloc1(nxi,&xi);
922:     PetscMalloc1(neta,&eta);
923:     PetscMalloc1(nzeta,&zeta);
924:     for (li=0; li<nxi; li++) xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
925:     for (lj=0; lj<neta; lj++) eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
926:     for (lk=0; lk<nzeta; lk++) zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));

928:     for (l=l_start; l<l_start+p_f; l++) {
929:       for (j=j_start; j<j_start+n_f; j++) {
930:         for (i=i_start; i<i_start+m_f; i++) {
931:           /* convert to local "natural" numbering and then to PETSc global numbering */
932:           row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];

934:           i_c = (i/ratioi);
935:           j_c = (j/ratioj);
936:           l_c = (l/ratiok);

938:           /* remainders */
939:           li = i - ratioi * (i/ratioi);
940:           if (i==mx-1) li = nxi-1;
941:           lj = j - ratioj * (j/ratioj);
942:           if (j==my-1) lj = neta-1;
943:           lk = l - ratiok * (l/ratiok);
944:           if (l==mz-1) lk = nzeta-1;

946:           /* corners */
947:           col     = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));
948:           cols[0] = idx_c[col];
949:           Ni[0]   = 1.0;
950:           if ((li==0) || (li==nxi-1)) {
951:             if ((lj==0) || (lj==neta-1)) {
952:               if ((lk==0) || (lk==nzeta-1)) {
953:                 MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);
954:                 continue;
955:               }
956:             }
957:           }

959:           /* edges + interior */
960:           /* remainders */
961:           if (i==mx-1) i_c--;
962:           if (j==my-1) j_c--;
963:           if (l==mz-1) l_c--;

965:           col     = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
966:           cols[0] = idx_c[col]; /* one left and below; or we are right on it */
967:           cols[1] = idx_c[col+1]; /* one right and below */
968:           cols[2] = idx_c[col+m_ghost_c];  /* one left and above */
969:           cols[3] = idx_c[col+(m_ghost_c+1)]; /* one right and above */

971:           cols[4] = idx_c[col+m_ghost_c*n_ghost_c]; /* one left and below and front; or we are right on it */
972:           cols[5] = idx_c[col+(m_ghost_c*n_ghost_c+1)]; /* one right and below, and front */
973:           cols[6] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)]; /* one left and above and front*/
974:           cols[7] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)]; /* one right and above and front */

976:           Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
977:           Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
978:           Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
979:           Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);

981:           Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
982:           Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
983:           Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
984:           Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);

986:           for (n=0; n<8; n++) {
987:             if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1;
988:           }
989:           MatSetValues(mat,1,&row,8,cols,Ni,INSERT_VALUES);

991:         }
992:       }
993:     }
994:     PetscFree(xi);
995:     PetscFree(eta);
996:     PetscFree(zeta);
997:   }
998:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
999:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);
1000:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
1001:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

1003:   MatCreateMAIJ(mat,dof,A);
1004:   MatDestroy(&mat);
1005:   return(0);
1006: }

1008: PetscErrorCode  DMCreateInterpolation_DA(DM dac,DM daf,Mat *A,Vec *scale)
1009: {
1010:   PetscErrorCode   ierr;
1011:   PetscInt         dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1012:   DMBoundaryType   bxc,byc,bzc,bxf,byf,bzf;
1013:   DMDAStencilType  stc,stf;
1014:   DM_DA            *ddc = (DM_DA*)dac->data;


1022:   DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1023:   DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1024:   if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1025:   if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1026:   if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1027:   if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1028:   if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1029:   if (Mc < 2 && Mf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
1030:   if (dimc > 1 && Nc < 2 && Nf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction");
1031:   if (dimc > 2 && Pc < 2 && Pf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction");

1033:   if (ddc->interptype == DMDA_Q1) {
1034:     if (dimc == 1) {
1035:       DMCreateInterpolation_DA_1D_Q1(dac,daf,A);
1036:     } else if (dimc == 2) {
1037:       DMCreateInterpolation_DA_2D_Q1(dac,daf,A);
1038:     } else if (dimc == 3) {
1039:       DMCreateInterpolation_DA_3D_Q1(dac,daf,A);
1040:     } else SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype);
1041:   } else if (ddc->interptype == DMDA_Q0) {
1042:     if (dimc == 1) {
1043:       DMCreateInterpolation_DA_1D_Q0(dac,daf,A);
1044:     } else if (dimc == 2) {
1045:       DMCreateInterpolation_DA_2D_Q0(dac,daf,A);
1046:     } else if (dimc == 3) {
1047:       DMCreateInterpolation_DA_3D_Q0(dac,daf,A);
1048:     } else SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype);
1049:   }
1050:   if (scale) {
1051:     DMCreateInterpolationScale((DM)dac,(DM)daf,*A,scale);
1052:   }
1053:   return(0);
1054: }

1056: PetscErrorCode DMCreateInjection_DA_1D(DM dac,DM daf,VecScatter *inject)
1057: {
1058:   PetscErrorCode         ierr;
1059:   PetscInt               i,i_start,m_f,Mx,dof;
1060:   const PetscInt         *idx_f;
1061:   ISLocalToGlobalMapping ltog_f;
1062:   PetscInt               m_ghost,m_ghost_c;
1063:   PetscInt               row,i_start_ghost,mx,m_c,nc,ratioi;
1064:   PetscInt               i_start_c,i_start_ghost_c;
1065:   PetscInt               *cols;
1066:   DMBoundaryType         bx;
1067:   Vec                    vecf,vecc;
1068:   IS                     isf;

1071:   DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);
1072:   DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);
1073:   if (bx == DM_BOUNDARY_PERIODIC) {
1074:     ratioi = mx/Mx;
1075:     if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
1076:   } else {
1077:     ratioi = (mx-1)/(Mx-1);
1078:     if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
1079:   }

1081:   DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);
1082:   DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
1083:   DMGetLocalToGlobalMapping(daf,&ltog_f);
1084:   ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);

1086:   DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
1087:   DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);


1090:   /* loop over local fine grid nodes setting interpolation for those*/
1091:   nc   = 0;
1092:   PetscMalloc1(m_f,&cols);


1095:   for (i=i_start_c; i<i_start_c+m_c; i++) {
1096:     PetscInt i_f = i*ratioi;

1098:     if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\ni_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);

1100:     row        = idx_f[(i_f-i_start_ghost)];
1101:     cols[nc++] = row;
1102:   }

1104:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
1105:   ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);
1106:   DMGetGlobalVector(dac,&vecc);
1107:   DMGetGlobalVector(daf,&vecf);
1108:   VecScatterCreate(vecf,isf,vecc,NULL,inject);
1109:   DMRestoreGlobalVector(dac,&vecc);
1110:   DMRestoreGlobalVector(daf,&vecf);
1111:   ISDestroy(&isf);
1112:   return(0);
1113: }

1115: PetscErrorCode DMCreateInjection_DA_2D(DM dac,DM daf,VecScatter *inject)
1116: {
1117:   PetscErrorCode         ierr;
1118:   PetscInt               i,j,i_start,j_start,m_f,n_f,Mx,My,dof;
1119:   const PetscInt         *idx_c,*idx_f;
1120:   ISLocalToGlobalMapping ltog_f,ltog_c;
1121:   PetscInt               m_ghost,n_ghost,m_ghost_c,n_ghost_c;
1122:   PetscInt               row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj;
1123:   PetscInt               i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
1124:   PetscInt               *cols;
1125:   DMBoundaryType         bx,by;
1126:   Vec                    vecf,vecc;
1127:   IS                     isf;

1130:   DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);
1131:   DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);
1132:   if (bx == DM_BOUNDARY_PERIODIC) {
1133:     ratioi = mx/Mx;
1134:     if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
1135:   } else {
1136:     ratioi = (mx-1)/(Mx-1);
1137:     if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
1138:   }
1139:   if (by == DM_BOUNDARY_PERIODIC) {
1140:     ratioj = my/My;
1141:     if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My  must be integer: my %D My %D",my,My);
1142:   } else {
1143:     ratioj = (my-1)/(My-1);
1144:     if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
1145:   }

1147:   DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
1148:   DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
1149:   DMGetLocalToGlobalMapping(daf,&ltog_f);
1150:   ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);

1152:   DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
1153:   DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
1154:   DMGetLocalToGlobalMapping(dac,&ltog_c);
1155:   ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);

1157:   /* loop over local fine grid nodes setting interpolation for those*/
1158:   nc   = 0;
1159:   PetscMalloc1(n_f*m_f,&cols);
1160:   for (j=j_start_c; j<j_start_c+n_c; j++) {
1161:     for (i=i_start_c; i<i_start_c+m_c; i++) {
1162:       PetscInt i_f = i*ratioi,j_f = j*ratioj;
1163:       if (j_f < j_start_ghost || j_f >= j_start_ghost+n_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
1164:     j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost);
1165:       if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\
1166:     i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1167:       row        = idx_f[(m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))];
1168:       cols[nc++] = row;
1169:     }
1170:   }
1171:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
1172:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);

1174:   ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);
1175:   DMGetGlobalVector(dac,&vecc);
1176:   DMGetGlobalVector(daf,&vecf);
1177:   VecScatterCreate(vecf,isf,vecc,NULL,inject);
1178:   DMRestoreGlobalVector(dac,&vecc);
1179:   DMRestoreGlobalVector(daf,&vecf);
1180:   ISDestroy(&isf);
1181:   return(0);
1182: }

1184: PetscErrorCode DMCreateInjection_DA_3D(DM dac,DM daf,VecScatter *inject)
1185: {
1186:   PetscErrorCode         ierr;
1187:   PetscInt               i,j,k,i_start,j_start,k_start,m_f,n_f,p_f,Mx,My,Mz;
1188:   PetscInt               m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c;
1189:   PetscInt               i_start_ghost,j_start_ghost,k_start_ghost;
1190:   PetscInt               mx,my,mz,ratioi,ratioj,ratiok;
1191:   PetscInt               i_start_c,j_start_c,k_start_c;
1192:   PetscInt               m_c,n_c,p_c;
1193:   PetscInt               i_start_ghost_c,j_start_ghost_c,k_start_ghost_c;
1194:   PetscInt               row,nc,dof;
1195:   const PetscInt         *idx_c,*idx_f;
1196:   ISLocalToGlobalMapping ltog_f,ltog_c;
1197:   PetscInt               *cols;
1198:   DMBoundaryType         bx,by,bz;
1199:   Vec                    vecf,vecc;
1200:   IS                     isf;

1203:   DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);
1204:   DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);

1206:   if (bx == DM_BOUNDARY_PERIODIC) {
1207:     ratioi = mx/Mx;
1208:     if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
1209:   } else {
1210:     ratioi = (mx-1)/(Mx-1);
1211:     if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
1212:   }
1213:   if (by == DM_BOUNDARY_PERIODIC) {
1214:     ratioj = my/My;
1215:     if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My  must be integer: my %D My %D",my,My);
1216:   } else {
1217:     ratioj = (my-1)/(My-1);
1218:     if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
1219:   }
1220:   if (bz == DM_BOUNDARY_PERIODIC) {
1221:     ratiok = mz/Mz;
1222:     if (ratiok*Mz != mz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mz/Mz  must be integer: mz %D My %D",mz,Mz);
1223:   } else {
1224:     ratiok = (mz-1)/(Mz-1);
1225:     if (ratiok*(Mz-1) != mz-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %D Mz %D",mz,Mz);
1226:   }

1228:   DMDAGetCorners(daf,&i_start,&j_start,&k_start,&m_f,&n_f,&p_f);
1229:   DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&k_start_ghost,&m_ghost,&n_ghost,&p_ghost);
1230:   DMGetLocalToGlobalMapping(daf,&ltog_f);
1231:   ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);

1233:   DMDAGetCorners(dac,&i_start_c,&j_start_c,&k_start_c,&m_c,&n_c,&p_c);
1234:   DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&k_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
1235:   DMGetLocalToGlobalMapping(dac,&ltog_c);
1236:   ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);


1239:   /* loop over local fine grid nodes setting interpolation for those*/
1240:   nc   = 0;
1241:   PetscMalloc1(n_f*m_f*p_f,&cols);
1242:   for (k=k_start_c; k<k_start_c+p_c; k++) {
1243:     for (j=j_start_c; j<j_start_c+n_c; j++) {
1244:       for (i=i_start_c; i<i_start_c+m_c; i++) {
1245:         PetscInt i_f = i*ratioi,j_f = j*ratioj,k_f = k*ratiok;
1246:         if (k_f < k_start_ghost || k_f >= k_start_ghost+p_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA  "
1247:                                                                           "k_c %D k_f %D fine ghost range [%D,%D]",k,k_f,k_start_ghost,k_start_ghost+p_ghost);
1248:         if (j_f < j_start_ghost || j_f >= j_start_ghost+n_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA  "
1249:                                                                           "j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost);
1250:         if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA  "
1251:                                                                           "i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost);
1252:         row        = idx_f[(m_ghost*n_ghost*(k_f-k_start_ghost) + m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))];
1253:         cols[nc++] = row;
1254:       }
1255:     }
1256:   }
1257:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);
1258:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);

1260:   ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);
1261:   DMGetGlobalVector(dac,&vecc);
1262:   DMGetGlobalVector(daf,&vecf);
1263:   VecScatterCreate(vecf,isf,vecc,NULL,inject);
1264:   DMRestoreGlobalVector(dac,&vecc);
1265:   DMRestoreGlobalVector(daf,&vecf);
1266:   ISDestroy(&isf);
1267:   return(0);
1268: }

1270: PetscErrorCode  DMCreateInjection_DA(DM dac,DM daf,Mat *mat)
1271: {
1272:   PetscErrorCode  ierr;
1273:   PetscInt        dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1274:   DMBoundaryType  bxc,byc,bzc,bxf,byf,bzf;
1275:   DMDAStencilType stc,stf;
1276:   VecScatter      inject = NULL;


1283:   DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1284:   DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1285:   if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1286:   if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1287:   if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1288:   if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1289:   if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");
1290:   if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
1291:   if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction");
1292:   if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction");

1294:   if (dimc == 1) {
1295:     DMCreateInjection_DA_1D(dac,daf,&inject);
1296:   } else if (dimc == 2) {
1297:     DMCreateInjection_DA_2D(dac,daf,&inject);
1298:   } else if (dimc == 3) {
1299:     DMCreateInjection_DA_3D(dac,daf,&inject);
1300:   }
1301:   MatCreateScatter(PetscObjectComm((PetscObject)inject), inject, mat);
1302:   VecScatterDestroy(&inject);
1303:   return(0);
1304: }

1306: PetscErrorCode  DMCreateAggregates_DA(DM dac,DM daf,Mat *rest)
1307: {
1308:   PetscErrorCode         ierr;
1309:   PetscInt               dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc;
1310:   PetscInt               dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
1311:   DMBoundaryType         bxc,byc,bzc,bxf,byf,bzf;
1312:   DMDAStencilType        stc,stf;
1313:   PetscInt               i,j,l;
1314:   PetscInt               i_start,j_start,l_start, m_f,n_f,p_f;
1315:   PetscInt               i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost;
1316:   const PetscInt         *idx_f;
1317:   PetscInt               i_c,j_c,l_c;
1318:   PetscInt               i_start_c,j_start_c,l_start_c, m_c,n_c,p_c;
1319:   PetscInt               i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c;
1320:   const PetscInt         *idx_c;
1321:   PetscInt               d;
1322:   PetscInt               a;
1323:   PetscInt               max_agg_size;
1324:   PetscInt               *fine_nodes;
1325:   PetscScalar            *one_vec;
1326:   PetscInt               fn_idx;
1327:   ISLocalToGlobalMapping ltogmf,ltogmc;


1334:   DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);
1335:   DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);
1336:   if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);
1337:   if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);
1338:   if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);
1339:   if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");
1340:   if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");

1342:   if (Mf < Mc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Mc %D, Mf %D", Mc, Mf);
1343:   if (Nf < Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Nc %D, Nf %D", Nc, Nf);
1344:   if (Pf < Pc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Pc %D, Pf %D", Pc, Pf);

1346:   if (Pc < 0) Pc = 1;
1347:   if (Pf < 0) Pf = 1;
1348:   if (Nc < 0) Nc = 1;
1349:   if (Nf < 0) Nf = 1;

1351:   DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
1352:   DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);

1354:   DMGetLocalToGlobalMapping(daf,&ltogmf);
1355:   ISLocalToGlobalMappingGetIndices(ltogmf,&idx_f);

1357:   DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
1358:   DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);

1360:   DMGetLocalToGlobalMapping(dac,&ltogmc);
1361:   ISLocalToGlobalMappingGetIndices(ltogmc,&idx_c);

1363:   /*
1364:      Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios
1365:      for dimension 1 and 2 respectively.
1366:      Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1)
1367:      and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate.
1368:      Each specific dof on the fine grid is mapped to one dof on the coarse grid.
1369:   */

1371:   max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1);

1373:   /* create the matrix that will contain the restriction operator */
1374:   MatCreateAIJ(PetscObjectComm((PetscObject)daf), m_c*n_c*p_c*dofc, m_f*n_f*p_f*doff, Mc*Nc*Pc*dofc, Mf*Nf*Pf*doff,
1375:                       max_agg_size, NULL, max_agg_size, NULL, rest);

1377:   /* store nodes in the fine grid here */
1378:   PetscMalloc2(max_agg_size, &one_vec,max_agg_size, &fine_nodes);
1379:   for (i=0; i<max_agg_size; i++) one_vec[i] = 1.0;

1381:   /* loop over all coarse nodes */
1382:   for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) {
1383:     for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) {
1384:       for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) {
1385:         for (d=0; d<dofc; d++) {
1386:           /* convert to local "natural" numbering and then to PETSc global numbering */
1387:           a = idx_c[dofc*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c))] + d;

1389:           fn_idx = 0;
1390:           /* Corresponding fine points are all points (i_f, j_f, l_f) such that
1391:              i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc
1392:              (same for other dimensions)
1393:           */
1394:           for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) {
1395:             for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) {
1396:               for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) {
1397:                 fine_nodes[fn_idx] = idx_f[doff*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))] + d;
1398:                 fn_idx++;
1399:               }
1400:             }
1401:           }
1402:           /* add all these points to one aggregate */
1403:           MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);
1404:         }
1405:       }
1406:     }
1407:   }
1408:   ISLocalToGlobalMappingRestoreIndices(ltogmf,&idx_f);
1409:   ISLocalToGlobalMappingRestoreIndices(ltogmc,&idx_c);
1410:   PetscFree2(one_vec,fine_nodes);
1411:   MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);
1412:   MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);
1413:   return(0);
1414: }