Actual source code: dainterp.c


  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:    Since the interpolation uses MATMAIJ for dof > 0 we convert request for non-MATAIJ baseded matrices to MATAIJ.
 18:    This is a bit of a hack, the reason for it is partially because -dm_mat_type defines the
 19:    matrix type for both the operator matrices and the interpolation matrices so that users
 20:    can select matrix types of base MATAIJ for accelerators
 21: */
 22: static PetscErrorCode ConvertToAIJ(MatType intype, MatType *outtype)
 23: {
 24:   PetscInt    i;
 25:   char const *types[3] = {MATAIJ, MATSEQAIJ, MATMPIAIJ};
 26:   PetscBool   flg;

 28:   *outtype = MATAIJ;
 29:   for (i = 0; i < 3; i++) {
 30:     PetscStrbeginswith(intype, types[i], &flg);
 31:     if (flg) {
 32:       *outtype = intype;
 33:       return 0;
 34:     }
 35:   }
 36:   return 0;
 37: }

 39: PetscErrorCode DMCreateInterpolation_DA_1D_Q1(DM dac, DM daf, Mat *A)
 40: {
 41:   PetscInt               i, i_start, m_f, Mx;
 42:   const PetscInt        *idx_f, *idx_c;
 43:   PetscInt               m_ghost, m_ghost_c;
 44:   PetscInt               row, col, i_start_ghost, mx, m_c, nc, ratio;
 45:   PetscInt               i_c, i_start_c, i_start_ghost_c, cols[2], dof;
 46:   PetscScalar            v[2], x;
 47:   Mat                    mat;
 48:   DMBoundaryType         bx;
 49:   ISLocalToGlobalMapping ltog_f, ltog_c;
 50:   MatType                mattype;

 52:   DMDAGetInfo(dac, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bx, NULL, NULL, NULL);
 53:   DMDAGetInfo(daf, NULL, &mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
 54:   if (bx == DM_BOUNDARY_PERIODIC) {
 55:     ratio = mx / Mx;
 57:   } else {
 58:     ratio = (mx - 1) / (Mx - 1);
 60:   }

 62:   DMDAGetCorners(daf, &i_start, NULL, NULL, &m_f, NULL, NULL);
 63:   DMDAGetGhostCorners(daf, &i_start_ghost, NULL, NULL, &m_ghost, NULL, NULL);
 64:   DMGetLocalToGlobalMapping(daf, &ltog_f);
 65:   ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f);

 67:   DMDAGetCorners(dac, &i_start_c, NULL, NULL, &m_c, NULL, NULL);
 68:   DMDAGetGhostCorners(dac, &i_start_ghost_c, NULL, NULL, &m_ghost_c, NULL, NULL);
 69:   DMGetLocalToGlobalMapping(dac, &ltog_c);
 70:   ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c);

 72:   /* create interpolation matrix */
 73:   MatCreate(PetscObjectComm((PetscObject)dac), &mat);
 74: #if defined(PETSC_HAVE_CUDA)
 75:   /*
 76:      Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
 77:      we don't want the original unconverted matrix copied to the GPU
 78:    */
 79:   if (dof > 1) MatBindToCPU(mat, PETSC_TRUE);
 80: #endif
 81:   MatSetSizes(mat, m_f, m_c, mx, Mx);
 82:   ConvertToAIJ(dac->mattype, &mattype);
 83:   MatSetType(mat, mattype);
 84:   MatSeqAIJSetPreallocation(mat, 2, NULL);
 85:   MatMPIAIJSetPreallocation(mat, 2, NULL, 1, NULL);

 87:   /* loop over local fine grid nodes setting interpolation for those*/
 88:   if (!NEWVERSION) {
 89:     for (i = i_start; i < i_start + m_f; i++) {
 90:       /* convert to local "natural" numbering and then to PETSc global numbering */
 91:       row = idx_f[i - i_start_ghost];

 93:       i_c = (i / ratio); /* coarse grid node to left of fine grid node */
 95:                                           i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT,
 96:                  i_start, i_c, i_start_ghost_c);

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

117:   } else {
118:     PetscScalar *xi;
119:     PetscInt     li, nxi, n;
120:     PetscScalar  Ni[2];

122:     /* compute local coordinate arrays */
123:     nxi = ratio + 1;
124:     PetscMalloc1(nxi, &xi);
125:     for (li = 0; li < nxi; li++) xi[li] = -1.0 + (PetscScalar)li * (2.0 / (PetscScalar)(nxi - 1));

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

131:       i_c = (i / ratio); /* coarse grid node to left of fine grid node */
133:                                           i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT,
134:                  i_start, i_c, i_start_ghost_c);

136:       /* remainders */
137:       li = i - ratio * (i / ratio);
138:       if (i == mx - 1) li = nxi - 1;

140:       /* corners */
141:       col     = (i_c - i_start_ghost_c);
142:       cols[0] = idx_c[col];
143:       Ni[0]   = 1.0;
144:       if ((li == 0) || (li == nxi - 1)) {
145:         MatSetValue(mat, row, cols[0], Ni[0], INSERT_VALUES);
146:         continue;
147:       }

149:       /* edges + interior */
150:       /* remainders */
151:       if (i == mx - 1) i_c--;

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

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

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

188:   DMDAGetInfo(dac, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bx, NULL, NULL, NULL);
189:   DMDAGetInfo(daf, NULL, &mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
190:   if (bx == DM_BOUNDARY_PERIODIC) {
192:     ratio = mx / Mx;
194:   } else {
196:     ratio = (mx - 1) / (Mx - 1);
198:   }

200:   DMDAGetCorners(daf, &i_start, NULL, NULL, &m_f, NULL, NULL);
201:   DMDAGetGhostCorners(daf, &i_start_ghost, NULL, NULL, &m_ghost, NULL, NULL);
202:   DMGetLocalToGlobalMapping(daf, &ltog_f);
203:   ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f);

205:   DMDAGetCorners(dac, &i_start_c, NULL, NULL, &m_c, NULL, NULL);
206:   DMDAGetGhostCorners(dac, &i_start_ghost_c, NULL, NULL, &m_ghost_c, NULL, NULL);
207:   DMGetLocalToGlobalMapping(dac, &ltog_c);
208:   ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c);

210:   /* create interpolation matrix */
211:   MatCreate(PetscObjectComm((PetscObject)dac), &mat);
212: #if defined(PETSC_HAVE_CUDA)
213:   /*
214:      Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
215:      we don't want the original unconverted matrix copied to the GPU
216:    */
217:   if (dof > 1) MatBindToCPU(mat, PETSC_TRUE);
218: #endif
219:   MatSetSizes(mat, m_f, m_c, mx, Mx);
220:   ConvertToAIJ(dac->mattype, &mattype);
221:   MatSetType(mat, mattype);
222:   MatSeqAIJSetPreallocation(mat, 2, NULL);
223:   MatMPIAIJSetPreallocation(mat, 2, NULL, 0, NULL);

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

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

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

260: PetscErrorCode DMCreateInterpolation_DA_2D_Q1(DM dac, DM daf, Mat *A)
261: {
262:   PetscInt               i, j, i_start, j_start, m_f, n_f, Mx, My, dof;
263:   const PetscInt        *idx_c, *idx_f;
264:   ISLocalToGlobalMapping ltog_f, ltog_c;
265:   PetscInt               m_ghost, n_ghost, m_ghost_c, n_ghost_c, *dnz, *onz;
266:   PetscInt               row, col, i_start_ghost, j_start_ghost, cols[4], mx, m_c, my, nc, ratioi, ratioj;
267:   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;
268:   PetscMPIInt            size_c, size_f, rank_f;
269:   PetscScalar            v[4], x, y;
270:   Mat                    mat;
271:   DMBoundaryType         bx, by;
272:   MatType                mattype;

274:   DMDAGetInfo(dac, NULL, &Mx, &My, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL);
275:   DMDAGetInfo(daf, NULL, &mx, &my, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
276:   if (bx == DM_BOUNDARY_PERIODIC) {
278:     ratioi = mx / Mx;
280:   } else {
282:     ratioi = (mx - 1) / (Mx - 1);
284:   }
285:   if (by == DM_BOUNDARY_PERIODIC) {
287:     ratioj = my / My;
289:   } else {
291:     ratioj = (my - 1) / (My - 1);
293:   }

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

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

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

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

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

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

331:                                           j_start %" PetscInt_FMT " j_c %" PetscInt_FMT " j_start_ghost_c %" PetscInt_FMT,
332:                  j_start, j_c, j_start_ghost_c);
334:                                           i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT,
335:                  i_start, i_c, i_start_ghost_c);

337:       /*
338:        Only include those interpolation points that are truly
339:        nonzero. Note this is very important for final grid lines
340:        in x and y directions; since they have no right/top neighbors
341:        */
342:       nc = 0;
343:       /* one left and below; or we are right on it */
344:       col        = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c));
345:       cols[nc++] = col_shift + idx_c[col];
346:       /* one right and below */
347:       if (i_c * ratioi != i) cols[nc++] = col_shift + idx_c[col + 1];
348:       /* one left and above */
349:       if (j_c * ratioj != j) cols[nc++] = col_shift + idx_c[col + m_ghost_c];
350:       /* one right and above */
351:       if (i_c * ratioi != i && j_c * ratioj != j) cols[nc++] = col_shift + idx_c[col + (m_ghost_c + 1)];
352:       MatPreallocateSet(row, nc, cols, dnz, onz);
353:     }
354:   }
355:   MatCreate(PetscObjectComm((PetscObject)daf), &mat);
356: #if defined(PETSC_HAVE_CUDA)
357:   /*
358:      Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
359:      we don't want the original unconverted matrix copied to the GPU
360:   */
361:   if (dof > 1) MatBindToCPU(mat, PETSC_TRUE);
362: #endif
363:   MatSetSizes(mat, m_f * n_f, col_scale * m_c * n_c, mx * my, col_scale * Mx * My);
364:   ConvertToAIJ(dac->mattype, &mattype);
365:   MatSetType(mat, mattype);
366:   MatSeqAIJSetPreallocation(mat, 0, dnz);
367:   MatMPIAIJSetPreallocation(mat, 0, dnz, 0, onz);
368:   MatPreallocateEnd(dnz, onz);

370:   /* loop over local fine grid nodes setting interpolation for those*/
371:   if (!NEWVERSION) {
372:     for (j = j_start; j < j_start + n_f; j++) {
373:       for (i = i_start; i < i_start + m_f; i++) {
374:         /* convert to local "natural" numbering and then to PETSc global numbering */
375:         row = idx_f[(m_ghost * (j - j_start_ghost) + (i - i_start_ghost))];

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

380:         /*
381:          Only include those interpolation points that are truly
382:          nonzero. Note this is very important for final grid lines
383:          in x and y directions; since they have no right/top neighbors
384:          */
385:         x = ((PetscReal)(i - i_c * ratioi)) / ((PetscReal)ratioi);
386:         y = ((PetscReal)(j - j_c * ratioj)) / ((PetscReal)ratioj);

388:         nc = 0;
389:         /* one left and below; or we are right on it */
390:         col      = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c));
391:         cols[nc] = col_shift + idx_c[col];
392:         v[nc++]  = x * y - x - y + 1.0;
393:         /* one right and below */
394:         if (i_c * ratioi != i) {
395:           cols[nc] = col_shift + idx_c[col + 1];
396:           v[nc++]  = -x * y + x;
397:         }
398:         /* one left and above */
399:         if (j_c * ratioj != j) {
400:           cols[nc] = col_shift + idx_c[col + m_ghost_c];
401:           v[nc++]  = -x * y + y;
402:         }
403:         /* one right and above */
404:         if (j_c * ratioj != j && i_c * ratioi != i) {
405:           cols[nc] = col_shift + idx_c[col + (m_ghost_c + 1)];
406:           v[nc++]  = x * y;
407:         }
408:         MatSetValues(mat, 1, &row, nc, cols, v, INSERT_VALUES);
409:       }
410:     }

412:   } else {
413:     PetscScalar  Ni[4];
414:     PetscScalar *xi, *eta;
415:     PetscInt     li, nxi, lj, neta;

417:     /* compute local coordinate arrays */
418:     nxi  = ratioi + 1;
419:     neta = ratioj + 1;
420:     PetscMalloc1(nxi, &xi);
421:     PetscMalloc1(neta, &eta);
422:     for (li = 0; li < nxi; li++) xi[li] = -1.0 + (PetscScalar)li * (2.0 / (PetscScalar)(nxi - 1));
423:     for (lj = 0; lj < neta; lj++) eta[lj] = -1.0 + (PetscScalar)lj * (2.0 / (PetscScalar)(neta - 1));

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

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

434:         /* remainders */
435:         li = i - ratioi * (i / ratioi);
436:         if (i == mx - 1) li = nxi - 1;
437:         lj = j - ratioj * (j / ratioj);
438:         if (j == my - 1) lj = neta - 1;

440:         /* corners */
441:         col     = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c));
442:         cols[0] = col_shift + idx_c[col]; /* left, below */
443:         Ni[0]   = 1.0;
444:         if ((li == 0) || (li == nxi - 1)) {
445:           if ((lj == 0) || (lj == neta - 1)) {
446:             MatSetValue(mat, row, cols[0], Ni[0], INSERT_VALUES);
447:             continue;
448:           }
449:         }

451:         /* edges + interior */
452:         /* remainders */
453:         if (i == mx - 1) i_c--;
454:         if (j == my - 1) j_c--;

456:         col     = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c));
457:         cols[0] = col_shift + idx_c[col];                   /* left, below */
458:         cols[1] = col_shift + idx_c[col + 1];               /* right, below */
459:         cols[2] = col_shift + idx_c[col + m_ghost_c];       /* left, above */
460:         cols[3] = col_shift + idx_c[col + (m_ghost_c + 1)]; /* right, above */

462:         Ni[0] = 0.25 * (1.0 - xi[li]) * (1.0 - eta[lj]);
463:         Ni[1] = 0.25 * (1.0 + xi[li]) * (1.0 - eta[lj]);
464:         Ni[2] = 0.25 * (1.0 - xi[li]) * (1.0 + eta[lj]);
465:         Ni[3] = 0.25 * (1.0 + xi[li]) * (1.0 + eta[lj]);

467:         nc = 0;
468:         if (PetscAbsScalar(Ni[0]) < 1.0e-32) cols[0] = -1;
469:         if (PetscAbsScalar(Ni[1]) < 1.0e-32) cols[1] = -1;
470:         if (PetscAbsScalar(Ni[2]) < 1.0e-32) cols[2] = -1;
471:         if (PetscAbsScalar(Ni[3]) < 1.0e-32) cols[3] = -1;

473:         MatSetValues(mat, 1, &row, 4, cols, Ni, INSERT_VALUES);
474:       }
475:     }
476:     PetscFree(xi);
477:     PetscFree(eta);
478:   }
479:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f);
480:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c);
481:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
482:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
483:   MatCreateMAIJ(mat, dof, A);
484:   MatDestroy(&mat);
485:   return 0;
486: }

488: /*
489:        Contributed by Andrei Draganescu <aidraga@sandia.gov>
490: */
491: PetscErrorCode DMCreateInterpolation_DA_2D_Q0(DM dac, DM daf, Mat *A)
492: {
493:   PetscInt               i, j, i_start, j_start, m_f, n_f, Mx, My, dof;
494:   const PetscInt        *idx_c, *idx_f;
495:   ISLocalToGlobalMapping ltog_f, ltog_c;
496:   PetscInt               m_ghost, n_ghost, m_ghost_c, n_ghost_c, *dnz, *onz;
497:   PetscInt               row, col, i_start_ghost, j_start_ghost, cols[4], mx, m_c, my, nc, ratioi, ratioj;
498:   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;
499:   PetscMPIInt            size_c, size_f, rank_f;
500:   PetscScalar            v[4];
501:   Mat                    mat;
502:   DMBoundaryType         bx, by;
503:   MatType                mattype;

505:   DMDAGetInfo(dac, NULL, &Mx, &My, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL);
506:   DMDAGetInfo(daf, NULL, &mx, &my, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
509:   ratioi = mx / Mx;
510:   ratioj = my / My;

516:   DMDAGetCorners(daf, &i_start, &j_start, NULL, &m_f, &n_f, NULL);
517:   DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, NULL, &m_ghost, &n_ghost, NULL);
518:   DMGetLocalToGlobalMapping(daf, &ltog_f);
519:   ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f);

521:   DMDAGetCorners(dac, &i_start_c, &j_start_c, NULL, &m_c, &n_c, NULL);
522:   DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, NULL, &m_ghost_c, &n_ghost_c, NULL);
523:   DMGetLocalToGlobalMapping(dac, &ltog_c);
524:   ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c);

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

534:      In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
535:   */
536:   MPI_Comm_size(PetscObjectComm((PetscObject)dac), &size_c);
537:   MPI_Comm_size(PetscObjectComm((PetscObject)daf), &size_f);
538:   MPI_Comm_rank(PetscObjectComm((PetscObject)daf), &rank_f);
539:   col_scale = size_f / size_c;
540:   col_shift = Mx * My * (rank_f / size_c);

542:   MatPreallocateBegin(PetscObjectComm((PetscObject)daf), m_f * n_f, col_scale * m_c * n_c, dnz, onz);
543:   for (j = j_start; j < j_start + n_f; j++) {
544:     for (i = i_start; i < i_start + m_f; i++) {
545:       /* convert to local "natural" numbering and then to PETSc global numbering */
546:       row = idx_f[(m_ghost * (j - j_start_ghost) + (i - i_start_ghost))];

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

552:     j_start %" PetscInt_FMT " j_c %" PetscInt_FMT " j_start_ghost_c %" PetscInt_FMT,
553:                  j_start, j_c, j_start_ghost_c);
555:     i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT,
556:                  i_start, i_c, i_start_ghost_c);

558:       /*
559:          Only include those interpolation points that are truly
560:          nonzero. Note this is very important for final grid lines
561:          in x and y directions; since they have no right/top neighbors
562:       */
563:       nc = 0;
564:       /* one left and below; or we are right on it */
565:       col        = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c));
566:       cols[nc++] = col_shift + idx_c[col];
567:       MatPreallocateSet(row, nc, cols, dnz, onz);
568:     }
569:   }
570:   MatCreate(PetscObjectComm((PetscObject)daf), &mat);
571: #if defined(PETSC_HAVE_CUDA)
572:   /*
573:      Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
574:      we don't want the original unconverted matrix copied to the GPU
575:   */
576:   if (dof > 1) MatBindToCPU(mat, PETSC_TRUE);
577: #endif
578:   MatSetSizes(mat, m_f * n_f, col_scale * m_c * n_c, mx * my, col_scale * Mx * My);
579:   ConvertToAIJ(dac->mattype, &mattype);
580:   MatSetType(mat, mattype);
581:   MatSeqAIJSetPreallocation(mat, 0, dnz);
582:   MatMPIAIJSetPreallocation(mat, 0, dnz, 0, onz);
583:   MatPreallocateEnd(dnz, onz);

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

591:       i_c = (i / ratioi); /* coarse grid node to left of fine grid node */
592:       j_c = (j / ratioj); /* coarse grid node below fine grid node */
593:       nc  = 0;
594:       /* one left and below; or we are right on it */
595:       col      = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c));
596:       cols[nc] = col_shift + idx_c[col];
597:       v[nc++]  = 1.0;

599:       MatSetValues(mat, 1, &row, nc, cols, v, INSERT_VALUES);
600:     }
601:   }
602:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f);
603:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c);
604:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
605:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
606:   MatCreateMAIJ(mat, dof, A);
607:   MatDestroy(&mat);
608:   PetscLogFlops(13.0 * m_f * n_f);
609:   return 0;
610: }

612: /*
613:        Contributed by Jianming Yang <jianming-yang@uiowa.edu>
614: */
615: PetscErrorCode DMCreateInterpolation_DA_3D_Q0(DM dac, DM daf, Mat *A)
616: {
617:   PetscInt               i, j, l, i_start, j_start, l_start, m_f, n_f, p_f, Mx, My, Mz, dof;
618:   const PetscInt        *idx_c, *idx_f;
619:   ISLocalToGlobalMapping ltog_f, ltog_c;
620:   PetscInt               m_ghost, n_ghost, p_ghost, m_ghost_c, n_ghost_c, p_ghost_c, nc, *dnz, *onz;
621:   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;
622:   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;
623:   PetscMPIInt            size_c, size_f, rank_f;
624:   PetscScalar            v[8];
625:   Mat                    mat;
626:   DMBoundaryType         bx, by, bz;
627:   MatType                mattype;

629:   DMDAGetInfo(dac, NULL, &Mx, &My, &Mz, NULL, NULL, NULL, NULL, NULL, &bx, &by, &bz, NULL);
633:   DMDAGetInfo(daf, NULL, &mx, &my, &mz, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
634:   ratioi = mx / Mx;
635:   ratioj = my / My;
636:   ratiol = mz / Mz;

644:   DMDAGetCorners(daf, &i_start, &j_start, &l_start, &m_f, &n_f, &p_f);
645:   DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, &l_start_ghost, &m_ghost, &n_ghost, &p_ghost);
646:   DMGetLocalToGlobalMapping(daf, &ltog_f);
647:   ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f);

649:   DMDAGetCorners(dac, &i_start_c, &j_start_c, &l_start_c, &m_c, &n_c, &p_c);
650:   DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, &l_start_ghost_c, &m_ghost_c, &n_ghost_c, &p_ghost_c);
651:   DMGetLocalToGlobalMapping(dac, &ltog_c);
652:   ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c);

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

662:      In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
663:   */
664:   MPI_Comm_size(PetscObjectComm((PetscObject)dac), &size_c);
665:   MPI_Comm_size(PetscObjectComm((PetscObject)daf), &size_f);
666:   MPI_Comm_rank(PetscObjectComm((PetscObject)daf), &rank_f);
667:   col_scale = size_f / size_c;
668:   col_shift = Mx * My * Mz * (rank_f / size_c);

670:   MatPreallocateBegin(PetscObjectComm((PetscObject)daf), m_f * n_f * p_f, col_scale * m_c * n_c * p_c, dnz, onz);
671:   for (l = l_start; l < l_start + p_f; l++) {
672:     for (j = j_start; j < j_start + n_f; j++) {
673:       for (i = i_start; i < i_start + m_f; i++) {
674:         /* convert to local "natural" numbering and then to PETSc global numbering */
675:         row = idx_f[(m_ghost * n_ghost * (l - l_start_ghost) + m_ghost * (j - j_start_ghost) + (i - i_start_ghost))];

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

682:     l_start %" PetscInt_FMT " l_c %" PetscInt_FMT " l_start_ghost_c %" PetscInt_FMT,
683:                    l_start, l_c, l_start_ghost_c);
685:     j_start %" PetscInt_FMT " j_c %" PetscInt_FMT " j_start_ghost_c %" PetscInt_FMT,
686:                    j_start, j_c, j_start_ghost_c);
688:     i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT,
689:                    i_start, i_c, i_start_ghost_c);

691:         /*
692:            Only include those interpolation points that are truly
693:            nonzero. Note this is very important for final grid lines
694:            in x and y directions; since they have no right/top neighbors
695:         */
696:         nc = 0;
697:         /* one left and below; or we are right on it */
698:         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));
699:         cols[nc++] = col_shift + idx_c[col];
700:         MatPreallocateSet(row, nc, cols, dnz, onz);
701:       }
702:     }
703:   }
704:   MatCreate(PetscObjectComm((PetscObject)daf), &mat);
705: #if defined(PETSC_HAVE_CUDA)
706:   /*
707:      Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
708:      we don't want the original unconverted matrix copied to the GPU
709:   */
710:   if (dof > 1) MatBindToCPU(mat, PETSC_TRUE);
711: #endif
712:   MatSetSizes(mat, m_f * n_f * p_f, col_scale * m_c * n_c * p_c, mx * my * mz, col_scale * Mx * My * Mz);
713:   ConvertToAIJ(dac->mattype, &mattype);
714:   MatSetType(mat, mattype);
715:   MatSeqAIJSetPreallocation(mat, 0, dnz);
716:   MatMPIAIJSetPreallocation(mat, 0, dnz, 0, onz);
717:   MatPreallocateEnd(dnz, onz);

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

726:         i_c = (i / ratioi); /* coarse grid node to left of fine grid node */
727:         j_c = (j / ratioj); /* coarse grid node below fine grid node */
728:         l_c = (l / ratiol);
729:         nc  = 0;
730:         /* one left and below; or we are right on it */
731:         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));
732:         cols[nc] = col_shift + idx_c[col];
733:         v[nc++]  = 1.0;

735:         MatSetValues(mat, 1, &row, nc, cols, v, INSERT_VALUES);
736:       }
737:     }
738:   }
739:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f);
740:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c);
741:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
742:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
743:   MatCreateMAIJ(mat, dof, A);
744:   MatDestroy(&mat);
745:   PetscLogFlops(13.0 * m_f * n_f * p_f);
746:   return 0;
747: }

749: PetscErrorCode DMCreateInterpolation_DA_3D_Q1(DM dac, DM daf, Mat *A)
750: {
751:   PetscInt               i, j, i_start, j_start, m_f, n_f, Mx, My, dof, l;
752:   const PetscInt        *idx_c, *idx_f;
753:   ISLocalToGlobalMapping ltog_f, ltog_c;
754:   PetscInt               m_ghost, n_ghost, m_ghost_c, n_ghost_c, Mz, mz;
755:   PetscInt               row, col, i_start_ghost, j_start_ghost, cols[8], mx, m_c, my, nc, ratioi, ratioj, ratiok;
756:   PetscInt               i_c, j_c, i_start_c, j_start_c, n_c, i_start_ghost_c, j_start_ghost_c;
757:   PetscInt               l_start, p_f, l_start_ghost, p_ghost, l_start_c, p_c;
758:   PetscInt               l_start_ghost_c, p_ghost_c, l_c, *dnz, *onz;
759:   PetscScalar            v[8], x, y, z;
760:   Mat                    mat;
761:   DMBoundaryType         bx, by, bz;
762:   MatType                mattype;

764:   DMDAGetInfo(dac, NULL, &Mx, &My, &Mz, NULL, NULL, NULL, NULL, NULL, &bx, &by, &bz, NULL);
765:   DMDAGetInfo(daf, NULL, &mx, &my, &mz, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
766:   if (mx == Mx) {
767:     ratioi = 1;
768:   } else if (bx == DM_BOUNDARY_PERIODIC) {
770:     ratioi = mx / Mx;
772:   } else {
774:     ratioi = (mx - 1) / (Mx - 1);
776:   }
777:   if (my == My) {
778:     ratioj = 1;
779:   } else if (by == DM_BOUNDARY_PERIODIC) {
781:     ratioj = my / My;
783:   } else {
785:     ratioj = (my - 1) / (My - 1);
787:   }
788:   if (mz == Mz) {
789:     ratiok = 1;
790:   } else if (bz == DM_BOUNDARY_PERIODIC) {
792:     ratiok = mz / Mz;
794:   } else {
796:     ratiok = (mz - 1) / (Mz - 1);
798:   }

800:   DMDAGetCorners(daf, &i_start, &j_start, &l_start, &m_f, &n_f, &p_f);
801:   DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, &l_start_ghost, &m_ghost, &n_ghost, &p_ghost);
802:   DMGetLocalToGlobalMapping(daf, &ltog_f);
803:   ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f);

805:   DMDAGetCorners(dac, &i_start_c, &j_start_c, &l_start_c, &m_c, &n_c, &p_c);
806:   DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, &l_start_ghost_c, &m_ghost_c, &n_ghost_c, &p_ghost_c);
807:   DMGetLocalToGlobalMapping(dac, &ltog_c);
808:   ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c);

810:   /* create interpolation matrix, determining exact preallocation */
811:   MatPreallocateBegin(PetscObjectComm((PetscObject)dac), m_f * n_f * p_f, m_c * n_c * p_c, dnz, onz);
812:   /* loop over local fine grid nodes counting interpolating points */
813:   for (l = l_start; l < l_start + p_f; l++) {
814:     for (j = j_start; j < j_start + n_f; j++) {
815:       for (i = i_start; i < i_start + m_f; i++) {
816:         /* convert to local "natural" numbering and then to PETSc global numbering */
817:         row = idx_f[(m_ghost * n_ghost * (l - l_start_ghost) + m_ghost * (j - j_start_ghost) + (i - i_start_ghost))];
818:         i_c = (i / ratioi);
819:         j_c = (j / ratioj);
820:         l_c = (l / ratiok);
822:                                             l_start %" PetscInt_FMT " l_c %" PetscInt_FMT " l_start_ghost_c %" PetscInt_FMT,
823:                    l_start, l_c, l_start_ghost_c);
825:                                             j_start %" PetscInt_FMT " j_c %" PetscInt_FMT " j_start_ghost_c %" PetscInt_FMT,
826:                    j_start, j_c, j_start_ghost_c);
828:                                             i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT,
829:                    i_start, i_c, i_start_ghost_c);

831:         /*
832:          Only include those interpolation points that are truly
833:          nonzero. Note this is very important for final grid lines
834:          in x and y directions; since they have no right/top neighbors
835:          */
836:         nc         = 0;
837:         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));
838:         cols[nc++] = idx_c[col];
839:         if (i_c * ratioi != i) cols[nc++] = idx_c[col + 1];
840:         if (j_c * ratioj != j) cols[nc++] = idx_c[col + m_ghost_c];
841:         if (l_c * ratiok != l) cols[nc++] = idx_c[col + m_ghost_c * n_ghost_c];
842:         if (j_c * ratioj != j && i_c * ratioi != i) cols[nc++] = idx_c[col + (m_ghost_c + 1)];
843:         if (j_c * ratioj != j && l_c * ratiok != l) cols[nc++] = idx_c[col + (m_ghost_c * n_ghost_c + m_ghost_c)];
844:         if (i_c * ratioi != i && l_c * ratiok != l) cols[nc++] = idx_c[col + (m_ghost_c * n_ghost_c + 1)];
845:         if (i_c * ratioi != i && l_c * ratiok != l && j_c * ratioj != j) cols[nc++] = idx_c[col + (m_ghost_c * n_ghost_c + m_ghost_c + 1)];
846:         MatPreallocateSet(row, nc, cols, dnz, onz);
847:       }
848:     }
849:   }
850:   MatCreate(PetscObjectComm((PetscObject)dac), &mat);
851: #if defined(PETSC_HAVE_CUDA)
852:   /*
853:      Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU
854:      we don't want the original unconverted matrix copied to the GPU
855:   */
856:   if (dof > 1) MatBindToCPU(mat, PETSC_TRUE);
857: #endif
858:   MatSetSizes(mat, m_f * n_f * p_f, m_c * n_c * p_c, mx * my * mz, Mx * My * Mz);
859:   ConvertToAIJ(dac->mattype, &mattype);
860:   MatSetType(mat, mattype);
861:   MatSeqAIJSetPreallocation(mat, 0, dnz);
862:   MatMPIAIJSetPreallocation(mat, 0, dnz, 0, onz);
863:   MatPreallocateEnd(dnz, onz);

865:   /* loop over local fine grid nodes setting interpolation for those*/
866:   if (!NEWVERSION) {
867:     for (l = l_start; l < l_start + p_f; l++) {
868:       for (j = j_start; j < j_start + n_f; j++) {
869:         for (i = i_start; i < i_start + m_f; i++) {
870:           /* convert to local "natural" numbering and then to PETSc global numbering */
871:           row = idx_f[(m_ghost * n_ghost * (l - l_start_ghost) + m_ghost * (j - j_start_ghost) + (i - i_start_ghost))];

873:           i_c = (i / ratioi);
874:           j_c = (j / ratioj);
875:           l_c = (l / ratiok);

877:           /*
878:            Only include those interpolation points that are truly
879:            nonzero. Note this is very important for final grid lines
880:            in x and y directions; since they have no right/top neighbors
881:            */
882:           x = ((PetscReal)(i - i_c * ratioi)) / ((PetscReal)ratioi);
883:           y = ((PetscReal)(j - j_c * ratioj)) / ((PetscReal)ratioj);
884:           z = ((PetscReal)(l - l_c * ratiok)) / ((PetscReal)ratiok);

886:           nc = 0;
887:           /* one left and below; or we are right on it */
888:           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));

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

893:           if (i_c * ratioi != i) {
894:             cols[nc] = idx_c[col + 1];
895:             v[nc++]  = .125 * (1. + (2.0 * x - 1.)) * (1. - (2.0 * y - 1.)) * (1. - (2.0 * z - 1.));
896:           }

898:           if (j_c * ratioj != j) {
899:             cols[nc] = idx_c[col + m_ghost_c];
900:             v[nc++]  = .125 * (1. - (2.0 * x - 1.)) * (1. + (2.0 * y - 1.)) * (1. - (2.0 * z - 1.));
901:           }

903:           if (l_c * ratiok != l) {
904:             cols[nc] = idx_c[col + m_ghost_c * n_ghost_c];
905:             v[nc++]  = .125 * (1. - (2.0 * x - 1.)) * (1. - (2.0 * y - 1.)) * (1. + (2.0 * z - 1.));
906:           }

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

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

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

923:           if (i_c * ratioi != i && l_c * ratiok != l && j_c * ratioj != j) {
924:             cols[nc] = idx_c[col + (m_ghost_c * n_ghost_c + m_ghost_c + 1)];
925:             v[nc++]  = .125 * (1. + (2.0 * x - 1.)) * (1. + (2.0 * y - 1.)) * (1. + (2.0 * z - 1.));
926:           }
927:           MatSetValues(mat, 1, &row, nc, cols, v, INSERT_VALUES);
928:         }
929:       }
930:     }

932:   } else {
933:     PetscScalar *xi, *eta, *zeta;
934:     PetscInt     li, nxi, lj, neta, lk, nzeta, n;
935:     PetscScalar  Ni[8];

937:     /* compute local coordinate arrays */
938:     nxi   = ratioi + 1;
939:     neta  = ratioj + 1;
940:     nzeta = ratiok + 1;
941:     PetscMalloc1(nxi, &xi);
942:     PetscMalloc1(neta, &eta);
943:     PetscMalloc1(nzeta, &zeta);
944:     for (li = 0; li < nxi; li++) xi[li] = -1.0 + (PetscScalar)li * (2.0 / (PetscScalar)(nxi - 1));
945:     for (lj = 0; lj < neta; lj++) eta[lj] = -1.0 + (PetscScalar)lj * (2.0 / (PetscScalar)(neta - 1));
946:     for (lk = 0; lk < nzeta; lk++) zeta[lk] = -1.0 + (PetscScalar)lk * (2.0 / (PetscScalar)(nzeta - 1));

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

954:           i_c = (i / ratioi);
955:           j_c = (j / ratioj);
956:           l_c = (l / ratiok);

958:           /* remainders */
959:           li = i - ratioi * (i / ratioi);
960:           if (i == mx - 1) li = nxi - 1;
961:           lj = j - ratioj * (j / ratioj);
962:           if (j == my - 1) lj = neta - 1;
963:           lk = l - ratiok * (l / ratiok);
964:           if (l == mz - 1) lk = nzeta - 1;

966:           /* corners */
967:           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));
968:           cols[0] = idx_c[col];
969:           Ni[0]   = 1.0;
970:           if ((li == 0) || (li == nxi - 1)) {
971:             if ((lj == 0) || (lj == neta - 1)) {
972:               if ((lk == 0) || (lk == nzeta - 1)) {
973:                 MatSetValue(mat, row, cols[0], Ni[0], INSERT_VALUES);
974:                 continue;
975:               }
976:             }
977:           }

979:           /* edges + interior */
980:           /* remainders */
981:           if (i == mx - 1) i_c--;
982:           if (j == my - 1) j_c--;
983:           if (l == mz - 1) l_c--;

985:           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));
986:           cols[0] = idx_c[col];                   /* one left and below; or we are right on it */
987:           cols[1] = idx_c[col + 1];               /* one right and below */
988:           cols[2] = idx_c[col + m_ghost_c];       /* one left and above */
989:           cols[3] = idx_c[col + (m_ghost_c + 1)]; /* one right and above */

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

996:           Ni[0] = 0.125 * (1.0 - xi[li]) * (1.0 - eta[lj]) * (1.0 - zeta[lk]);
997:           Ni[1] = 0.125 * (1.0 + xi[li]) * (1.0 - eta[lj]) * (1.0 - zeta[lk]);
998:           Ni[2] = 0.125 * (1.0 - xi[li]) * (1.0 + eta[lj]) * (1.0 - zeta[lk]);
999:           Ni[3] = 0.125 * (1.0 + xi[li]) * (1.0 + eta[lj]) * (1.0 - zeta[lk]);

1001:           Ni[4] = 0.125 * (1.0 - xi[li]) * (1.0 - eta[lj]) * (1.0 + zeta[lk]);
1002:           Ni[5] = 0.125 * (1.0 + xi[li]) * (1.0 - eta[lj]) * (1.0 + zeta[lk]);
1003:           Ni[6] = 0.125 * (1.0 - xi[li]) * (1.0 + eta[lj]) * (1.0 + zeta[lk]);
1004:           Ni[7] = 0.125 * (1.0 + xi[li]) * (1.0 + eta[lj]) * (1.0 + zeta[lk]);

1006:           for (n = 0; n < 8; n++) {
1007:             if (PetscAbsScalar(Ni[n]) < 1.0e-32) cols[n] = -1;
1008:           }
1009:           MatSetValues(mat, 1, &row, 8, cols, Ni, INSERT_VALUES);
1010:         }
1011:       }
1012:     }
1013:     PetscFree(xi);
1014:     PetscFree(eta);
1015:     PetscFree(zeta);
1016:   }
1017:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f);
1018:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c);
1019:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
1020:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);

1022:   MatCreateMAIJ(mat, dof, A);
1023:   MatDestroy(&mat);
1024:   return 0;
1025: }

1027: PetscErrorCode DMCreateInterpolation_DA(DM dac, DM daf, Mat *A, Vec *scale)
1028: {
1029:   PetscInt        dimc, Mc, Nc, Pc, mc, nc, pc, dofc, sc, dimf, Mf, Nf, Pf, mf, nf, pf, doff, sf;
1030:   DMBoundaryType  bxc, byc, bzc, bxf, byf, bzf;
1031:   DMDAStencilType stc, stf;
1032:   DM_DA          *ddc = (DM_DA *)dac->data;


1039:   DMDAGetInfo(dac, &dimc, &Mc, &Nc, &Pc, &mc, &nc, &pc, &dofc, &sc, &bxc, &byc, &bzc, &stc);
1040:   DMDAGetInfo(daf, &dimf, &Mf, &Nf, &Pf, &mf, &nf, &pf, &doff, &sf, &bxf, &byf, &bzf, &stf);

1050:   if (ddc->interptype == DMDA_Q1) {
1051:     if (dimc == 1) {
1052:       DMCreateInterpolation_DA_1D_Q1(dac, daf, A);
1053:     } else if (dimc == 2) {
1054:       DMCreateInterpolation_DA_2D_Q1(dac, daf, A);
1055:     } else if (dimc == 3) {
1056:       DMCreateInterpolation_DA_3D_Q1(dac, daf, A);
1057:     } else SETERRQ(PetscObjectComm((PetscObject)daf), PETSC_ERR_SUP, "No support for this DMDA dimension %" PetscInt_FMT " for interpolation type %d", dimc, (int)ddc->interptype);
1058:   } else if (ddc->interptype == DMDA_Q0) {
1059:     if (dimc == 1) {
1060:       DMCreateInterpolation_DA_1D_Q0(dac, daf, A);
1061:     } else if (dimc == 2) {
1062:       DMCreateInterpolation_DA_2D_Q0(dac, daf, A);
1063:     } else if (dimc == 3) {
1064:       DMCreateInterpolation_DA_3D_Q0(dac, daf, A);
1065:     } else SETERRQ(PetscObjectComm((PetscObject)daf), PETSC_ERR_SUP, "No support for this DMDA dimension %" PetscInt_FMT " for interpolation type %d", dimc, (int)ddc->interptype);
1066:   }
1067:   if (scale) DMCreateInterpolationScale((DM)dac, (DM)daf, *A, scale);
1068:   return 0;
1069: }

1071: PetscErrorCode DMCreateInjection_DA_1D(DM dac, DM daf, VecScatter *inject)
1072: {
1073:   PetscInt               i, i_start, m_f, Mx, dof;
1074:   const PetscInt        *idx_f;
1075:   ISLocalToGlobalMapping ltog_f;
1076:   PetscInt               m_ghost, m_ghost_c;
1077:   PetscInt               row, i_start_ghost, mx, m_c, nc, ratioi;
1078:   PetscInt               i_start_c, i_start_ghost_c;
1079:   PetscInt              *cols;
1080:   DMBoundaryType         bx;
1081:   Vec                    vecf, vecc;
1082:   IS                     isf;

1084:   DMDAGetInfo(dac, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bx, NULL, NULL, NULL);
1085:   DMDAGetInfo(daf, NULL, &mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
1086:   if (bx == DM_BOUNDARY_PERIODIC) {
1087:     ratioi = mx / Mx;
1089:   } else {
1090:     ratioi = (mx - 1) / (Mx - 1);
1092:   }

1094:   DMDAGetCorners(daf, &i_start, NULL, NULL, &m_f, NULL, NULL);
1095:   DMDAGetGhostCorners(daf, &i_start_ghost, NULL, NULL, &m_ghost, NULL, NULL);
1096:   DMGetLocalToGlobalMapping(daf, &ltog_f);
1097:   ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f);

1099:   DMDAGetCorners(dac, &i_start_c, NULL, NULL, &m_c, NULL, NULL);
1100:   DMDAGetGhostCorners(dac, &i_start_ghost_c, NULL, NULL, &m_ghost_c, NULL, NULL);

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

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


1111:     row        = idx_f[(i_f - i_start_ghost)];
1112:     cols[nc++] = row;
1113:   }

1115:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f);
1116:   ISCreateBlock(PetscObjectComm((PetscObject)daf), dof, nc, cols, PETSC_OWN_POINTER, &isf);
1117:   DMGetGlobalVector(dac, &vecc);
1118:   DMGetGlobalVector(daf, &vecf);
1119:   VecScatterCreate(vecf, isf, vecc, NULL, inject);
1120:   DMRestoreGlobalVector(dac, &vecc);
1121:   DMRestoreGlobalVector(daf, &vecf);
1122:   ISDestroy(&isf);
1123:   return 0;
1124: }

1126: PetscErrorCode DMCreateInjection_DA_2D(DM dac, DM daf, VecScatter *inject)
1127: {
1128:   PetscInt               i, j, i_start, j_start, m_f, n_f, Mx, My, dof;
1129:   const PetscInt        *idx_c, *idx_f;
1130:   ISLocalToGlobalMapping ltog_f, ltog_c;
1131:   PetscInt               m_ghost, n_ghost, m_ghost_c, n_ghost_c;
1132:   PetscInt               row, i_start_ghost, j_start_ghost, mx, m_c, my, nc, ratioi, ratioj;
1133:   PetscInt               i_start_c, j_start_c, n_c, i_start_ghost_c, j_start_ghost_c;
1134:   PetscInt              *cols;
1135:   DMBoundaryType         bx, by;
1136:   Vec                    vecf, vecc;
1137:   IS                     isf;

1139:   DMDAGetInfo(dac, NULL, &Mx, &My, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL);
1140:   DMDAGetInfo(daf, NULL, &mx, &my, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
1141:   if (bx == DM_BOUNDARY_PERIODIC) {
1142:     ratioi = mx / Mx;
1144:   } else {
1145:     ratioi = (mx - 1) / (Mx - 1);
1147:   }
1148:   if (by == DM_BOUNDARY_PERIODIC) {
1149:     ratioj = my / My;
1151:   } else {
1152:     ratioj = (my - 1) / (My - 1);
1154:   }

1156:   DMDAGetCorners(daf, &i_start, &j_start, NULL, &m_f, &n_f, NULL);
1157:   DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, NULL, &m_ghost, &n_ghost, NULL);
1158:   DMGetLocalToGlobalMapping(daf, &ltog_f);
1159:   ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f);

1161:   DMDAGetCorners(dac, &i_start_c, &j_start_c, NULL, &m_c, &n_c, NULL);
1162:   DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, NULL, &m_ghost_c, &n_ghost_c, NULL);
1163:   DMGetLocalToGlobalMapping(dac, &ltog_c);
1164:   ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c);

1166:   /* loop over local fine grid nodes setting interpolation for those*/
1167:   nc = 0;
1168:   PetscMalloc1(n_f * m_f, &cols);
1169:   for (j = j_start_c; j < j_start_c + n_c; j++) {
1170:     for (i = i_start_c; i < i_start_c + m_c; i++) {
1171:       PetscInt i_f = i * ratioi, j_f = j * ratioj;
1173:     j_c %" PetscInt_FMT " j_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]",
1174:                  j, j_f, j_start_ghost, j_start_ghost + n_ghost);
1176:     i_c %" PetscInt_FMT " i_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]",
1177:                  i, i_f, i_start_ghost, i_start_ghost + m_ghost);
1178:       row        = idx_f[(m_ghost * (j_f - j_start_ghost) + (i_f - i_start_ghost))];
1179:       cols[nc++] = row;
1180:     }
1181:   }
1182:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f);
1183:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c);

1185:   ISCreateBlock(PetscObjectComm((PetscObject)daf), dof, nc, cols, PETSC_OWN_POINTER, &isf);
1186:   DMGetGlobalVector(dac, &vecc);
1187:   DMGetGlobalVector(daf, &vecf);
1188:   VecScatterCreate(vecf, isf, vecc, NULL, inject);
1189:   DMRestoreGlobalVector(dac, &vecc);
1190:   DMRestoreGlobalVector(daf, &vecf);
1191:   ISDestroy(&isf);
1192:   return 0;
1193: }

1195: PetscErrorCode DMCreateInjection_DA_3D(DM dac, DM daf, VecScatter *inject)
1196: {
1197:   PetscInt               i, j, k, i_start, j_start, k_start, m_f, n_f, p_f, Mx, My, Mz;
1198:   PetscInt               m_ghost, n_ghost, p_ghost, m_ghost_c, n_ghost_c, p_ghost_c;
1199:   PetscInt               i_start_ghost, j_start_ghost, k_start_ghost;
1200:   PetscInt               mx, my, mz, ratioi, ratioj, ratiok;
1201:   PetscInt               i_start_c, j_start_c, k_start_c;
1202:   PetscInt               m_c, n_c, p_c;
1203:   PetscInt               i_start_ghost_c, j_start_ghost_c, k_start_ghost_c;
1204:   PetscInt               row, nc, dof;
1205:   const PetscInt        *idx_c, *idx_f;
1206:   ISLocalToGlobalMapping ltog_f, ltog_c;
1207:   PetscInt              *cols;
1208:   DMBoundaryType         bx, by, bz;
1209:   Vec                    vecf, vecc;
1210:   IS                     isf;

1212:   DMDAGetInfo(dac, NULL, &Mx, &My, &Mz, NULL, NULL, NULL, NULL, NULL, &bx, &by, &bz, NULL);
1213:   DMDAGetInfo(daf, NULL, &mx, &my, &mz, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);

1215:   if (bx == DM_BOUNDARY_PERIODIC) {
1216:     ratioi = mx / Mx;
1218:   } else {
1219:     ratioi = (mx - 1) / (Mx - 1);
1221:   }
1222:   if (by == DM_BOUNDARY_PERIODIC) {
1223:     ratioj = my / My;
1225:   } else {
1226:     ratioj = (my - 1) / (My - 1);
1228:   }
1229:   if (bz == DM_BOUNDARY_PERIODIC) {
1230:     ratiok = mz / Mz;
1232:   } else {
1233:     ratiok = (mz - 1) / (Mz - 1);
1235:   }

1237:   DMDAGetCorners(daf, &i_start, &j_start, &k_start, &m_f, &n_f, &p_f);
1238:   DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, &k_start_ghost, &m_ghost, &n_ghost, &p_ghost);
1239:   DMGetLocalToGlobalMapping(daf, &ltog_f);
1240:   ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f);

1242:   DMDAGetCorners(dac, &i_start_c, &j_start_c, &k_start_c, &m_c, &n_c, &p_c);
1243:   DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, &k_start_ghost_c, &m_ghost_c, &n_ghost_c, &p_ghost_c);
1244:   DMGetLocalToGlobalMapping(dac, &ltog_c);
1245:   ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c);

1247:   /* loop over local fine grid nodes setting interpolation for those*/
1248:   nc = 0;
1249:   PetscMalloc1(n_f * m_f * p_f, &cols);
1250:   for (k = k_start_c; k < k_start_c + p_c; k++) {
1251:     for (j = j_start_c; j < j_start_c + n_c; j++) {
1252:       for (i = i_start_c; i < i_start_c + m_c; i++) {
1253:         PetscInt i_f = i * ratioi, j_f = j * ratioj, k_f = k * ratiok;
1255:                    "Processor's coarse DMDA must lie over fine DMDA  "
1256:                    "k_c %" PetscInt_FMT " k_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]",
1257:                    k, k_f, k_start_ghost, k_start_ghost + p_ghost);
1259:                    "Processor's coarse DMDA must lie over fine DMDA  "
1260:                    "j_c %" PetscInt_FMT " j_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]",
1261:                    j, j_f, j_start_ghost, j_start_ghost + n_ghost);
1263:                    "Processor's coarse DMDA must lie over fine DMDA  "
1264:                    "i_c %" PetscInt_FMT " i_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]",
1265:                    i, i_f, i_start_ghost, i_start_ghost + m_ghost);
1266:         row        = idx_f[(m_ghost * n_ghost * (k_f - k_start_ghost) + m_ghost * (j_f - j_start_ghost) + (i_f - i_start_ghost))];
1267:         cols[nc++] = row;
1268:       }
1269:     }
1270:   }
1271:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f);
1272:   ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c);

1274:   ISCreateBlock(PetscObjectComm((PetscObject)daf), dof, nc, cols, PETSC_OWN_POINTER, &isf);
1275:   DMGetGlobalVector(dac, &vecc);
1276:   DMGetGlobalVector(daf, &vecf);
1277:   VecScatterCreate(vecf, isf, vecc, NULL, inject);
1278:   DMRestoreGlobalVector(dac, &vecc);
1279:   DMRestoreGlobalVector(daf, &vecf);
1280:   ISDestroy(&isf);
1281:   return 0;
1282: }

1284: PetscErrorCode DMCreateInjection_DA(DM dac, DM daf, Mat *mat)
1285: {
1286:   PetscInt        dimc, Mc, Nc, Pc, mc, nc, pc, dofc, sc, dimf, Mf, Nf, Pf, mf, nf, pf, doff, sf;
1287:   DMBoundaryType  bxc, byc, bzc, bxf, byf, bzf;
1288:   DMDAStencilType stc, stf;
1289:   VecScatter      inject = NULL;


1295:   DMDAGetInfo(dac, &dimc, &Mc, &Nc, &Pc, &mc, &nc, &pc, &dofc, &sc, &bxc, &byc, &bzc, &stc);
1296:   DMDAGetInfo(daf, &dimf, &Mf, &Nf, &Pf, &mf, &nf, &pf, &doff, &sf, &bxf, &byf, &bzf, &stf);

1306:   if (dimc == 1) {
1307:     DMCreateInjection_DA_1D(dac, daf, &inject);
1308:   } else if (dimc == 2) {
1309:     DMCreateInjection_DA_2D(dac, daf, &inject);
1310:   } else if (dimc == 3) {
1311:     DMCreateInjection_DA_3D(dac, daf, &inject);
1312:   }
1313:   MatCreateScatter(PetscObjectComm((PetscObject)inject), inject, mat);
1314:   VecScatterDestroy(&inject);
1315:   return 0;
1316: }

1318: /*@
1319:    DMCreateAggregates - Deprecated, see DMDACreateAggregates.

1321:    Level: intermediate
1322: @*/
1323: PetscErrorCode DMCreateAggregates(DM dac, DM daf, Mat *mat)
1324: {
1325:   return DMDACreateAggregates(dac, daf, mat);
1326: }

1328: /*@
1329:    DMDACreateAggregates - Gets the aggregates that map between
1330:    grids associated with two `DMDA`

1332:    Collective on dmc

1334:    Input Parameters:
1335: +  dmc - the coarse grid `DMDA`
1336: -  dmf - the fine grid `DMDA`

1338:    Output Parameters:
1339: .  rest - the restriction matrix (transpose of the projection matrix)

1341:    Level: intermediate

1343:    Note:
1344:    This routine is not used by PETSc.
1345:    It is not clear what its use case is and it may be removed in a future release.
1346:    Users should contact petsc-maint@mcs.anl.gov if they plan to use it.

1348: .seealso: `DMRefine()`, `DMCreateInjection()`, `DMCreateInterpolation()`
1349: @*/
1350: PetscErrorCode DMDACreateAggregates(DM dac, DM daf, Mat *rest)
1351: {
1352:   PetscInt               dimc, Mc, Nc, Pc, mc, nc, pc, dofc, sc;
1353:   PetscInt               dimf, Mf, Nf, Pf, mf, nf, pf, doff, sf;
1354:   DMBoundaryType         bxc, byc, bzc, bxf, byf, bzf;
1355:   DMDAStencilType        stc, stf;
1356:   PetscInt               i, j, l;
1357:   PetscInt               i_start, j_start, l_start, m_f, n_f, p_f;
1358:   PetscInt               i_start_ghost, j_start_ghost, l_start_ghost, m_ghost, n_ghost, p_ghost;
1359:   const PetscInt        *idx_f;
1360:   PetscInt               i_c, j_c, l_c;
1361:   PetscInt               i_start_c, j_start_c, l_start_c, m_c, n_c, p_c;
1362:   PetscInt               i_start_ghost_c, j_start_ghost_c, l_start_ghost_c, m_ghost_c, n_ghost_c, p_ghost_c;
1363:   const PetscInt        *idx_c;
1364:   PetscInt               d;
1365:   PetscInt               a;
1366:   PetscInt               max_agg_size;
1367:   PetscInt              *fine_nodes;
1368:   PetscScalar           *one_vec;
1369:   PetscInt               fn_idx;
1370:   ISLocalToGlobalMapping ltogmf, ltogmc;


1376:   DMDAGetInfo(dac, &dimc, &Mc, &Nc, &Pc, &mc, &nc, &pc, &dofc, &sc, &bxc, &byc, &bzc, &stc);
1377:   DMDAGetInfo(daf, &dimf, &Mf, &Nf, &Pf, &mf, &nf, &pf, &doff, &sf, &bxf, &byf, &bzf, &stf);


1388:   if (Pc < 0) Pc = 1;
1389:   if (Pf < 0) Pf = 1;
1390:   if (Nc < 0) Nc = 1;
1391:   if (Nf < 0) Nf = 1;

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

1396:   DMGetLocalToGlobalMapping(daf, &ltogmf);
1397:   ISLocalToGlobalMappingGetIndices(ltogmf, &idx_f);

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

1402:   DMGetLocalToGlobalMapping(dac, &ltogmc);
1403:   ISLocalToGlobalMappingGetIndices(ltogmc, &idx_c);

1405:   /*
1406:      Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios
1407:      for dimension 1 and 2 respectively.
1408:      Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1)
1409:      and r_y*j and r_y*(j+1) will be grouped into the same coarse grid aggregate.
1410:      Each specific dof on the fine grid is mapped to one dof on the coarse grid.
1411:   */

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

1415:   /* create the matrix that will contain the restriction operator */
1416:   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, max_agg_size, NULL, max_agg_size, NULL, rest);

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

1422:   /* loop over all coarse nodes */
1423:   for (l_c = l_start_c; l_c < l_start_c + p_c; l_c++) {
1424:     for (j_c = j_start_c; j_c < j_start_c + n_c; j_c++) {
1425:       for (i_c = i_start_c; i_c < i_start_c + m_c; i_c++) {
1426:         for (d = 0; d < dofc; d++) {
1427:           /* convert to local "natural" numbering and then to PETSc global numbering */
1428:           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;

1430:           fn_idx = 0;
1431:           /* Corresponding fine points are all points (i_f, j_f, l_f) such that
1432:              i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc
1433:              (same for other dimensions)
1434:           */
1435:           for (l = l_c * Pf / Pc; l < PetscMin((l_c + 1) * Pf / Pc, Pf); l++) {
1436:             for (j = j_c * Nf / Nc; j < PetscMin((j_c + 1) * Nf / Nc, Nf); j++) {
1437:               for (i = i_c * Mf / Mc; i < PetscMin((i_c + 1) * Mf / Mc, Mf); i++) {
1438:                 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;
1439:                 fn_idx++;
1440:               }
1441:             }
1442:           }
1443:           /* add all these points to one aggregate */
1444:           MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);
1445:         }
1446:       }
1447:     }
1448:   }
1449:   ISLocalToGlobalMappingRestoreIndices(ltogmf, &idx_f);
1450:   ISLocalToGlobalMappingRestoreIndices(ltogmc, &idx_c);
1451:   PetscFree2(one_vec, fine_nodes);
1452:   MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);
1453:   MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);
1454:   return 0;
1455: }