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