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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <og_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, <ogmf);
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, <ogmc);
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: }