Actual source code: dainterp.c

      1: #define PETSCDM_DLL

      3: /*
      4:   Code for interpolating between grids represented by DAs
      5: */

     7:  #include src/dm/da/daimpl.h
     8:  #include petscmg.h

     12: PetscErrorCode PETSCDM_DLLEXPORT DMGetInterpolationScale(DM dac,DM daf,Mat mat,Vec *scale)
     13: {
     15:   Vec            fine;
     16:   PetscScalar    one = 1.0;

     19:   DMCreateGlobalVector(daf,&fine);
     20:   DMCreateGlobalVector(dac,scale);
     21:   VecSet(fine,one);
     22:   MatRestrict(mat,fine,*scale);
     23:   VecDestroy(fine);
     24:   VecReciprocal(*scale);
     25:   return(0);
     26: }

     30: PetscErrorCode DAGetInterpolation_1D_Q1(DA dac,DA daf,Mat *A)
     31: {
     33:   PetscInt       i,i_start,m_f,Mx,*idx_f;
     34:   PetscInt       m_ghost,*idx_c,m_ghost_c;
     35:   PetscInt       row,col,i_start_ghost,mx,m_c,nc,ratio;
     36:   PetscInt       i_c,i_start_c,i_start_ghost_c,cols[2],dof;
     37:   PetscScalar    v[2],x,*coors = 0,*ccoors;
     38:   Mat            mat;
     39:   DAPeriodicType pt;
     40:   Vec            vcoors,cvcoors;

     43:   DAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&pt,0);
     44:   DAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0);
     45:   if (pt == DA_XPERIODIC) {
     46:     ratio = mx/Mx;
     47:     if (ratio*Mx != mx) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
     48:   } else {
     49:     ratio = (mx-1)/(Mx-1);
     50:     if (ratio*(Mx-1) != mx-1) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
     51:   }

     53:   DAGetCorners(daf,&i_start,0,0,&m_f,0,0);
     54:   DAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
     55:   DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);

     57:   DAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
     58:   DAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
     59:   DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);

     61:   /* create interpolation matrix */
     62:   MatCreate(dac->comm,&mat);
     63:   MatSetSizes(mat,m_f,m_c,mx,Mx);
     64:   MatSetType(mat,MATAIJ);
     65:   MatSeqAIJSetPreallocation(mat,2,PETSC_NULL);
     66:   MatMPIAIJSetPreallocation(mat,2,PETSC_NULL,0,PETSC_NULL);
     67:   if (!DAXPeriodic(pt)){MatSetOption(mat,MAT_COLUMNS_SORTED);}

     69:   DAGetCoordinates(daf,&vcoors);
     70:   if (vcoors) {
     71:     DAGetGhostedCoordinates(dac,&cvcoors);
     72:     DAVecGetArray(daf->da_coordinates,vcoors,&coors);
     73:     DAVecGetArray(dac->da_coordinates,cvcoors,&ccoors);
     74:   }
     75:   /* loop over local fine grid nodes setting interpolation for those*/
     76:   for (i=i_start; i<i_start+m_f; i++) {
     77:     /* convert to local "natural" numbering and then to PETSc global numbering */
     78:     row    = idx_f[dof*(i-i_start_ghost)]/dof;

     80:     i_c = (i/ratio);    /* coarse grid node to left of fine grid node */
     81:     if (i_c < i_start_ghost_c) SETERRQ3(PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\
     82:     i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);

     84:     /*
     85:          Only include those interpolation points that are truly
     86:          nonzero. Note this is very important for final grid lines
     87:          in x direction; since they have no right neighbor
     88:     */
     89:     if (coors) {
     90:       x = (coors[i] - ccoors[i_c]);
     91:       /* only access the next coors point if we know there is one */
     92:       /* note this is dangerous because x may not exactly equal ZERO */
     93:       if (PetscAbsScalar(x) != 0.0) x = x/(ccoors[i_c+1] - ccoors[i_c]);
     94:     } else {
     95:       x  = ((double)(i - i_c*ratio))/((double)ratio);
     96:     }
     97:     nc = 0;
     98:       /* one left and below; or we are right on it */
     99:     col      = dof*(i_c-i_start_ghost_c);
    100:     cols[nc] = idx_c[col]/dof;
    101:     v[nc++]  = - x + 1.0;
    102:     /* one right? */
    103:     if (i_c*ratio != i) {
    104:       cols[nc] = idx_c[col+dof]/dof;
    105:       v[nc++]  = x;
    106:     }
    107:     MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
    108:   }
    109:   if (vcoors) {
    110:     DAVecRestoreArray(daf->da_coordinates,vcoors,&coors);
    111:     DAVecRestoreArray(dac->da_coordinates,cvcoors,&ccoors);
    112:   }
    113:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
    114:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
    115:   MatCreateMAIJ(mat,dof,A);
    116:   MatDestroy(mat);
    117:   PetscLogFlops(5*m_f);
    118:   return(0);
    119: }

    123: PetscErrorCode DAGetInterpolation_1D_Q0(DA dac,DA daf,Mat *A)
    124: {
    126:   PetscInt       i,i_start,m_f,Mx,*idx_f;
    127:   PetscInt       m_ghost,*idx_c,m_ghost_c;
    128:   PetscInt       row,col,i_start_ghost,mx,m_c,nc,ratio;
    129:   PetscInt       i_c,i_start_c,i_start_ghost_c,cols[2],dof;
    130:   PetscScalar    v[2],x;
    131:   Mat            mat;
    132:   DAPeriodicType pt;

    135:   DAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&pt,0);
    136:   DAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0);
    137:   if (pt == DA_XPERIODIC) {
    138:     ratio = mx/Mx;
    139:     if (ratio*Mx != mx) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
    140:   } else {
    141:     ratio = (mx-1)/(Mx-1);
    142:     if (ratio*(Mx-1) != mx-1) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
    143:   }

    145:   DAGetCorners(daf,&i_start,0,0,&m_f,0,0);
    146:   DAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);
    147:   DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);

    149:   DAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);
    150:   DAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);
    151:   DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);

    153:   /* create interpolation matrix */
    154:   MatCreate(dac->comm,&mat);
    155:   MatSetSizes(mat,m_f,m_c,mx,Mx);
    156:   MatSetType(mat,MATAIJ);
    157:   MatSeqAIJSetPreallocation(mat,2,PETSC_NULL);
    158:   MatMPIAIJSetPreallocation(mat,2,PETSC_NULL,0,PETSC_NULL);
    159:   if (!DAXPeriodic(pt)) {MatSetOption(mat,MAT_COLUMNS_SORTED);}

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

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

    168:     /*
    169:          Only include those interpolation points that are truly
    170:          nonzero. Note this is very important for final grid lines
    171:          in x direction; since they have no right neighbor
    172:     */
    173:     x  = ((double)(i - i_c*ratio))/((double)ratio);
    174:     nc = 0;
    175:       /* one left and below; or we are right on it */
    176:     col      = dof*(i_c-i_start_ghost_c);
    177:     cols[nc] = idx_c[col]/dof;
    178:     v[nc++]  = - x + 1.0;
    179:     /* one right? */
    180:     if (i_c*ratio != i) {
    181:       cols[nc] = idx_c[col+dof]/dof;
    182:       v[nc++]  = x;
    183:     }
    184:     MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
    185:   }
    186:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
    187:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
    188:   MatCreateMAIJ(mat,dof,A);
    189:   MatDestroy(mat);
    190:   PetscLogFlops(5*m_f);
    191:   return(0);
    192: }

    196: PetscErrorCode DAGetInterpolation_2D_Q1(DA dac,DA daf,Mat *A)
    197: {
    199:   PetscInt       i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof;
    200:   PetscInt       m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,*dnz,*onz;
    201:   PetscInt       row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj;
    202:   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;
    203:   PetscMPIInt    size_c,size_f,rank_f;
    204:   PetscScalar    v[4],x,y;
    205:   Mat            mat;
    206:   DAPeriodicType pt;
    207:   DACoor2d       **coors = 0,**ccoors;
    208:   Vec            vcoors,cvcoors;

    211:   DAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&pt,0);
    212:   DAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0);
    213:   if (DAXPeriodic(pt)){
    214:     ratioi = mx/Mx;
    215:     if (ratioi*Mx != mx) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
    216:   } else {
    217:     ratioi = (mx-1)/(Mx-1);
    218:     if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
    219:   }
    220:   if (DAYPeriodic(pt)){
    221:     ratioj = my/My;
    222:     if (ratioj*My != my) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My  must be integer: my %D My %D",my,My);
    223:   } else {
    224:     ratioj = (my-1)/(My-1);
    225:     if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
    226:   }

    229:   DAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
    230:   DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
    231:   DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);

    233:   DAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
    234:   DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
    235:   DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);

    237:   /*
    238:      Used for handling a coarse DA that lives on 1/4 the processors of the fine DA.
    239:      The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the
    240:      processors). It's effective length is hence 4 times its normal length, this is
    241:      why the col_scale is multiplied by the interpolation matrix column sizes.
    242:      sol_shift allows each set of 1/4 processors do its own interpolation using ITS
    243:      copy of the coarse vector. A bit of a hack but you do better.

    245:      In the standard case when size_f == size_c col_scale == 1 and col_shift == 0
    246:   */
    247:   MPI_Comm_size(dac->comm,&size_c);
    248:   MPI_Comm_size(daf->comm,&size_f);
    249:   MPI_Comm_rank(daf->comm,&rank_f);
    250:   col_scale = size_f/size_c;
    251:   col_shift = Mx*My*(rank_f/size_c);

    253:   MatPreallocateInitialize(daf->comm,m_f*n_f,col_scale*m_c*n_c,dnz,onz);
    254:   for (j=j_start; j<j_start+n_f; j++) {
    255:     for (i=i_start; i<i_start+m_f; i++) {
    256:       /* convert to local "natural" numbering and then to PETSc global numbering */
    257:       row    = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;

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

    262:       if (j_c < j_start_ghost_c) SETERRQ3(PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\
    263:     j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
    264:       if (i_c < i_start_ghost_c) SETERRQ3(PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\
    265:     i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);

    267:       /*
    268:          Only include those interpolation points that are truly
    269:          nonzero. Note this is very important for final grid lines
    270:          in x and y directions; since they have no right/top neighbors
    271:       */
    272:       nc = 0;
    273:       /* one left and below; or we are right on it */
    274:       col        = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
    275:       cols[nc++] = col_shift + idx_c[col]/dof;
    276:       /* one right and below */
    277:       if (i_c*ratioi != i) {
    278:         cols[nc++] = col_shift + idx_c[col+dof]/dof;
    279:       }
    280:       /* one left and above */
    281:       if (j_c*ratioj != j) {
    282:         cols[nc++] = col_shift + idx_c[col+m_ghost_c*dof]/dof;
    283:       }
    284:       /* one right and above */
    285:       if (j_c*ratioi != j && i_c*ratioj != i) {
    286:         cols[nc++] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof;
    287:       }
    288:       MatPreallocateSet(row,nc,cols,dnz,onz);
    289:     }
    290:   }
    291:   MatCreate(daf->comm,&mat);
    292:   MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);
    293:   MatSetType(mat,MATAIJ);
    294:   MatSeqAIJSetPreallocation(mat,0,dnz);
    295:   MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
    296:   MatPreallocateFinalize(dnz,onz);
    297:   if (!DAXPeriodic(pt) && !DAYPeriodic(pt)) {MatSetOption(mat,MAT_COLUMNS_SORTED);}

    299:   DAGetCoordinates(daf,&vcoors);
    300:   if (vcoors) {
    301:     DAGetGhostedCoordinates(dac,&cvcoors);
    302:     DAVecGetArray(daf->da_coordinates,vcoors,&coors);
    303:     DAVecGetArray(dac->da_coordinates,cvcoors,&ccoors);
    304:   }

    306:   /* loop over local fine grid nodes setting interpolation for those*/
    307:   for (j=j_start; j<j_start+n_f; j++) {
    308:     for (i=i_start; i<i_start+m_f; i++) {
    309:       /* convert to local "natural" numbering and then to PETSc global numbering */
    310:       row    = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;

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

    315:       /*
    316:          Only include those interpolation points that are truly
    317:          nonzero. Note this is very important for final grid lines
    318:          in x and y directions; since they have no right/top neighbors
    319:       */
    320:       if (coors) {
    321:         /* only access the next coors point if we know there is one */
    322:         /* note this is dangerous because x may not exactly equal ZERO */
    323:         x = (coors[j][i].x - ccoors[j_c][i_c].x);
    324:         if (PetscAbsScalar(x) != 0.0) x = x/(ccoors[j_c][i_c+1].x - ccoors[j_c][i_c].x);
    325:         y = (coors[j][i].y - ccoors[j_c][i_c].y);
    326:         if (PetscAbsScalar(y) != 0.0) y = y/(ccoors[j_c+1][i_c].y - ccoors[j_c][i_c].y);
    327:       } else {
    328:         x  = ((double)(i - i_c*ratioi))/((double)ratioi);
    329:         y  = ((double)(j - j_c*ratioj))/((double)ratioj);
    330:       }
    331:       nc = 0;
    332:       /* one left and below; or we are right on it */
    333:       col      = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
    334:       cols[nc] = col_shift + idx_c[col]/dof;
    335:       v[nc++]  = x*y - x - y + 1.0;
    336:       /* one right and below */
    337:       if (i_c*ratioi != i) {
    338:         cols[nc] = col_shift + idx_c[col+dof]/dof;
    339:         v[nc++]  = -x*y + x;
    340:       }
    341:       /* one left and above */
    342:       if (j_c*ratioj != j) {
    343:         cols[nc] = col_shift + idx_c[col+m_ghost_c*dof]/dof;
    344:         v[nc++]  = -x*y + y;
    345:       }
    346:       /* one right and above */
    347:       if (j_c*ratioj != j && i_c*ratioi != i) {
    348:         cols[nc] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof;
    349:         v[nc++]  = x*y;
    350:       }
    351:       MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
    352:     }
    353:   }
    354:   if (vcoors) {
    355:     DAVecRestoreArray(daf->da_coordinates,vcoors,&coors);
    356:     DAVecRestoreArray(dac->da_coordinates,cvcoors,&ccoors);
    357:   }
    358:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
    359:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
    360:   MatCreateMAIJ(mat,dof,A);
    361:   MatDestroy(mat);
    362:   PetscLogFlops(13*m_f*n_f);
    363:   return(0);
    364: }

    367: /*   dof degree of freedom per node, nonperiodic */
    370: PetscErrorCode DAGetInterpolation_3D_Q1(DA dac,DA daf,Mat *A)
    371: {
    373:   PetscInt       i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof,l;
    374:   PetscInt       m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,Mz,mz;
    375:   PetscInt       row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok;
    376:   PetscInt       i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
    377:   PetscInt       l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c;
    378:   PetscInt       l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz;
    379:   PetscScalar    v[8],x,y,z;
    380:   Mat            mat;
    381:   DAPeriodicType pt;
    382:   DACoor3d       ***coors = 0,***ccoors;
    383:   Vec            vcoors,cvcoors;

    386:   DAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&pt,0);
    387:   DAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);
    388:   if (mx == Mx) {
    389:     ratioi = 1;
    390:   } else if (DAXPeriodic(pt)){
    391:     ratioi = mx/Mx;
    392:     if (ratioi*Mx != mx) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
    393:   } else {
    394:     ratioi = (mx-1)/(Mx-1);
    395:     if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
    396:   }
    397:   if (my == My) {
    398:     ratioj = 1;
    399:   } else if (DAYPeriodic(pt)){
    400:     ratioj = my/My;
    401:     if (ratioj*My != my) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My  must be integer: my %D My %D",my,My);
    402:   } else {
    403:     ratioj = (my-1)/(My-1);
    404:     if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
    405:   }
    406:   if (mz == Mz) {
    407:     ratiok = 1;
    408:   } else if (DAZPeriodic(pt)){
    409:     ratiok = mz/Mz;
    410:     if (ratiok*Mz != mz) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: mz/Mz  must be integer: mz %D Mz %D",mz,Mz);
    411:   } else {
    412:     ratiok = (mz-1)/(Mz-1);
    413:     if (ratiok*(Mz-1) != mz-1) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %D Mz %D",mz,Mz);
    414:   }

    416:   DAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);
    417:   DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);
    418:   DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);

    420:   DAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);
    421:   DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);
    422:   DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);

    424:   /* create interpolation matrix, determining exact preallocation */
    425:   MatPreallocateInitialize(dac->comm,m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);
    426:   /* loop over local fine grid nodes counting interpolating points */
    427:   for (l=l_start; l<l_start+p_f; l++) {
    428:     for (j=j_start; j<j_start+n_f; j++) {
    429:       for (i=i_start; i<i_start+m_f; i++) {
    430:         /* convert to local "natural" numbering and then to PETSc global numbering */
    431:         row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
    432:         i_c = (i/ratioi);
    433:         j_c = (j/ratioj);
    434:         l_c = (l/ratiok);
    435:         if (l_c < l_start_ghost_c) SETERRQ3(PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\
    436:           l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c);
    437:         if (j_c < j_start_ghost_c) SETERRQ3(PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\
    438:           j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
    439:         if (i_c < i_start_ghost_c) SETERRQ3(PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\
    440:           i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);

    442:         /*
    443:          Only include those interpolation points that are truly
    444:          nonzero. Note this is very important for final grid lines
    445:          in x and y directions; since they have no right/top neighbors
    446:         */
    447:         nc       = 0;
    448:         col      = dof*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
    449:         cols[nc++] = idx_c[col]/dof;
    450:         if (i_c*ratioi != i) {
    451:           cols[nc++] = idx_c[col+dof]/dof;
    452:         }
    453:         if (j_c*ratioj != j) {
    454:           cols[nc++] = idx_c[col+m_ghost_c*dof]/dof;
    455:         }
    456:         if (l_c*ratiok != l) {
    457:           cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof;
    458:         }
    459:         if (j_c*ratioj != j && i_c*ratioi != i) {
    460:           cols[nc++] = idx_c[col+(m_ghost_c+1)*dof]/dof;
    461:         }
    462:         if (j_c*ratioj != j && l_c*ratiok != l) {
    463:           cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;
    464:         }
    465:         if (i_c*ratioi != i && l_c*ratiok != l) {
    466:           cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof;
    467:         }
    468:         if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
    469:           cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof;
    470:         }
    471:         MatPreallocateSet(row,nc,cols,dnz,onz);
    472:       }
    473:     }
    474:   }
    475:   MatCreate(dac->comm,&mat);
    476:   MatSetSizes(mat,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz);
    477:   MatSetType(mat,MATAIJ);
    478:   MatSeqAIJSetPreallocation(mat,0,dnz);
    479:   MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);
    480:   MatPreallocateFinalize(dnz,onz);
    481:   if (!DAXPeriodic(pt) && !DAYPeriodic(pt) && !DAZPeriodic(pt)) {MatSetOption(mat,MAT_COLUMNS_SORTED);}

    483:   DAGetCoordinates(daf,&vcoors);
    484:   if (vcoors) {
    485:     DAGetGhostedCoordinates(dac,&cvcoors);
    486:     DAVecGetArray(daf->da_coordinates,vcoors,&coors);
    487:     DAVecGetArray(dac->da_coordinates,cvcoors,&ccoors);
    488:   }

    490:   /* loop over local fine grid nodes setting interpolation for those*/
    491:   for (l=l_start; l<l_start+p_f; l++) {
    492:     for (j=j_start; j<j_start+n_f; j++) {
    493:       for (i=i_start; i<i_start+m_f; i++) {
    494:         /* convert to local "natural" numbering and then to PETSc global numbering */
    495:         row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;

    497:         i_c = (i/ratioi);
    498:         j_c = (j/ratioj);
    499:         l_c = (l/ratiok);

    501:         /*
    502:            Only include those interpolation points that are truly
    503:            nonzero. Note this is very important for final grid lines
    504:            in x and y directions; since they have no right/top neighbors
    505:         */
    506:         if (coors) {
    507:           /* only access the next coors point if we know there is one */
    508:           /* note this is dangerous because x may not exactly equal ZERO */
    509:           x = (coors[l][j][i].x - ccoors[l_c][j_c][i_c].x);
    510:           if (PetscAbsScalar(x) != 0.0) x = x/(ccoors[l_c][j_c][i_c+1].x - ccoors[l_c][j_c][i_c].x);
    511:           y = (coors[l][j][i].y - ccoors[l_c][j_c][i_c].y);
    512:           if (PetscAbsScalar(y) != 0.0) y = y/(ccoors[l_c][j_c+1][i_c].y - ccoors[l_c][j_c][i_c].y);
    513:           z = (coors[l][j][i].z - ccoors[l_c][j_c][i_c].z);
    514:           if (PetscAbsScalar(z) != 0.0) z = z/(ccoors[l_c+1][j_c][i_c].z - ccoors[l_c][j_c][i_c].z);
    515:         } else {
    516:           x  = ((double)(i - i_c*ratioi))/((double)ratioi);
    517:           y  = ((double)(j - j_c*ratioj))/((double)ratioj);
    518:           z  = ((double)(l - l_c*ratiok))/((double)ratiok);
    519:         }
    520:         nc = 0;
    521:         /* one left and below; or we are right on it */
    522:         col      = dof*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));

    524:         cols[nc] = idx_c[col]/dof;
    525:         v[nc++]  = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));

    527:         if (i_c*ratioi != i) {
    528:           cols[nc] = idx_c[col+dof]/dof;
    529:           v[nc++]  = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
    530:         }

    532:         if (j_c*ratioj != j) {
    533:           cols[nc] = idx_c[col+m_ghost_c*dof]/dof;
    534:           v[nc++]  = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
    535:         }

    537:         if (l_c*ratiok != l) {
    538:           cols[nc] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof;
    539:           v[nc++]  = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
    540:         }

    542:         if (j_c*ratioj != j && i_c*ratioi != i) {
    543:           cols[nc] = idx_c[col+(m_ghost_c+1)*dof]/dof;
    544:           v[nc++]  = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
    545:         }

    547:         if (j_c*ratioj != j && l_c*ratiok != l) {
    548:           cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;
    549:           v[nc++]  = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
    550:         }

    552:         if (i_c*ratioi != i && l_c*ratiok != l) {
    553:           cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof;
    554:           v[nc++]  = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
    555:         }

    557:         if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
    558:           cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof;
    559:           v[nc++]  = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
    560:         }
    561:         MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
    562:       }
    563:     }
    564:   }
    565:   if (vcoors) {
    566:     DAVecRestoreArray(daf->da_coordinates,vcoors,&coors);
    567:     DAVecRestoreArray(dac->da_coordinates,cvcoors,&ccoors);
    568:   }
    569:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
    570:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

    572:   MatCreateMAIJ(mat,dof,A);
    573:   MatDestroy(mat);
    574:   PetscLogFlops(13*m_f*n_f);
    575:   return(0);
    576: }

    580: /*@C
    581:    DAGetInterpolation - Gets an interpolation matrix that maps between
    582:    grids associated with two DAs.

    584:    Collective on DA

    586:    Input Parameters:
    587: +  dac - the coarse grid DA
    588: -  daf - the fine grid DA

    590:    Output Parameters:
    591: +  A - the interpolation matrix
    592: -  scale - a scaling vector used to scale the coarse grid restricted vector before applying the
    593:            grid function or grid Jacobian to it.

    595:    Level: intermediate

    597: .keywords: interpolation, restriction, multigrid

    599: .seealso: DARefine(), DAGetInjection()
    600: @*/
    601: PetscErrorCode PETSCDM_DLLEXPORT DAGetInterpolation(DA dac,DA daf,Mat *A,Vec *scale)
    602: {
    604:   PetscInt       dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
    605:   DAPeriodicType wrapc,wrapf;
    606:   DAStencilType  stc,stf;

    614:   DAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&wrapc,&stc);
    615:   DAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&wrapf,&stf);
    616:   if (dimc != dimf) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Dimensions of DA do not match %D %D",dimc,dimf);
    617:   if (dofc != doff) SETERRQ2(PETSC_ERR_ARG_INCOMP,"DOF of DA do not match %D %D",dofc,doff);
    618:   if (sc != sf) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Stencil width of DA do not match %D %D",sc,sf);
    619:   if (wrapc != wrapf) SETERRQ(PETSC_ERR_ARG_INCOMP,"Periodic type different in two DAs");
    620:   if (stc != stf) SETERRQ(PETSC_ERR_ARG_INCOMP,"Stencil type different in two DAs");
    621:   if (Mc < 2 && Mf > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
    622:   if (dimc > 1 && Nc < 2 && Nf > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction");
    623:   if (dimc > 2 && Pc < 2 && Pf > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction");

    625:   if (dac->interptype == DA_Q1){
    626:     if (dimc == 1){
    627:       DAGetInterpolation_1D_Q1(dac,daf,A);
    628:     } else if (dimc == 2){
    629:       DAGetInterpolation_2D_Q1(dac,daf,A);
    630:     } else if (dimc == 3){
    631:       DAGetInterpolation_3D_Q1(dac,daf,A);
    632:     } else {
    633:       SETERRQ2(PETSC_ERR_SUP,"No support for this DA dimension %D for interpolation type %d",dimc,(int)dac->interptype);
    634:     }
    635:   } else if (dac->interptype == DA_Q0){
    636:     if (dimc == 1){
    637:       DAGetInterpolation_1D_Q0(dac,daf,A);
    638:     } else {
    639:       SETERRQ2(PETSC_ERR_SUP,"No support for this DA dimension %D for interpolation type %d",dimc,(int)dac->interptype);
    640:     }
    641:   }
    642:   if (scale) {
    643:     DMGetInterpolationScale((DM)dac,(DM)daf,*A,scale);
    644:   }
    645:   return(0);
    646: }

    650: PetscErrorCode DAGetInjection_2D(DA dac,DA daf,VecScatter *inject)
    651: {
    653:   PetscInt       i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof;
    654:   PetscInt       m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c;
    655:   PetscInt       row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj;
    656:   PetscInt       i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c;
    657:   PetscInt       *cols;
    658:   DAPeriodicType pt;
    659:   Vec            vecf,vecc;
    660:   IS             isf;

    664:   DAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&pt,0);
    665:   DAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0);
    666:   if (DAXPeriodic(pt)){
    667:     ratioi = mx/Mx;
    668:     if (ratioi*Mx != mx) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx  must be integer: mx %D Mx %D",mx,Mx);
    669:   } else {
    670:     ratioi = (mx-1)/(Mx-1);
    671:     if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx);
    672:   }
    673:   if (DAYPeriodic(pt)){
    674:     ratioj = my/My;
    675:     if (ratioj*My != my) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My  must be integer: my %D My %D",my,My);
    676:   } else {
    677:     ratioj = (my-1)/(My-1);
    678:     if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My);
    679:   }

    682:   DAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);
    683:   DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
    684:   DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);

    686:   DAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
    687:   DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
    688:   DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);

    691:   /* loop over local fine grid nodes setting interpolation for those*/
    692:   nc = 0;
    693:   PetscMalloc(n_f*m_f*sizeof(PetscInt),&cols);
    694:   for (j=j_start; j<j_start+n_f; j++) {
    695:     for (i=i_start; i<i_start+m_f; i++) {

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

    700:       if (j_c < j_start_ghost_c) SETERRQ3(PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\
    701:     j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c);
    702:       if (i_c < i_start_ghost_c) SETERRQ3(PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\
    703:     i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c);

    705:       if (i_c*ratioi == i && j_c*ratioj == j) {
    706:         /* convert to local "natural" numbering and then to PETSc global numbering */
    707:         row    = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))];
    708:         cols[nc++] = row;
    709:       }
    710:     }
    711:   }

    713:   ISCreateBlock(daf->comm,dof,nc,cols,&isf);
    714:   DAGetGlobalVector(dac,&vecc);
    715:   DAGetGlobalVector(daf,&vecf);
    716:   VecScatterCreate(vecf,isf,vecc,PETSC_NULL,inject);
    717:   DARestoreGlobalVector(dac,&vecc);
    718:   DARestoreGlobalVector(daf,&vecf);
    719:   ISDestroy(isf);
    720:   PetscFree(cols);
    721:   return(0);
    722: }

    726: /*@C
    727:    DAGetInjection - Gets an injection matrix that maps between
    728:    grids associated with two DAs.

    730:    Collective on DA

    732:    Input Parameters:
    733: +  dac - the coarse grid DA
    734: -  daf - the fine grid DA

    736:    Output Parameters:
    737: .  inject - the injection scatter

    739:    Level: intermediate

    741: .keywords: interpolation, restriction, multigrid, injection

    743: .seealso: DARefine(), DAGetInterpolation()
    744: @*/
    745: PetscErrorCode PETSCDM_DLLEXPORT DAGetInjection(DA dac,DA daf,VecScatter *inject)
    746: {
    748:   PetscInt       dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf;
    749:   DAPeriodicType wrapc,wrapf;
    750:   DAStencilType  stc,stf;

    757:   DAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&wrapc,&stc);
    758:   DAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&wrapf,&stf);
    759:   if (dimc != dimf) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Dimensions of DA do not match %D %D",dimc,dimf);
    760:   if (dofc != doff) SETERRQ2(PETSC_ERR_ARG_INCOMP,"DOF of DA do not match %D %D",dofc,doff);
    761:   if (sc != sf) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Stencil width of DA do not match %D %D",sc,sf);
    762:   if (wrapc != wrapf) SETERRQ(PETSC_ERR_ARG_INCOMP,"Periodic type different in two DAs");
    763:   if (stc != stf) SETERRQ(PETSC_ERR_ARG_INCOMP,"Stencil type different in two DAs");
    764:   if (Mc < 2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction");
    765:   if (dimc > 1 && Nc < 2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction");
    766:   if (dimc > 2 && Pc < 2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction");

    768:   if (dimc == 2){
    769:     DAGetInjection_2D(dac,daf,inject);
    770:   } else {
    771:     SETERRQ1(PETSC_ERR_SUP,"No support for this DA dimension %D",dimc);
    772:   }
    773:   return(0);
    774: }