1 2 /* 3 Code for interpolating between grids represented by DMDAs 4 */ 5 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 13 14 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 15 16 /*@ 17 DMCreateInterpolationScale - Forms L = R*1/diag(R*1) - L.*v is like a coarse grid average of the 18 nearby fine grid points. 19 20 Input Parameters: 21 + dac - DM that defines a coarse mesh 22 . daf - DM that defines a fine mesh 23 - mat - the restriction (or interpolation operator) from fine to coarse 24 25 Output Parameter: 26 . scale - the scaled vector 27 28 Level: developer 29 30 .seealso: DMCreateInterpolation() 31 32 @*/ 33 PetscErrorCode DMCreateInterpolationScale(DM dac,DM daf,Mat mat,Vec *scale) 34 { 35 PetscErrorCode ierr; 36 Vec fine; 37 PetscScalar one = 1.0; 38 39 PetscFunctionBegin; 40 ierr = DMCreateGlobalVector(daf,&fine);CHKERRQ(ierr); 41 ierr = DMCreateGlobalVector(dac,scale);CHKERRQ(ierr); 42 ierr = VecSet(fine,one);CHKERRQ(ierr); 43 ierr = MatRestrict(mat,fine,*scale);CHKERRQ(ierr); 44 ierr = VecDestroy(&fine);CHKERRQ(ierr); 45 ierr = VecReciprocal(*scale);CHKERRQ(ierr); 46 PetscFunctionReturn(0); 47 } 48 49 PetscErrorCode DMCreateInterpolation_DA_1D_Q1(DM dac,DM daf,Mat *A) 50 { 51 PetscErrorCode ierr; 52 PetscInt i,i_start,m_f,Mx; 53 const PetscInt *idx_f,*idx_c; 54 PetscInt m_ghost,m_ghost_c; 55 PetscInt row,col,i_start_ghost,mx,m_c,nc,ratio; 56 PetscInt i_c,i_start_c,i_start_ghost_c,cols[2],dof; 57 PetscScalar v[2],x; 58 Mat mat; 59 DMBoundaryType bx; 60 ISLocalToGlobalMapping ltog_f,ltog_c; 61 62 63 PetscFunctionBegin; 64 ierr = DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);CHKERRQ(ierr); 65 ierr = DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 66 if (bx == DM_BOUNDARY_PERIODIC) { 67 ratio = mx/Mx; 68 if (ratio*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx); 69 } else { 70 ratio = (mx-1)/(Mx-1); 71 if (ratio*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx); 72 } 73 74 ierr = DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);CHKERRQ(ierr); 75 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);CHKERRQ(ierr); 76 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 77 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 78 79 ierr = DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);CHKERRQ(ierr); 80 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);CHKERRQ(ierr); 81 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 82 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 83 84 /* create interpolation matrix */ 85 ierr = MatCreate(PetscObjectComm((PetscObject)dac),&mat);CHKERRQ(ierr); 86 ierr = MatSetSizes(mat,m_f,m_c,mx,Mx);CHKERRQ(ierr); 87 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 88 ierr = MatSeqAIJSetPreallocation(mat,2,NULL);CHKERRQ(ierr); 89 ierr = MatMPIAIJSetPreallocation(mat,2,NULL,1,NULL);CHKERRQ(ierr); 90 91 /* loop over local fine grid nodes setting interpolation for those*/ 92 if (!NEWVERSION) { 93 94 for (i=i_start; i<i_start+m_f; i++) { 95 /* convert to local "natural" numbering and then to PETSc global numbering */ 96 row = idx_f[i-i_start_ghost]; 97 98 i_c = (i/ratio); /* coarse grid node to left of fine grid node */ 99 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 100 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 101 102 /* 103 Only include those interpolation points that are truly 104 nonzero. Note this is very important for final grid lines 105 in x direction; since they have no right neighbor 106 */ 107 x = ((PetscReal)(i - i_c*ratio))/((PetscReal)ratio); 108 nc = 0; 109 /* one left and below; or we are right on it */ 110 col = (i_c-i_start_ghost_c); 111 cols[nc] = idx_c[col]; 112 v[nc++] = -x + 1.0; 113 /* one right? */ 114 if (i_c*ratio != i) { 115 cols[nc] = idx_c[col+1]; 116 v[nc++] = x; 117 } 118 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 119 } 120 121 } else { 122 PetscScalar *xi; 123 PetscInt li,nxi,n; 124 PetscScalar Ni[2]; 125 126 /* compute local coordinate arrays */ 127 nxi = ratio + 1; 128 ierr = PetscMalloc1(nxi,&xi);CHKERRQ(ierr); 129 for (li=0; li<nxi; li++) { 130 xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1)); 131 } 132 133 for (i=i_start; i<i_start+m_f; i++) { 134 /* convert to local "natural" numbering and then to PETSc global numbering */ 135 row = idx_f[(i-i_start_ghost)]; 136 137 i_c = (i/ratio); /* coarse grid node to left of fine grid node */ 138 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 139 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 140 141 /* remainders */ 142 li = i - ratio * (i/ratio); 143 if (i==mx-1) li = nxi-1; 144 145 /* corners */ 146 col = (i_c-i_start_ghost_c); 147 cols[0] = idx_c[col]; 148 Ni[0] = 1.0; 149 if ((li==0) || (li==nxi-1)) { 150 ierr = MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr); 151 continue; 152 } 153 154 /* edges + interior */ 155 /* remainders */ 156 if (i==mx-1) i_c--; 157 158 col = (i_c-i_start_ghost_c); 159 cols[0] = idx_c[col]; /* one left and below; or we are right on it */ 160 cols[1] = idx_c[col+1]; 161 162 Ni[0] = 0.5*(1.0-xi[li]); 163 Ni[1] = 0.5*(1.0+xi[li]); 164 for (n=0; n<2; n++) { 165 if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1; 166 } 167 ierr = MatSetValues(mat,1,&row,2,cols,Ni,INSERT_VALUES);CHKERRQ(ierr); 168 } 169 ierr = PetscFree(xi);CHKERRQ(ierr); 170 } 171 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 172 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 173 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 174 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 175 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 176 ierr = MatDestroy(&mat);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 PetscErrorCode DMCreateInterpolation_DA_1D_Q0(DM dac,DM daf,Mat *A) 181 { 182 PetscErrorCode ierr; 183 PetscInt i,i_start,m_f,Mx; 184 const PetscInt *idx_f,*idx_c; 185 ISLocalToGlobalMapping ltog_f,ltog_c; 186 PetscInt m_ghost,m_ghost_c; 187 PetscInt row,col,i_start_ghost,mx,m_c,nc,ratio; 188 PetscInt i_c,i_start_c,i_start_ghost_c,cols[2],dof; 189 PetscScalar v[2],x; 190 Mat mat; 191 DMBoundaryType bx; 192 193 PetscFunctionBegin; 194 ierr = DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);CHKERRQ(ierr); 195 ierr = DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 196 if (bx == DM_BOUNDARY_PERIODIC) { 197 if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx); 198 ratio = mx/Mx; 199 if (ratio*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx); 200 } else { 201 if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx); 202 ratio = (mx-1)/(Mx-1); 203 if (ratio*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx); 204 } 205 206 ierr = DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);CHKERRQ(ierr); 207 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);CHKERRQ(ierr); 208 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 209 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 210 211 ierr = DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);CHKERRQ(ierr); 212 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);CHKERRQ(ierr); 213 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 214 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 215 216 /* create interpolation matrix */ 217 ierr = MatCreate(PetscObjectComm((PetscObject)dac),&mat);CHKERRQ(ierr); 218 ierr = MatSetSizes(mat,m_f,m_c,mx,Mx);CHKERRQ(ierr); 219 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 220 ierr = MatSeqAIJSetPreallocation(mat,2,NULL);CHKERRQ(ierr); 221 ierr = MatMPIAIJSetPreallocation(mat,2,NULL,0,NULL);CHKERRQ(ierr); 222 223 /* loop over local fine grid nodes setting interpolation for those*/ 224 for (i=i_start; i<i_start+m_f; i++) { 225 /* convert to local "natural" numbering and then to PETSc global numbering */ 226 row = idx_f[(i-i_start_ghost)]; 227 228 i_c = (i/ratio); /* coarse grid node to left of fine grid node */ 229 230 /* 231 Only include those interpolation points that are truly 232 nonzero. Note this is very important for final grid lines 233 in x direction; since they have no right neighbor 234 */ 235 x = ((PetscReal)(i - i_c*ratio))/((PetscReal)ratio); 236 nc = 0; 237 /* one left and below; or we are right on it */ 238 col = (i_c-i_start_ghost_c); 239 cols[nc] = idx_c[col]; 240 v[nc++] = -x + 1.0; 241 /* one right? */ 242 if (i_c*ratio != i) { 243 cols[nc] = idx_c[col+1]; 244 v[nc++] = x; 245 } 246 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 247 } 248 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 249 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 250 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 251 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 252 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 253 ierr = MatDestroy(&mat);CHKERRQ(ierr); 254 ierr = PetscLogFlops(5.0*m_f);CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257 258 PetscErrorCode DMCreateInterpolation_DA_2D_Q1(DM dac,DM daf,Mat *A) 259 { 260 PetscErrorCode ierr; 261 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof; 262 const PetscInt *idx_c,*idx_f; 263 ISLocalToGlobalMapping ltog_f,ltog_c; 264 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c,*dnz,*onz; 265 PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj; 266 PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c,col_shift,col_scale; 267 PetscMPIInt size_c,size_f,rank_f; 268 PetscScalar v[4],x,y; 269 Mat mat; 270 DMBoundaryType bx,by; 271 272 PetscFunctionBegin; 273 ierr = DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);CHKERRQ(ierr); 274 ierr = DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 275 if (bx == DM_BOUNDARY_PERIODIC) { 276 if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx); 277 ratioi = mx/Mx; 278 if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx); 279 } else { 280 if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx); 281 ratioi = (mx-1)/(Mx-1); 282 if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx); 283 } 284 if (by == DM_BOUNDARY_PERIODIC) { 285 if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My); 286 ratioj = my/My; 287 if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My); 288 } else { 289 if (My < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be greater than 1",My); 290 ratioj = (my-1)/(My-1); 291 if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My); 292 } 293 294 ierr = DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 295 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 296 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 297 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 298 299 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 300 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 301 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 302 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 303 304 /* 305 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 306 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 307 processors). It's effective length is hence 4 times its normal length, this is 308 why the col_scale is multiplied by the interpolation matrix column sizes. 309 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 310 copy of the coarse vector. A bit of a hack but you do better. 311 312 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 313 */ 314 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);CHKERRQ(ierr); 315 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRQ(ierr); 316 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRQ(ierr); 317 col_scale = size_f/size_c; 318 col_shift = Mx*My*(rank_f/size_c); 319 320 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f,col_scale*m_c*n_c,dnz,onz);CHKERRQ(ierr); 321 for (j=j_start; j<j_start+n_f; j++) { 322 for (i=i_start; i<i_start+m_f; i++) { 323 /* convert to local "natural" numbering and then to PETSc global numbering */ 324 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 325 326 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 327 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 328 329 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 330 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 331 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 332 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 333 334 /* 335 Only include those interpolation points that are truly 336 nonzero. Note this is very important for final grid lines 337 in x and y directions; since they have no right/top neighbors 338 */ 339 nc = 0; 340 /* one left and below; or we are right on it */ 341 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 342 cols[nc++] = col_shift + idx_c[col]; 343 /* one right and below */ 344 if (i_c*ratioi != i) cols[nc++] = col_shift + idx_c[col+1]; 345 /* one left and above */ 346 if (j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+m_ghost_c]; 347 /* one right and above */ 348 if (i_c*ratioi != i && j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+(m_ghost_c+1)]; 349 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 350 } 351 } 352 ierr = MatCreate(PetscObjectComm((PetscObject)daf),&mat);CHKERRQ(ierr); 353 ierr = MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);CHKERRQ(ierr); 354 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 355 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 356 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 357 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 358 359 /* loop over local fine grid nodes setting interpolation for those*/ 360 if (!NEWVERSION) { 361 362 for (j=j_start; j<j_start+n_f; j++) { 363 for (i=i_start; i<i_start+m_f; i++) { 364 /* convert to local "natural" numbering and then to PETSc global numbering */ 365 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 366 367 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 368 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 369 370 /* 371 Only include those interpolation points that are truly 372 nonzero. Note this is very important for final grid lines 373 in x and y directions; since they have no right/top neighbors 374 */ 375 x = ((PetscReal)(i - i_c*ratioi))/((PetscReal)ratioi); 376 y = ((PetscReal)(j - j_c*ratioj))/((PetscReal)ratioj); 377 378 nc = 0; 379 /* one left and below; or we are right on it */ 380 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 381 cols[nc] = col_shift + idx_c[col]; 382 v[nc++] = x*y - x - y + 1.0; 383 /* one right and below */ 384 if (i_c*ratioi != i) { 385 cols[nc] = col_shift + idx_c[col+1]; 386 v[nc++] = -x*y + x; 387 } 388 /* one left and above */ 389 if (j_c*ratioj != j) { 390 cols[nc] = col_shift + idx_c[col+m_ghost_c]; 391 v[nc++] = -x*y + y; 392 } 393 /* one right and above */ 394 if (j_c*ratioj != j && i_c*ratioi != i) { 395 cols[nc] = col_shift + idx_c[col+(m_ghost_c+1)]; 396 v[nc++] = x*y; 397 } 398 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 399 } 400 } 401 402 } else { 403 PetscScalar Ni[4]; 404 PetscScalar *xi,*eta; 405 PetscInt li,nxi,lj,neta; 406 407 /* compute local coordinate arrays */ 408 nxi = ratioi + 1; 409 neta = ratioj + 1; 410 ierr = PetscMalloc1(nxi,&xi);CHKERRQ(ierr); 411 ierr = PetscMalloc1(neta,&eta);CHKERRQ(ierr); 412 for (li=0; li<nxi; li++) { 413 xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1)); 414 } 415 for (lj=0; lj<neta; lj++) { 416 eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1)); 417 } 418 419 /* loop over local fine grid nodes setting interpolation for those*/ 420 for (j=j_start; j<j_start+n_f; j++) { 421 for (i=i_start; i<i_start+m_f; i++) { 422 /* convert to local "natural" numbering and then to PETSc global numbering */ 423 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 424 425 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 426 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 427 428 /* remainders */ 429 li = i - ratioi * (i/ratioi); 430 if (i==mx-1) li = nxi-1; 431 lj = j - ratioj * (j/ratioj); 432 if (j==my-1) lj = neta-1; 433 434 /* corners */ 435 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 436 cols[0] = col_shift + idx_c[col]; /* left, below */ 437 Ni[0] = 1.0; 438 if ((li==0) || (li==nxi-1)) { 439 if ((lj==0) || (lj==neta-1)) { 440 ierr = MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr); 441 continue; 442 } 443 } 444 445 /* edges + interior */ 446 /* remainders */ 447 if (i==mx-1) i_c--; 448 if (j==my-1) j_c--; 449 450 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 451 cols[0] = col_shift + idx_c[col]; /* left, below */ 452 cols[1] = col_shift + idx_c[col+1]; /* right, below */ 453 cols[2] = col_shift + idx_c[col+m_ghost_c]; /* left, above */ 454 cols[3] = col_shift + idx_c[col+(m_ghost_c+1)]; /* right, above */ 455 456 Ni[0] = 0.25*(1.0-xi[li])*(1.0-eta[lj]); 457 Ni[1] = 0.25*(1.0+xi[li])*(1.0-eta[lj]); 458 Ni[2] = 0.25*(1.0-xi[li])*(1.0+eta[lj]); 459 Ni[3] = 0.25*(1.0+xi[li])*(1.0+eta[lj]); 460 461 nc = 0; 462 if (PetscAbsScalar(Ni[0])<1.0e-32) cols[0]=-1; 463 if (PetscAbsScalar(Ni[1])<1.0e-32) cols[1]=-1; 464 if (PetscAbsScalar(Ni[2])<1.0e-32) cols[2]=-1; 465 if (PetscAbsScalar(Ni[3])<1.0e-32) cols[3]=-1; 466 467 ierr = MatSetValues(mat,1,&row,4,cols,Ni,INSERT_VALUES);CHKERRQ(ierr); 468 } 469 } 470 ierr = PetscFree(xi);CHKERRQ(ierr); 471 ierr = PetscFree(eta);CHKERRQ(ierr); 472 } 473 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 474 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 475 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 476 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 477 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 478 ierr = MatDestroy(&mat);CHKERRQ(ierr); 479 PetscFunctionReturn(0); 480 } 481 482 /* 483 Contributed by Andrei Draganescu <aidraga@sandia.gov> 484 */ 485 PetscErrorCode DMCreateInterpolation_DA_2D_Q0(DM dac,DM daf,Mat *A) 486 { 487 PetscErrorCode ierr; 488 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof; 489 const PetscInt *idx_c,*idx_f; 490 ISLocalToGlobalMapping ltog_f,ltog_c; 491 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c,*dnz,*onz; 492 PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj; 493 PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c,col_shift,col_scale; 494 PetscMPIInt size_c,size_f,rank_f; 495 PetscScalar v[4]; 496 Mat mat; 497 DMBoundaryType bx,by; 498 499 PetscFunctionBegin; 500 ierr = DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);CHKERRQ(ierr); 501 ierr = DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 502 if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx); 503 if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My); 504 ratioi = mx/Mx; 505 ratioj = my/My; 506 if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x"); 507 if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y"); 508 if (ratioi != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 2"); 509 if (ratioj != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 2"); 510 511 ierr = DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 512 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 513 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 514 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 515 516 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 517 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 518 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 519 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 520 521 /* 522 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 523 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 524 processors). It's effective length is hence 4 times its normal length, this is 525 why the col_scale is multiplied by the interpolation matrix column sizes. 526 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 527 copy of the coarse vector. A bit of a hack but you do better. 528 529 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 530 */ 531 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);CHKERRQ(ierr); 532 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRQ(ierr); 533 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRQ(ierr); 534 col_scale = size_f/size_c; 535 col_shift = Mx*My*(rank_f/size_c); 536 537 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f,col_scale*m_c*n_c,dnz,onz);CHKERRQ(ierr); 538 for (j=j_start; j<j_start+n_f; j++) { 539 for (i=i_start; i<i_start+m_f; i++) { 540 /* convert to local "natural" numbering and then to PETSc global numbering */ 541 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 542 543 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 544 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 545 546 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 547 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 548 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 549 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 550 551 /* 552 Only include those interpolation points that are truly 553 nonzero. Note this is very important for final grid lines 554 in x and y directions; since they have no right/top neighbors 555 */ 556 nc = 0; 557 /* one left and below; or we are right on it */ 558 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 559 cols[nc++] = col_shift + idx_c[col]; 560 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 561 } 562 } 563 ierr = MatCreate(PetscObjectComm((PetscObject)daf),&mat);CHKERRQ(ierr); 564 ierr = MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);CHKERRQ(ierr); 565 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 566 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 567 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 568 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 569 570 /* loop over local fine grid nodes setting interpolation for those*/ 571 for (j=j_start; j<j_start+n_f; j++) { 572 for (i=i_start; i<i_start+m_f; i++) { 573 /* convert to local "natural" numbering and then to PETSc global numbering */ 574 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 575 576 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 577 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 578 nc = 0; 579 /* one left and below; or we are right on it */ 580 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 581 cols[nc] = col_shift + idx_c[col]; 582 v[nc++] = 1.0; 583 584 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 585 } 586 } 587 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 588 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 589 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 590 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 591 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 592 ierr = MatDestroy(&mat);CHKERRQ(ierr); 593 ierr = PetscLogFlops(13.0*m_f*n_f);CHKERRQ(ierr); 594 PetscFunctionReturn(0); 595 } 596 597 /* 598 Contributed by Jianming Yang <jianming-yang@uiowa.edu> 599 */ 600 PetscErrorCode DMCreateInterpolation_DA_3D_Q0(DM dac,DM daf,Mat *A) 601 { 602 PetscErrorCode ierr; 603 PetscInt i,j,l,i_start,j_start,l_start,m_f,n_f,p_f,Mx,My,Mz,dof; 604 const PetscInt *idx_c,*idx_f; 605 ISLocalToGlobalMapping ltog_f,ltog_c; 606 PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c,nc,*dnz,*onz; 607 PetscInt row,col,i_start_ghost,j_start_ghost,l_start_ghost,cols[8],mx,m_c,my,n_c,mz,p_c,ratioi,ratioj,ratiol; 608 PetscInt i_c,j_c,l_c,i_start_c,j_start_c,l_start_c,i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,col_shift,col_scale; 609 PetscMPIInt size_c,size_f,rank_f; 610 PetscScalar v[8]; 611 Mat mat; 612 DMBoundaryType bx,by,bz; 613 614 PetscFunctionBegin; 615 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 616 if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx); 617 if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My); 618 if (!Mz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be positive",Mz); 619 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 620 ratioi = mx/Mx; 621 ratioj = my/My; 622 ratiol = mz/Mz; 623 if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x"); 624 if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y"); 625 if (ratiol*Mz != mz) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in z"); 626 if (ratioi != 2 && ratioi != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 1 or 2"); 627 if (ratioj != 2 && ratioj != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 1 or 2"); 628 if (ratiol != 2 && ratiol != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in z must be 1 or 2"); 629 630 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 631 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 632 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 633 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 634 635 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 636 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr); 637 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 638 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 639 640 /* 641 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 642 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 643 processors). It's effective length is hence 4 times its normal length, this is 644 why the col_scale is multiplied by the interpolation matrix column sizes. 645 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 646 copy of the coarse vector. A bit of a hack but you do better. 647 648 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 649 */ 650 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);CHKERRQ(ierr); 651 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRQ(ierr); 652 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRQ(ierr); 653 col_scale = size_f/size_c; 654 col_shift = Mx*My*Mz*(rank_f/size_c); 655 656 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f*p_f,col_scale*m_c*n_c*p_c,dnz,onz);CHKERRQ(ierr); 657 for (l=l_start; l<l_start+p_f; l++) { 658 for (j=j_start; j<j_start+n_f; j++) { 659 for (i=i_start; i<i_start+m_f; i++) { 660 /* convert to local "natural" numbering and then to PETSc global numbering */ 661 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 662 663 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 664 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 665 l_c = (l/ratiol); 666 667 if (l_c < l_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 668 l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c); 669 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 670 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 671 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 672 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 673 674 /* 675 Only include those interpolation points that are truly 676 nonzero. Note this is very important for final grid lines 677 in x and y directions; since they have no right/top neighbors 678 */ 679 nc = 0; 680 /* one left and below; or we are right on it */ 681 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 682 cols[nc++] = col_shift + idx_c[col]; 683 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 684 } 685 } 686 } 687 ierr = MatCreate(PetscObjectComm((PetscObject)daf),&mat);CHKERRQ(ierr); 688 ierr = MatSetSizes(mat,m_f*n_f*p_f,col_scale*m_c*n_c*p_c,mx*my*mz,col_scale*Mx*My*Mz);CHKERRQ(ierr); 689 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 690 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 691 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 692 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 693 694 /* loop over local fine grid nodes setting interpolation for those*/ 695 for (l=l_start; l<l_start+p_f; l++) { 696 for (j=j_start; j<j_start+n_f; j++) { 697 for (i=i_start; i<i_start+m_f; i++) { 698 /* convert to local "natural" numbering and then to PETSc global numbering */ 699 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 700 701 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 702 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 703 l_c = (l/ratiol); 704 nc = 0; 705 /* one left and below; or we are right on it */ 706 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 707 cols[nc] = col_shift + idx_c[col]; 708 v[nc++] = 1.0; 709 710 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 711 } 712 } 713 } 714 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 715 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 716 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 717 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 718 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 719 ierr = MatDestroy(&mat);CHKERRQ(ierr); 720 ierr = PetscLogFlops(13.0*m_f*n_f*p_f);CHKERRQ(ierr); 721 PetscFunctionReturn(0); 722 } 723 724 PetscErrorCode DMCreateInterpolation_DA_3D_Q1(DM dac,DM daf,Mat *A) 725 { 726 PetscErrorCode ierr; 727 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof,l; 728 const PetscInt *idx_c,*idx_f; 729 ISLocalToGlobalMapping ltog_f,ltog_c; 730 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c,Mz,mz; 731 PetscInt row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok; 732 PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c; 733 PetscInt l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c; 734 PetscInt l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz; 735 PetscScalar v[8],x,y,z; 736 Mat mat; 737 DMBoundaryType bx,by,bz; 738 739 PetscFunctionBegin; 740 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 741 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 742 if (mx == Mx) { 743 ratioi = 1; 744 } else if (bx == DM_BOUNDARY_PERIODIC) { 745 if (!Mx) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be positive",Mx); 746 ratioi = mx/Mx; 747 if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx); 748 } else { 749 if (Mx < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of x coarse grid points %D must be greater than 1",Mx); 750 ratioi = (mx-1)/(Mx-1); 751 if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx); 752 } 753 if (my == My) { 754 ratioj = 1; 755 } else if (by == DM_BOUNDARY_PERIODIC) { 756 if (!My) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be positive",My); 757 ratioj = my/My; 758 if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My); 759 } else { 760 if (My < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of y coarse grid points %D must be greater than 1",My); 761 ratioj = (my-1)/(My-1); 762 if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My); 763 } 764 if (mz == Mz) { 765 ratiok = 1; 766 } else if (bz == DM_BOUNDARY_PERIODIC) { 767 if (!Mz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be positive",Mz); 768 ratiok = mz/Mz; 769 if (ratiok*Mz != mz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mz/Mz must be integer: mz %D Mz %D",mz,Mz); 770 } else { 771 if (Mz < 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of z coarse grid points %D must be greater than 1",Mz); 772 ratiok = (mz-1)/(Mz-1); 773 if (ratiok*(Mz-1) != mz-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %D Mz %D",mz,Mz); 774 } 775 776 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 777 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 778 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 779 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 780 781 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 782 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr); 783 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 784 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 785 786 /* create interpolation matrix, determining exact preallocation */ 787 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)dac),m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);CHKERRQ(ierr); 788 /* loop over local fine grid nodes counting interpolating points */ 789 for (l=l_start; l<l_start+p_f; l++) { 790 for (j=j_start; j<j_start+n_f; j++) { 791 for (i=i_start; i<i_start+m_f; i++) { 792 /* convert to local "natural" numbering and then to PETSc global numbering */ 793 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 794 i_c = (i/ratioi); 795 j_c = (j/ratioj); 796 l_c = (l/ratiok); 797 if (l_c < l_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 798 l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c); 799 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 800 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 801 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 802 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 803 804 /* 805 Only include those interpolation points that are truly 806 nonzero. Note this is very important for final grid lines 807 in x and y directions; since they have no right/top neighbors 808 */ 809 nc = 0; 810 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 811 cols[nc++] = idx_c[col]; 812 if (i_c*ratioi != i) { 813 cols[nc++] = idx_c[col+1]; 814 } 815 if (j_c*ratioj != j) { 816 cols[nc++] = idx_c[col+m_ghost_c]; 817 } 818 if (l_c*ratiok != l) { 819 cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c]; 820 } 821 if (j_c*ratioj != j && i_c*ratioi != i) { 822 cols[nc++] = idx_c[col+(m_ghost_c+1)]; 823 } 824 if (j_c*ratioj != j && l_c*ratiok != l) { 825 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)]; 826 } 827 if (i_c*ratioi != i && l_c*ratiok != l) { 828 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)]; 829 } 830 if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) { 831 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)]; 832 } 833 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 834 } 835 } 836 } 837 ierr = MatCreate(PetscObjectComm((PetscObject)dac),&mat);CHKERRQ(ierr); 838 ierr = MatSetSizes(mat,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz);CHKERRQ(ierr); 839 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 840 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 841 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 842 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 843 844 /* loop over local fine grid nodes setting interpolation for those*/ 845 if (!NEWVERSION) { 846 847 for (l=l_start; l<l_start+p_f; l++) { 848 for (j=j_start; j<j_start+n_f; j++) { 849 for (i=i_start; i<i_start+m_f; i++) { 850 /* convert to local "natural" numbering and then to PETSc global numbering */ 851 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 852 853 i_c = (i/ratioi); 854 j_c = (j/ratioj); 855 l_c = (l/ratiok); 856 857 /* 858 Only include those interpolation points that are truly 859 nonzero. Note this is very important for final grid lines 860 in x and y directions; since they have no right/top neighbors 861 */ 862 x = ((PetscReal)(i - i_c*ratioi))/((PetscReal)ratioi); 863 y = ((PetscReal)(j - j_c*ratioj))/((PetscReal)ratioj); 864 z = ((PetscReal)(l - l_c*ratiok))/((PetscReal)ratiok); 865 866 nc = 0; 867 /* one left and below; or we are right on it */ 868 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c)); 869 870 cols[nc] = idx_c[col]; 871 v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.)); 872 873 if (i_c*ratioi != i) { 874 cols[nc] = idx_c[col+1]; 875 v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.)); 876 } 877 878 if (j_c*ratioj != j) { 879 cols[nc] = idx_c[col+m_ghost_c]; 880 v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.)); 881 } 882 883 if (l_c*ratiok != l) { 884 cols[nc] = idx_c[col+m_ghost_c*n_ghost_c]; 885 v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.)); 886 } 887 888 if (j_c*ratioj != j && i_c*ratioi != i) { 889 cols[nc] = idx_c[col+(m_ghost_c+1)]; 890 v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.)); 891 } 892 893 if (j_c*ratioj != j && l_c*ratiok != l) { 894 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)]; 895 v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.)); 896 } 897 898 if (i_c*ratioi != i && l_c*ratiok != l) { 899 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)]; 900 v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.)); 901 } 902 903 if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) { 904 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)]; 905 v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.)); 906 } 907 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 908 } 909 } 910 } 911 912 } else { 913 PetscScalar *xi,*eta,*zeta; 914 PetscInt li,nxi,lj,neta,lk,nzeta,n; 915 PetscScalar Ni[8]; 916 917 /* compute local coordinate arrays */ 918 nxi = ratioi + 1; 919 neta = ratioj + 1; 920 nzeta = ratiok + 1; 921 ierr = PetscMalloc1(nxi,&xi);CHKERRQ(ierr); 922 ierr = PetscMalloc1(neta,&eta);CHKERRQ(ierr); 923 ierr = PetscMalloc1(nzeta,&zeta);CHKERRQ(ierr); 924 for (li=0; li<nxi; li++) xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1)); 925 for (lj=0; lj<neta; lj++) eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1)); 926 for (lk=0; lk<nzeta; lk++) zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1)); 927 928 for (l=l_start; l<l_start+p_f; l++) { 929 for (j=j_start; j<j_start+n_f; j++) { 930 for (i=i_start; i<i_start+m_f; i++) { 931 /* convert to local "natural" numbering and then to PETSc global numbering */ 932 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 933 934 i_c = (i/ratioi); 935 j_c = (j/ratioj); 936 l_c = (l/ratiok); 937 938 /* remainders */ 939 li = i - ratioi * (i/ratioi); 940 if (i==mx-1) li = nxi-1; 941 lj = j - ratioj * (j/ratioj); 942 if (j==my-1) lj = neta-1; 943 lk = l - ratiok * (l/ratiok); 944 if (l==mz-1) lk = nzeta-1; 945 946 /* corners */ 947 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c)); 948 cols[0] = idx_c[col]; 949 Ni[0] = 1.0; 950 if ((li==0) || (li==nxi-1)) { 951 if ((lj==0) || (lj==neta-1)) { 952 if ((lk==0) || (lk==nzeta-1)) { 953 ierr = MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr); 954 continue; 955 } 956 } 957 } 958 959 /* edges + interior */ 960 /* remainders */ 961 if (i==mx-1) i_c--; 962 if (j==my-1) j_c--; 963 if (l==mz-1) l_c--; 964 965 col = (m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 966 cols[0] = idx_c[col]; /* one left and below; or we are right on it */ 967 cols[1] = idx_c[col+1]; /* one right and below */ 968 cols[2] = idx_c[col+m_ghost_c]; /* one left and above */ 969 cols[3] = idx_c[col+(m_ghost_c+1)]; /* one right and above */ 970 971 cols[4] = idx_c[col+m_ghost_c*n_ghost_c]; /* one left and below and front; or we are right on it */ 972 cols[5] = idx_c[col+(m_ghost_c*n_ghost_c+1)]; /* one right and below, and front */ 973 cols[6] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)]; /* one left and above and front*/ 974 cols[7] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)]; /* one right and above and front */ 975 976 Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]); 977 Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]); 978 Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]); 979 Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]); 980 981 Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]); 982 Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]); 983 Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]); 984 Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]); 985 986 for (n=0; n<8; n++) { 987 if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1; 988 } 989 ierr = MatSetValues(mat,1,&row,8,cols,Ni,INSERT_VALUES);CHKERRQ(ierr); 990 991 } 992 } 993 } 994 ierr = PetscFree(xi);CHKERRQ(ierr); 995 ierr = PetscFree(eta);CHKERRQ(ierr); 996 ierr = PetscFree(zeta);CHKERRQ(ierr); 997 } 998 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 999 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1000 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1001 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1002 1003 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 1004 ierr = MatDestroy(&mat);CHKERRQ(ierr); 1005 PetscFunctionReturn(0); 1006 } 1007 1008 PetscErrorCode DMCreateInterpolation_DA(DM dac,DM daf,Mat *A,Vec *scale) 1009 { 1010 PetscErrorCode ierr; 1011 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1012 DMBoundaryType bxc,byc,bzc,bxf,byf,bzf; 1013 DMDAStencilType stc,stf; 1014 DM_DA *ddc = (DM_DA*)dac->data; 1015 1016 PetscFunctionBegin; 1017 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1018 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1019 PetscValidPointer(A,3); 1020 if (scale) PetscValidPointer(scale,4); 1021 1022 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr); 1023 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr); 1024 if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf); 1025 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff); 1026 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf); 1027 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs"); 1028 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs"); 1029 if (Mc < 2 && Mf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); 1030 if (dimc > 1 && Nc < 2 && Nf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction"); 1031 if (dimc > 2 && Pc < 2 && Pf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction"); 1032 1033 if (ddc->interptype == DMDA_Q1) { 1034 if (dimc == 1) { 1035 ierr = DMCreateInterpolation_DA_1D_Q1(dac,daf,A);CHKERRQ(ierr); 1036 } else if (dimc == 2) { 1037 ierr = DMCreateInterpolation_DA_2D_Q1(dac,daf,A);CHKERRQ(ierr); 1038 } else if (dimc == 3) { 1039 ierr = DMCreateInterpolation_DA_3D_Q1(dac,daf,A);CHKERRQ(ierr); 1040 } else SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype); 1041 } else if (ddc->interptype == DMDA_Q0) { 1042 if (dimc == 1) { 1043 ierr = DMCreateInterpolation_DA_1D_Q0(dac,daf,A);CHKERRQ(ierr); 1044 } else if (dimc == 2) { 1045 ierr = DMCreateInterpolation_DA_2D_Q0(dac,daf,A);CHKERRQ(ierr); 1046 } else if (dimc == 3) { 1047 ierr = DMCreateInterpolation_DA_3D_Q0(dac,daf,A);CHKERRQ(ierr); 1048 } else SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype); 1049 } 1050 if (scale) { 1051 ierr = DMCreateInterpolationScale((DM)dac,(DM)daf,*A,scale);CHKERRQ(ierr); 1052 } 1053 PetscFunctionReturn(0); 1054 } 1055 1056 PetscErrorCode DMCreateInjection_DA_1D(DM dac,DM daf,VecScatter *inject) 1057 { 1058 PetscErrorCode ierr; 1059 PetscInt i,i_start,m_f,Mx,dof; 1060 const PetscInt *idx_f; 1061 ISLocalToGlobalMapping ltog_f; 1062 PetscInt m_ghost,m_ghost_c; 1063 PetscInt row,i_start_ghost,mx,m_c,nc,ratioi; 1064 PetscInt i_start_c,i_start_ghost_c; 1065 PetscInt *cols; 1066 DMBoundaryType bx; 1067 Vec vecf,vecc; 1068 IS isf; 1069 1070 PetscFunctionBegin; 1071 ierr = DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);CHKERRQ(ierr); 1072 ierr = DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 1073 if (bx == DM_BOUNDARY_PERIODIC) { 1074 ratioi = mx/Mx; 1075 if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx); 1076 } else { 1077 ratioi = (mx-1)/(Mx-1); 1078 if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx); 1079 } 1080 1081 ierr = DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);CHKERRQ(ierr); 1082 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);CHKERRQ(ierr); 1083 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 1084 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1085 1086 ierr = DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);CHKERRQ(ierr); 1087 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);CHKERRQ(ierr); 1088 1089 1090 /* loop over local fine grid nodes setting interpolation for those*/ 1091 nc = 0; 1092 ierr = PetscMalloc1(m_f,&cols);CHKERRQ(ierr); 1093 1094 1095 for (i=i_start_c; i<i_start_c+m_c; i++) { 1096 PetscInt i_f = i*ratioi; 1097 1098 if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\ni_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1099 1100 row = idx_f[(i_f-i_start_ghost)]; 1101 cols[nc++] = row; 1102 } 1103 1104 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1105 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1106 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1107 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1108 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr); 1109 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1110 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1111 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1112 PetscFunctionReturn(0); 1113 } 1114 1115 PetscErrorCode DMCreateInjection_DA_2D(DM dac,DM daf,VecScatter *inject) 1116 { 1117 PetscErrorCode ierr; 1118 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof; 1119 const PetscInt *idx_c,*idx_f; 1120 ISLocalToGlobalMapping ltog_f,ltog_c; 1121 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c; 1122 PetscInt row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj; 1123 PetscInt i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c; 1124 PetscInt *cols; 1125 DMBoundaryType bx,by; 1126 Vec vecf,vecc; 1127 IS isf; 1128 1129 PetscFunctionBegin; 1130 ierr = DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);CHKERRQ(ierr); 1131 ierr = DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 1132 if (bx == DM_BOUNDARY_PERIODIC) { 1133 ratioi = mx/Mx; 1134 if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx); 1135 } else { 1136 ratioi = (mx-1)/(Mx-1); 1137 if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx); 1138 } 1139 if (by == DM_BOUNDARY_PERIODIC) { 1140 ratioj = my/My; 1141 if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My); 1142 } else { 1143 ratioj = (my-1)/(My-1); 1144 if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My); 1145 } 1146 1147 ierr = DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 1148 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 1149 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 1150 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1151 1152 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 1153 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 1154 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 1155 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1156 1157 /* loop over local fine grid nodes setting interpolation for those*/ 1158 nc = 0; 1159 ierr = PetscMalloc1(n_f*m_f,&cols);CHKERRQ(ierr); 1160 for (j=j_start_c; j<j_start_c+n_c; j++) { 1161 for (i=i_start_c; i<i_start_c+m_c; i++) { 1162 PetscInt i_f = i*ratioi,j_f = j*ratioj; 1163 if (j_f < j_start_ghost || j_f >= j_start_ghost+n_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 1164 j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost); 1165 if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA\n\ 1166 i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1167 row = idx_f[(m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))]; 1168 cols[nc++] = row; 1169 } 1170 } 1171 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1172 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1173 1174 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1175 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1176 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1177 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr); 1178 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1179 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1180 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1181 PetscFunctionReturn(0); 1182 } 1183 1184 PetscErrorCode DMCreateInjection_DA_3D(DM dac,DM daf,VecScatter *inject) 1185 { 1186 PetscErrorCode ierr; 1187 PetscInt i,j,k,i_start,j_start,k_start,m_f,n_f,p_f,Mx,My,Mz; 1188 PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c; 1189 PetscInt i_start_ghost,j_start_ghost,k_start_ghost; 1190 PetscInt mx,my,mz,ratioi,ratioj,ratiok; 1191 PetscInt i_start_c,j_start_c,k_start_c; 1192 PetscInt m_c,n_c,p_c; 1193 PetscInt i_start_ghost_c,j_start_ghost_c,k_start_ghost_c; 1194 PetscInt row,nc,dof; 1195 const PetscInt *idx_c,*idx_f; 1196 ISLocalToGlobalMapping ltog_f,ltog_c; 1197 PetscInt *cols; 1198 DMBoundaryType bx,by,bz; 1199 Vec vecf,vecc; 1200 IS isf; 1201 1202 PetscFunctionBegin; 1203 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 1204 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 1205 1206 if (bx == DM_BOUNDARY_PERIODIC) { 1207 ratioi = mx/Mx; 1208 if (ratioi*Mx != mx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mx/Mx must be integer: mx %D Mx %D",mx,Mx); 1209 } else { 1210 ratioi = (mx-1)/(Mx-1); 1211 if (ratioi*(Mx-1) != mx-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %D Mx %D",mx,Mx); 1212 } 1213 if (by == DM_BOUNDARY_PERIODIC) { 1214 ratioj = my/My; 1215 if (ratioj*My != my) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: my/My must be integer: my %D My %D",my,My); 1216 } else { 1217 ratioj = (my-1)/(My-1); 1218 if (ratioj*(My-1) != my-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (my - 1)/(My - 1) must be integer: my %D My %D",my,My); 1219 } 1220 if (bz == DM_BOUNDARY_PERIODIC) { 1221 ratiok = mz/Mz; 1222 if (ratiok*Mz != mz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: mz/Mz must be integer: mz %D My %D",mz,Mz); 1223 } else { 1224 ratiok = (mz-1)/(Mz-1); 1225 if (ratiok*(Mz-1) != mz-1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %D Mz %D",mz,Mz); 1226 } 1227 1228 ierr = DMDAGetCorners(daf,&i_start,&j_start,&k_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 1229 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&k_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 1230 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 1231 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1232 1233 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&k_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 1234 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&k_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr); 1235 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 1236 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1237 1238 1239 /* loop over local fine grid nodes setting interpolation for those*/ 1240 nc = 0; 1241 ierr = PetscMalloc1(n_f*m_f*p_f,&cols);CHKERRQ(ierr); 1242 for (k=k_start_c; k<k_start_c+p_c; k++) { 1243 for (j=j_start_c; j<j_start_c+n_c; j++) { 1244 for (i=i_start_c; i<i_start_c+m_c; i++) { 1245 PetscInt i_f = i*ratioi,j_f = j*ratioj,k_f = k*ratiok; 1246 if (k_f < k_start_ghost || k_f >= k_start_ghost+p_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA " 1247 "k_c %D k_f %D fine ghost range [%D,%D]",k,k_f,k_start_ghost,k_start_ghost+p_ghost); 1248 if (j_f < j_start_ghost || j_f >= j_start_ghost+n_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA " 1249 "j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost); 1250 if (i_f < i_start_ghost || i_f >= i_start_ghost+m_ghost) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DMDA must lie over fine DMDA " 1251 "i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1252 row = idx_f[(m_ghost*n_ghost*(k_f-k_start_ghost) + m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))]; 1253 cols[nc++] = row; 1254 } 1255 } 1256 } 1257 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1258 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1259 1260 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1261 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1262 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1263 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr); 1264 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1265 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1266 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1267 PetscFunctionReturn(0); 1268 } 1269 1270 PetscErrorCode DMCreateInjection_DA(DM dac,DM daf,Mat *mat) 1271 { 1272 PetscErrorCode ierr; 1273 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1274 DMBoundaryType bxc,byc,bzc,bxf,byf,bzf; 1275 DMDAStencilType stc,stf; 1276 VecScatter inject = NULL; 1277 1278 PetscFunctionBegin; 1279 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1280 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1281 PetscValidPointer(mat,3); 1282 1283 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr); 1284 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr); 1285 if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf); 1286 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff); 1287 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf); 1288 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs"); 1289 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs"); 1290 if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); 1291 if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction"); 1292 if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction"); 1293 1294 if (dimc == 1) { 1295 ierr = DMCreateInjection_DA_1D(dac,daf,&inject);CHKERRQ(ierr); 1296 } else if (dimc == 2) { 1297 ierr = DMCreateInjection_DA_2D(dac,daf,&inject);CHKERRQ(ierr); 1298 } else if (dimc == 3) { 1299 ierr = DMCreateInjection_DA_3D(dac,daf,&inject);CHKERRQ(ierr); 1300 } 1301 ierr = MatCreateScatter(PetscObjectComm((PetscObject)inject), inject, mat);CHKERRQ(ierr); 1302 ierr = VecScatterDestroy(&inject);CHKERRQ(ierr); 1303 PetscFunctionReturn(0); 1304 } 1305 1306 PetscErrorCode DMCreateAggregates_DA(DM dac,DM daf,Mat *rest) 1307 { 1308 PetscErrorCode ierr; 1309 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc; 1310 PetscInt dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1311 DMBoundaryType bxc,byc,bzc,bxf,byf,bzf; 1312 DMDAStencilType stc,stf; 1313 PetscInt i,j,l; 1314 PetscInt i_start,j_start,l_start, m_f,n_f,p_f; 1315 PetscInt i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost; 1316 const PetscInt *idx_f; 1317 PetscInt i_c,j_c,l_c; 1318 PetscInt i_start_c,j_start_c,l_start_c, m_c,n_c,p_c; 1319 PetscInt i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c; 1320 const PetscInt *idx_c; 1321 PetscInt d; 1322 PetscInt a; 1323 PetscInt max_agg_size; 1324 PetscInt *fine_nodes; 1325 PetscScalar *one_vec; 1326 PetscInt fn_idx; 1327 ISLocalToGlobalMapping ltogmf,ltogmc; 1328 1329 PetscFunctionBegin; 1330 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1331 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1332 PetscValidPointer(rest,3); 1333 1334 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr); 1335 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr); 1336 if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf); 1337 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff); 1338 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf); 1339 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs"); 1340 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs"); 1341 1342 if (Mf < Mc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Mc %D, Mf %D", Mc, Mf); 1343 if (Nf < Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Nc %D, Nf %D", Nc, Nf); 1344 if (Pf < Pc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Coarse grid has more points than fine grid, Pc %D, Pf %D", Pc, Pf); 1345 1346 if (Pc < 0) Pc = 1; 1347 if (Pf < 0) Pf = 1; 1348 if (Nc < 0) Nc = 1; 1349 if (Nf < 0) Nf = 1; 1350 1351 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 1352 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 1353 1354 ierr = DMGetLocalToGlobalMapping(daf,<ogmf);CHKERRQ(ierr); 1355 ierr = ISLocalToGlobalMappingGetIndices(ltogmf,&idx_f);CHKERRQ(ierr); 1356 1357 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 1358 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr); 1359 1360 ierr = DMGetLocalToGlobalMapping(dac,<ogmc);CHKERRQ(ierr); 1361 ierr = ISLocalToGlobalMappingGetIndices(ltogmc,&idx_c);CHKERRQ(ierr); 1362 1363 /* 1364 Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios 1365 for dimension 1 and 2 respectively. 1366 Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1) 1367 and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate. 1368 Each specific dof on the fine grid is mapped to one dof on the coarse grid. 1369 */ 1370 1371 max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1); 1372 1373 /* create the matrix that will contain the restriction operator */ 1374 ierr = MatCreateAIJ(PetscObjectComm((PetscObject)daf), m_c*n_c*p_c*dofc, m_f*n_f*p_f*doff, Mc*Nc*Pc*dofc, Mf*Nf*Pf*doff, 1375 max_agg_size, NULL, max_agg_size, NULL, rest);CHKERRQ(ierr); 1376 1377 /* store nodes in the fine grid here */ 1378 ierr = PetscMalloc2(max_agg_size, &one_vec,max_agg_size, &fine_nodes);CHKERRQ(ierr); 1379 for (i=0; i<max_agg_size; i++) one_vec[i] = 1.0; 1380 1381 /* loop over all coarse nodes */ 1382 for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) { 1383 for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) { 1384 for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) { 1385 for (d=0; d<dofc; d++) { 1386 /* convert to local "natural" numbering and then to PETSc global numbering */ 1387 a = idx_c[dofc*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c))] + d; 1388 1389 fn_idx = 0; 1390 /* Corresponding fine points are all points (i_f, j_f, l_f) such that 1391 i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc 1392 (same for other dimensions) 1393 */ 1394 for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) { 1395 for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) { 1396 for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) { 1397 fine_nodes[fn_idx] = idx_f[doff*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))] + d; 1398 fn_idx++; 1399 } 1400 } 1401 } 1402 /* add all these points to one aggregate */ 1403 ierr = MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);CHKERRQ(ierr); 1404 } 1405 } 1406 } 1407 } 1408 ierr = ISLocalToGlobalMappingRestoreIndices(ltogmf,&idx_f);CHKERRQ(ierr); 1409 ierr = ISLocalToGlobalMappingRestoreIndices(ltogmc,&idx_c);CHKERRQ(ierr); 1410 ierr = PetscFree2(one_vec,fine_nodes);CHKERRQ(ierr); 1411 ierr = MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1412 ierr = MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1413 PetscFunctionReturn(0); 1414 } 1415