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