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