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 ratioi = mx/Mx; 508 ratioj = my/My; 509 if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x"); 510 if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y"); 511 if (ratioi != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 2"); 512 if (ratioj != 2) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 2"); 513 514 ierr = DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 515 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 516 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 517 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 518 519 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 520 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 521 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 522 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 523 524 /* 525 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 526 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 527 processors). It's effective length is hence 4 times its normal length, this is 528 why the col_scale is multiplied by the interpolation matrix column sizes. 529 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 530 copy of the coarse vector. A bit of a hack but you do better. 531 532 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 533 */ 534 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);CHKERRQ(ierr); 535 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRQ(ierr); 536 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRQ(ierr); 537 col_scale = size_f/size_c; 538 col_shift = Mx*My*(rank_f/size_c); 539 540 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f,col_scale*m_c*n_c,dnz,onz);CHKERRQ(ierr); 541 for (j=j_start; j<j_start+n_f; j++) { 542 for (i=i_start; i<i_start+m_f; i++) { 543 /* convert to local "natural" numbering and then to PETSc global numbering */ 544 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 545 546 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 547 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 548 549 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\ 550 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 551 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\ 552 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 553 554 /* 555 Only include those interpolation points that are truly 556 nonzero. Note this is very important for final grid lines 557 in x and y directions; since they have no right/top neighbors 558 */ 559 nc = 0; 560 /* one left and below; or we are right on it */ 561 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 562 cols[nc++] = col_shift + idx_c[col]; 563 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 564 } 565 } 566 ierr = MatCreate(PetscObjectComm((PetscObject)daf),&mat);CHKERRQ(ierr); 567 ierr = MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);CHKERRQ(ierr); 568 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 569 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 570 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 571 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 572 573 /* loop over local fine grid nodes setting interpolation for those*/ 574 for (j=j_start; j<j_start+n_f; j++) { 575 for (i=i_start; i<i_start+m_f; i++) { 576 /* convert to local "natural" numbering and then to PETSc global numbering */ 577 row = idx_f[(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 578 579 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 580 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 581 nc = 0; 582 /* one left and below; or we are right on it */ 583 col = (m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 584 cols[nc] = col_shift + idx_c[col]; 585 v[nc++] = 1.0; 586 587 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 588 } 589 } 590 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 591 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 592 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 593 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 594 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 595 ierr = MatDestroy(&mat);CHKERRQ(ierr); 596 ierr = PetscLogFlops(13.0*m_f*n_f);CHKERRQ(ierr); 597 PetscFunctionReturn(0); 598 } 599 600 /* 601 Contributed by Jianming Yang <jianming-yang@uiowa.edu> 602 */ 603 #undef __FUNCT__ 604 #define __FUNCT__ "DMCreateInterpolation_DA_3D_Q0" 605 PetscErrorCode DMCreateInterpolation_DA_3D_Q0(DM dac,DM daf,Mat *A) 606 { 607 PetscErrorCode ierr; 608 PetscInt i,j,l,i_start,j_start,l_start,m_f,n_f,p_f,Mx,My,Mz,dof; 609 const PetscInt *idx_c,*idx_f; 610 ISLocalToGlobalMapping ltog_f,ltog_c; 611 PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c,nc,*dnz,*onz; 612 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; 613 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; 614 PetscMPIInt size_c,size_f,rank_f; 615 PetscScalar v[8]; 616 Mat mat; 617 DMBoundaryType bx,by,bz; 618 619 PetscFunctionBegin; 620 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 621 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 622 ratioi = mx/Mx; 623 ratioj = my/My; 624 ratiol = mz/Mz; 625 if (ratioi*Mx != mx) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x"); 626 if (ratioj*My != my) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y"); 627 if (ratiol*Mz != mz) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in z"); 628 if (ratioi != 2 && ratioi != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 1 or 2"); 629 if (ratioj != 2 && ratioj != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 1 or 2"); 630 if (ratiol != 2 && ratiol != 1) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_WRONG,"Coarsening factor in z must be 1 or 2"); 631 632 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 633 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 634 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 635 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 636 637 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 638 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); 639 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 640 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 641 642 /* 643 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 644 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 645 processors). It's effective length is hence 4 times its normal length, this is 646 why the col_scale is multiplied by the interpolation matrix column sizes. 647 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 648 copy of the coarse vector. A bit of a hack but you do better. 649 650 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 651 */ 652 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dac),&size_c);CHKERRQ(ierr); 653 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)daf),&size_f);CHKERRQ(ierr); 654 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)daf),&rank_f);CHKERRQ(ierr); 655 col_scale = size_f/size_c; 656 col_shift = Mx*My*Mz*(rank_f/size_c); 657 658 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)daf),m_f*n_f*p_f,col_scale*m_c*n_c*p_c,dnz,onz);CHKERRQ(ierr); 659 for (l=l_start; l<l_start+p_f; l++) { 660 for (j=j_start; j<j_start+n_f; j++) { 661 for (i=i_start; i<i_start+m_f; i++) { 662 /* convert to local "natural" numbering and then to PETSc global numbering */ 663 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 664 665 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 666 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 667 l_c = (l/ratiol); 668 669 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\ 670 l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c); 671 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\ 672 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 673 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\ 674 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 675 676 /* 677 Only include those interpolation points that are truly 678 nonzero. Note this is very important for final grid lines 679 in x and y directions; since they have no right/top neighbors 680 */ 681 nc = 0; 682 /* one left and below; or we are right on it */ 683 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)); 684 cols[nc++] = col_shift + idx_c[col]; 685 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 686 } 687 } 688 } 689 ierr = MatCreate(PetscObjectComm((PetscObject)daf),&mat);CHKERRQ(ierr); 690 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); 691 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 692 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 693 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 694 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 695 696 /* loop over local fine grid nodes setting interpolation for those*/ 697 for (l=l_start; l<l_start+p_f; l++) { 698 for (j=j_start; j<j_start+n_f; j++) { 699 for (i=i_start; i<i_start+m_f; i++) { 700 /* convert to local "natural" numbering and then to PETSc global numbering */ 701 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 702 703 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 704 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 705 l_c = (l/ratiol); 706 nc = 0; 707 /* one left and below; or we are right on it */ 708 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)); 709 cols[nc] = col_shift + idx_c[col]; 710 v[nc++] = 1.0; 711 712 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 713 } 714 } 715 } 716 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 717 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 718 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 719 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 720 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 721 ierr = MatDestroy(&mat);CHKERRQ(ierr); 722 ierr = PetscLogFlops(13.0*m_f*n_f*p_f);CHKERRQ(ierr); 723 PetscFunctionReturn(0); 724 } 725 726 #undef __FUNCT__ 727 #define __FUNCT__ "DMCreateInterpolation_DA_3D_Q1" 728 PetscErrorCode DMCreateInterpolation_DA_3D_Q1(DM dac,DM daf,Mat *A) 729 { 730 PetscErrorCode ierr; 731 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof,l; 732 const PetscInt *idx_c,*idx_f; 733 ISLocalToGlobalMapping ltog_f,ltog_c; 734 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c,Mz,mz; 735 PetscInt row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok; 736 PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c; 737 PetscInt l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c; 738 PetscInt l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz; 739 PetscScalar v[8],x,y,z; 740 Mat mat; 741 DMBoundaryType bx,by,bz; 742 743 PetscFunctionBegin; 744 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 745 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 746 if (mx == Mx) { 747 ratioi = 1; 748 } else if (bx == DM_BOUNDARY_PERIODIC) { 749 ratioi = mx/Mx; 750 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); 751 } else { 752 ratioi = (mx-1)/(Mx-1); 753 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); 754 } 755 if (my == My) { 756 ratioj = 1; 757 } else if (by == DM_BOUNDARY_PERIODIC) { 758 ratioj = my/My; 759 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); 760 } else { 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 ratiok = mz/Mz; 768 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); 769 } else { 770 ratiok = (mz-1)/(Mz-1); 771 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); 772 } 773 774 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 775 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 776 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 777 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 778 779 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 780 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); 781 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 782 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 783 784 /* create interpolation matrix, determining exact preallocation */ 785 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)dac),m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);CHKERRQ(ierr); 786 /* loop over local fine grid nodes counting interpolating points */ 787 for (l=l_start; l<l_start+p_f; l++) { 788 for (j=j_start; j<j_start+n_f; j++) { 789 for (i=i_start; i<i_start+m_f; i++) { 790 /* convert to local "natural" numbering and then to PETSc global numbering */ 791 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 792 i_c = (i/ratioi); 793 j_c = (j/ratioj); 794 l_c = (l/ratiok); 795 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\ 796 l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c); 797 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\ 798 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 799 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\ 800 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 801 802 /* 803 Only include those interpolation points that are truly 804 nonzero. Note this is very important for final grid lines 805 in x and y directions; since they have no right/top neighbors 806 */ 807 nc = 0; 808 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)); 809 cols[nc++] = idx_c[col]; 810 if (i_c*ratioi != i) { 811 cols[nc++] = idx_c[col+1]; 812 } 813 if (j_c*ratioj != j) { 814 cols[nc++] = idx_c[col+m_ghost_c]; 815 } 816 if (l_c*ratiok != l) { 817 cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c]; 818 } 819 if (j_c*ratioj != j && i_c*ratioi != i) { 820 cols[nc++] = idx_c[col+(m_ghost_c+1)]; 821 } 822 if (j_c*ratioj != j && l_c*ratiok != l) { 823 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)]; 824 } 825 if (i_c*ratioi != i && l_c*ratiok != l) { 826 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)]; 827 } 828 if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) { 829 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)]; 830 } 831 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 832 } 833 } 834 } 835 ierr = MatCreate(PetscObjectComm((PetscObject)dac),&mat);CHKERRQ(ierr); 836 ierr = MatSetSizes(mat,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz);CHKERRQ(ierr); 837 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 838 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 839 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 840 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 841 842 /* loop over local fine grid nodes setting interpolation for those*/ 843 if (!NEWVERSION) { 844 845 for (l=l_start; l<l_start+p_f; l++) { 846 for (j=j_start; j<j_start+n_f; j++) { 847 for (i=i_start; i<i_start+m_f; i++) { 848 /* convert to local "natural" numbering and then to PETSc global numbering */ 849 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 850 851 i_c = (i/ratioi); 852 j_c = (j/ratioj); 853 l_c = (l/ratiok); 854 855 /* 856 Only include those interpolation points that are truly 857 nonzero. Note this is very important for final grid lines 858 in x and y directions; since they have no right/top neighbors 859 */ 860 x = ((PetscReal)(i - i_c*ratioi))/((PetscReal)ratioi); 861 y = ((PetscReal)(j - j_c*ratioj))/((PetscReal)ratioj); 862 z = ((PetscReal)(l - l_c*ratiok))/((PetscReal)ratiok); 863 864 nc = 0; 865 /* one left and below; or we are right on it */ 866 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)); 867 868 cols[nc] = idx_c[col]; 869 v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.)); 870 871 if (i_c*ratioi != i) { 872 cols[nc] = idx_c[col+1]; 873 v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.)); 874 } 875 876 if (j_c*ratioj != j) { 877 cols[nc] = idx_c[col+m_ghost_c]; 878 v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.)); 879 } 880 881 if (l_c*ratiok != l) { 882 cols[nc] = idx_c[col+m_ghost_c*n_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 (j_c*ratioj != j && i_c*ratioi != i) { 887 cols[nc] = idx_c[col+(m_ghost_c+1)]; 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 && l_c*ratiok != l) { 892 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)]; 893 v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.)); 894 } 895 896 if (i_c*ratioi != i && l_c*ratiok != l) { 897 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)]; 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 && j_c*ratioj != j) { 902 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)]; 903 v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.)); 904 } 905 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 906 } 907 } 908 } 909 910 } else { 911 PetscScalar *xi,*eta,*zeta; 912 PetscInt li,nxi,lj,neta,lk,nzeta,n; 913 PetscScalar Ni[8]; 914 915 /* compute local coordinate arrays */ 916 nxi = ratioi + 1; 917 neta = ratioj + 1; 918 nzeta = ratiok + 1; 919 ierr = PetscMalloc1(nxi,&xi);CHKERRQ(ierr); 920 ierr = PetscMalloc1(neta,&eta);CHKERRQ(ierr); 921 ierr = PetscMalloc1(nzeta,&zeta);CHKERRQ(ierr); 922 for (li=0; li<nxi; li++) xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1)); 923 for (lj=0; lj<neta; lj++) eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1)); 924 for (lk=0; lk<nzeta; lk++) zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1)); 925 926 for (l=l_start; l<l_start+p_f; l++) { 927 for (j=j_start; j<j_start+n_f; j++) { 928 for (i=i_start; i<i_start+m_f; i++) { 929 /* convert to local "natural" numbering and then to PETSc global numbering */ 930 row = idx_f[(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 931 932 i_c = (i/ratioi); 933 j_c = (j/ratioj); 934 l_c = (l/ratiok); 935 936 /* remainders */ 937 li = i - ratioi * (i/ratioi); 938 if (i==mx-1) li = nxi-1; 939 lj = j - ratioj * (j/ratioj); 940 if (j==my-1) lj = neta-1; 941 lk = l - ratiok * (l/ratiok); 942 if (l==mz-1) lk = nzeta-1; 943 944 /* corners */ 945 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)); 946 cols[0] = idx_c[col]; 947 Ni[0] = 1.0; 948 if ((li==0) || (li==nxi-1)) { 949 if ((lj==0) || (lj==neta-1)) { 950 if ((lk==0) || (lk==nzeta-1)) { 951 ierr = MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr); 952 continue; 953 } 954 } 955 } 956 957 /* edges + interior */ 958 /* remainders */ 959 if (i==mx-1) i_c--; 960 if (j==my-1) j_c--; 961 if (l==mz-1) l_c--; 962 963 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)); 964 cols[0] = idx_c[col]; /* one left and below; or we are right on it */ 965 cols[1] = idx_c[col+1]; /* one right and below */ 966 cols[2] = idx_c[col+m_ghost_c]; /* one left and above */ 967 cols[3] = idx_c[col+(m_ghost_c+1)]; /* one right and above */ 968 969 cols[4] = idx_c[col+m_ghost_c*n_ghost_c]; /* one left and below and front; or we are right on it */ 970 cols[5] = idx_c[col+(m_ghost_c*n_ghost_c+1)]; /* one right and below, and front */ 971 cols[6] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)]; /* one left and above and front*/ 972 cols[7] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)]; /* one right and above and front */ 973 974 Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]); 975 Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]); 976 Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]); 977 Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]); 978 979 Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]); 980 Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]); 981 Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]); 982 Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]); 983 984 for (n=0; n<8; n++) { 985 if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1; 986 } 987 ierr = MatSetValues(mat,1,&row,8,cols,Ni,INSERT_VALUES);CHKERRQ(ierr); 988 989 } 990 } 991 } 992 ierr = PetscFree(xi);CHKERRQ(ierr); 993 ierr = PetscFree(eta);CHKERRQ(ierr); 994 ierr = PetscFree(zeta);CHKERRQ(ierr); 995 } 996 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 997 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 998 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 999 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1000 1001 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 1002 ierr = MatDestroy(&mat);CHKERRQ(ierr); 1003 PetscFunctionReturn(0); 1004 } 1005 1006 #undef __FUNCT__ 1007 #define __FUNCT__ "DMCreateInterpolation_DA" 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);CHKERRQ(ierr); 1025 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);CHKERRQ(ierr); 1026 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);CHKERRQ(ierr); 1027 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");CHKERRQ(ierr); 1028 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");CHKERRQ(ierr); 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 #undef __FUNCT__ 1057 #define __FUNCT__ "DMCreateInjection_DA_1D" 1058 PetscErrorCode DMCreateInjection_DA_1D(DM dac,DM daf,VecScatter *inject) 1059 { 1060 PetscErrorCode ierr; 1061 PetscInt i,i_start,m_f,Mx,dof; 1062 const PetscInt *idx_f; 1063 ISLocalToGlobalMapping ltog_f; 1064 PetscInt m_ghost,m_ghost_c; 1065 PetscInt row,i_start_ghost,mx,m_c,nc,ratioi; 1066 PetscInt i_start_c,i_start_ghost_c; 1067 PetscInt *cols; 1068 DMBoundaryType bx; 1069 Vec vecf,vecc; 1070 IS isf; 1071 1072 PetscFunctionBegin; 1073 ierr = DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);CHKERRQ(ierr); 1074 ierr = DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 1075 if (bx == DM_BOUNDARY_PERIODIC) { 1076 ratioi = mx/Mx; 1077 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); 1078 } else { 1079 ratioi = (mx-1)/(Mx-1); 1080 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); 1081 } 1082 1083 ierr = DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);CHKERRQ(ierr); 1084 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);CHKERRQ(ierr); 1085 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 1086 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1087 1088 ierr = DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);CHKERRQ(ierr); 1089 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);CHKERRQ(ierr); 1090 1091 1092 /* loop over local fine grid nodes setting interpolation for those*/ 1093 nc = 0; 1094 ierr = PetscMalloc1(m_f,&cols);CHKERRQ(ierr); 1095 1096 1097 for (i=i_start_c; i<i_start_c+m_c; i++) { 1098 PetscInt i_f = i*ratioi; 1099 1100 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); 1101 1102 row = idx_f[(i_f-i_start_ghost)]; 1103 cols[nc++] = row; 1104 } 1105 1106 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1107 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1108 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1109 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1110 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr); 1111 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1112 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1113 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1114 PetscFunctionReturn(0); 1115 } 1116 1117 #undef __FUNCT__ 1118 #define __FUNCT__ "DMCreateInjection_DA_2D" 1119 PetscErrorCode DMCreateInjection_DA_2D(DM dac,DM daf,VecScatter *inject) 1120 { 1121 PetscErrorCode ierr; 1122 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,dof; 1123 const PetscInt *idx_c,*idx_f; 1124 ISLocalToGlobalMapping ltog_f,ltog_c; 1125 PetscInt m_ghost,n_ghost,m_ghost_c,n_ghost_c; 1126 PetscInt row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj; 1127 PetscInt i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c; 1128 PetscInt *cols; 1129 DMBoundaryType bx,by; 1130 Vec vecf,vecc; 1131 IS isf; 1132 1133 PetscFunctionBegin; 1134 ierr = DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);CHKERRQ(ierr); 1135 ierr = DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 1136 if (bx == DM_BOUNDARY_PERIODIC) { 1137 ratioi = mx/Mx; 1138 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); 1139 } else { 1140 ratioi = (mx-1)/(Mx-1); 1141 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); 1142 } 1143 if (by == DM_BOUNDARY_PERIODIC) { 1144 ratioj = my/My; 1145 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); 1146 } else { 1147 ratioj = (my-1)/(My-1); 1148 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); 1149 } 1150 1151 ierr = DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 1152 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 1153 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 1154 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1155 1156 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 1157 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 1158 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 1159 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1160 1161 /* loop over local fine grid nodes setting interpolation for those*/ 1162 nc = 0; 1163 ierr = PetscMalloc1(n_f*m_f,&cols);CHKERRQ(ierr); 1164 for (j=j_start_c; j<j_start_c+n_c; j++) { 1165 for (i=i_start_c; i<i_start_c+m_c; i++) { 1166 PetscInt i_f = i*ratioi,j_f = j*ratioj; 1167 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\ 1168 j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost); 1169 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\ 1170 i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1171 row = idx_f[(m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))]; 1172 cols[nc++] = row; 1173 } 1174 } 1175 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1176 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1177 1178 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1179 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1180 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1181 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr); 1182 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1183 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1184 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1185 PetscFunctionReturn(0); 1186 } 1187 1188 #undef __FUNCT__ 1189 #define __FUNCT__ "DMCreateInjection_DA_3D" 1190 PetscErrorCode DMCreateInjection_DA_3D(DM dac,DM daf,VecScatter *inject) 1191 { 1192 PetscErrorCode ierr; 1193 PetscInt i,j,k,i_start,j_start,k_start,m_f,n_f,p_f,Mx,My,Mz; 1194 PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c; 1195 PetscInt i_start_ghost,j_start_ghost,k_start_ghost; 1196 PetscInt mx,my,mz,ratioi,ratioj,ratiok; 1197 PetscInt i_start_c,j_start_c,k_start_c; 1198 PetscInt m_c,n_c,p_c; 1199 PetscInt i_start_ghost_c,j_start_ghost_c,k_start_ghost_c; 1200 PetscInt row,nc,dof; 1201 const PetscInt *idx_c,*idx_f; 1202 ISLocalToGlobalMapping ltog_f,ltog_c; 1203 PetscInt *cols; 1204 DMBoundaryType bx,by,bz; 1205 Vec vecf,vecc; 1206 IS isf; 1207 1208 PetscFunctionBegin; 1209 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 1210 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 1211 1212 if (bx == DM_BOUNDARY_PERIODIC) { 1213 ratioi = mx/Mx; 1214 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); 1215 } else { 1216 ratioi = (mx-1)/(Mx-1); 1217 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); 1218 } 1219 if (by == DM_BOUNDARY_PERIODIC) { 1220 ratioj = my/My; 1221 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); 1222 } else { 1223 ratioj = (my-1)/(My-1); 1224 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); 1225 } 1226 if (bz == DM_BOUNDARY_PERIODIC) { 1227 ratiok = mz/Mz; 1228 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); 1229 } else { 1230 ratiok = (mz-1)/(Mz-1); 1231 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); 1232 } 1233 1234 ierr = DMDAGetCorners(daf,&i_start,&j_start,&k_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 1235 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&k_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 1236 ierr = DMGetLocalToGlobalMapping(daf,<og_f);CHKERRQ(ierr); 1237 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1238 1239 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&k_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 1240 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); 1241 ierr = DMGetLocalToGlobalMapping(dac,<og_c);CHKERRQ(ierr); 1242 ierr = ISLocalToGlobalMappingGetBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1243 1244 1245 /* loop over local fine grid nodes setting interpolation for those*/ 1246 nc = 0; 1247 ierr = PetscMalloc1(n_f*m_f*p_f,&cols);CHKERRQ(ierr); 1248 for (k=k_start_c; k<k_start_c+p_c; k++) { 1249 for (j=j_start_c; j<j_start_c+n_c; j++) { 1250 for (i=i_start_c; i<i_start_c+m_c; i++) { 1251 PetscInt i_f = i*ratioi,j_f = j*ratioj,k_f = k*ratiok; 1252 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 " 1253 "k_c %D k_f %D fine ghost range [%D,%D]",k,k_f,k_start_ghost,k_start_ghost+p_ghost); 1254 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 " 1255 "j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost); 1256 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 " 1257 "i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1258 row = idx_f[(m_ghost*n_ghost*(k_f-k_start_ghost) + m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))]; 1259 cols[nc++] = row; 1260 } 1261 } 1262 } 1263 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_f,&idx_f);CHKERRQ(ierr); 1264 ierr = ISLocalToGlobalMappingRestoreBlockIndices(ltog_c,&idx_c);CHKERRQ(ierr); 1265 1266 ierr = ISCreateBlock(PetscObjectComm((PetscObject)daf),dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1267 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1268 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1269 ierr = VecScatterCreate(vecf,isf,vecc,NULL,inject);CHKERRQ(ierr); 1270 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1271 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1272 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1273 PetscFunctionReturn(0); 1274 } 1275 1276 #undef __FUNCT__ 1277 #define __FUNCT__ "DMCreateInjection_DA" 1278 PetscErrorCode DMCreateInjection_DA(DM dac,DM daf,Mat *mat) 1279 { 1280 PetscErrorCode ierr; 1281 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1282 DMBoundaryType bxc,byc,bzc,bxf,byf,bzf; 1283 DMDAStencilType stc,stf; 1284 VecScatter inject = NULL; 1285 1286 PetscFunctionBegin; 1287 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1288 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1289 PetscValidPointer(mat,3); 1290 1291 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr); 1292 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr); 1293 if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 1294 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);CHKERRQ(ierr); 1295 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);CHKERRQ(ierr); 1296 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");CHKERRQ(ierr); 1297 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");CHKERRQ(ierr); 1298 if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); 1299 if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction"); 1300 if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction"); 1301 1302 if (dimc == 1) { 1303 ierr = DMCreateInjection_DA_1D(dac,daf,&inject);CHKERRQ(ierr); 1304 } else if (dimc == 2) { 1305 ierr = DMCreateInjection_DA_2D(dac,daf,&inject);CHKERRQ(ierr); 1306 } else if (dimc == 3) { 1307 ierr = DMCreateInjection_DA_3D(dac,daf,&inject);CHKERRQ(ierr); 1308 } 1309 ierr = MatCreateScatter(PetscObjectComm((PetscObject)inject), inject, mat);CHKERRQ(ierr); 1310 ierr = VecScatterDestroy(&inject);CHKERRQ(ierr); 1311 PetscFunctionReturn(0); 1312 } 1313 1314 #undef __FUNCT__ 1315 #define __FUNCT__ "DMCreateAggregates_DA" 1316 PetscErrorCode DMCreateAggregates_DA(DM dac,DM daf,Mat *rest) 1317 { 1318 PetscErrorCode ierr; 1319 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc; 1320 PetscInt dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1321 DMBoundaryType bxc,byc,bzc,bxf,byf,bzf; 1322 DMDAStencilType stc,stf; 1323 PetscInt i,j,l; 1324 PetscInt i_start,j_start,l_start, m_f,n_f,p_f; 1325 PetscInt i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost; 1326 const PetscInt *idx_f; 1327 PetscInt i_c,j_c,l_c; 1328 PetscInt i_start_c,j_start_c,l_start_c, m_c,n_c,p_c; 1329 PetscInt i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c; 1330 const PetscInt *idx_c; 1331 PetscInt d; 1332 PetscInt a; 1333 PetscInt max_agg_size; 1334 PetscInt *fine_nodes; 1335 PetscScalar *one_vec; 1336 PetscInt fn_idx; 1337 ISLocalToGlobalMapping ltogmf,ltogmc; 1338 1339 PetscFunctionBegin; 1340 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1341 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1342 PetscValidPointer(rest,3); 1343 1344 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr); 1345 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr); 1346 if (dimc != dimf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 1347 if (dofc != doff) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);CHKERRQ(ierr); 1348 if (sc != sf) SETERRQ2(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);CHKERRQ(ierr); 1349 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");CHKERRQ(ierr); 1350 if (stc != stf) SETERRQ(PetscObjectComm((PetscObject)daf),PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");CHKERRQ(ierr); 1351 1352 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); 1353 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); 1354 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); 1355 1356 if (Pc < 0) Pc = 1; 1357 if (Pf < 0) Pf = 1; 1358 if (Nc < 0) Nc = 1; 1359 if (Nf < 0) Nf = 1; 1360 1361 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 1362 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 1363 1364 ierr = DMGetLocalToGlobalMapping(daf,<ogmf);CHKERRQ(ierr); 1365 ierr = ISLocalToGlobalMappingGetIndices(ltogmf,&idx_f);CHKERRQ(ierr); 1366 1367 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 1368 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); 1369 1370 ierr = DMGetLocalToGlobalMapping(dac,<ogmc);CHKERRQ(ierr); 1371 ierr = ISLocalToGlobalMappingGetIndices(ltogmc,&idx_c);CHKERRQ(ierr); 1372 1373 /* 1374 Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios 1375 for dimension 1 and 2 respectively. 1376 Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1) 1377 and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate. 1378 Each specific dof on the fine grid is mapped to one dof on the coarse grid. 1379 */ 1380 1381 max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1); 1382 1383 /* create the matrix that will contain the restriction operator */ 1384 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, 1385 max_agg_size, NULL, max_agg_size, NULL, rest);CHKERRQ(ierr); 1386 1387 /* store nodes in the fine grid here */ 1388 ierr = PetscMalloc2(max_agg_size, &one_vec,max_agg_size, &fine_nodes);CHKERRQ(ierr); 1389 for (i=0; i<max_agg_size; i++) one_vec[i] = 1.0; 1390 1391 /* loop over all coarse nodes */ 1392 for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) { 1393 for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) { 1394 for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) { 1395 for (d=0; d<dofc; d++) { 1396 /* convert to local "natural" numbering and then to PETSc global numbering */ 1397 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; 1398 1399 fn_idx = 0; 1400 /* Corresponding fine points are all points (i_f, j_f, l_f) such that 1401 i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc 1402 (same for other dimensions) 1403 */ 1404 for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) { 1405 for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) { 1406 for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) { 1407 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; 1408 fn_idx++; 1409 } 1410 } 1411 } 1412 /* add all these points to one aggregate */ 1413 ierr = MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);CHKERRQ(ierr); 1414 } 1415 } 1416 } 1417 } 1418 ierr = ISLocalToGlobalMappingRestoreIndices(ltogmf,&idx_f);CHKERRQ(ierr); 1419 ierr = ISLocalToGlobalMappingRestoreIndices(ltogmc,&idx_c);CHKERRQ(ierr); 1420 ierr = PetscFree2(one_vec,fine_nodes);CHKERRQ(ierr); 1421 ierr = MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1422 ierr = MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1423 PetscFunctionReturn(0); 1424 } 1425