1 #define PETSCDM_DLL 2 3 /* 4 Code for interpolating between grids represented by DAs 5 */ 6 7 #include "private/daimpl.h" /*I "petscda.h" I*/ 8 #include "petscmg.h" 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "DMGetInterpolationScale" 12 /*@ 13 DMGetInterpolationScale - Forms L = R*1/diag(R*1) - L.*v is like a coarse grid average of the 14 nearby fine grid points. 15 16 Input Parameters: 17 + dac - DM that defines a coarse mesh 18 . daf - DM that defines a fine mesh 19 - mat - the restriction (or interpolation operator) from fine to coarse 20 21 Output Parameter: 22 . scale - the scaled vector 23 24 Level: developer 25 26 .seealso: DMGetInterpolation() 27 28 @*/ 29 PetscErrorCode PETSCDM_DLLEXPORT DMGetInterpolationScale(DM dac,DM daf,Mat mat,Vec *scale) 30 { 31 PetscErrorCode ierr; 32 Vec fine; 33 PetscScalar one = 1.0; 34 35 PetscFunctionBegin; 36 ierr = DMCreateGlobalVector(daf,&fine);CHKERRQ(ierr); 37 ierr = DMCreateGlobalVector(dac,scale);CHKERRQ(ierr); 38 ierr = VecSet(fine,one);CHKERRQ(ierr); 39 ierr = MatRestrict(mat,fine,*scale);CHKERRQ(ierr); 40 ierr = VecDestroy(fine);CHKERRQ(ierr); 41 ierr = VecReciprocal(*scale);CHKERRQ(ierr); 42 PetscFunctionReturn(0); 43 } 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "DAGetInterpolation_1D_Q1" 47 PetscErrorCode DAGetInterpolation_1D_Q1(DM dac,DM daf,Mat *A) 48 { 49 PetscErrorCode ierr; 50 PetscInt i,i_start,m_f,Mx,*idx_f; 51 PetscInt m_ghost,*idx_c,m_ghost_c; 52 PetscInt row,col,i_start_ghost,mx,m_c,nc,ratio; 53 PetscInt i_c,i_start_c,i_start_ghost_c,cols[2],dof; 54 PetscScalar v[2],x,*coors = 0,*ccoors; 55 Mat mat; 56 DAPeriodicType pt; 57 Vec vcoors,cvcoors; 58 DM_DA *ddc = (DM_DA*)dac->data, *ddf = (DM_DA*)daf->data; 59 60 PetscFunctionBegin; 61 ierr = DAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 62 ierr = DAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 63 if (pt == DA_XPERIODIC) { 64 ratio = mx/Mx; 65 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); 66 } else { 67 ratio = (mx-1)/(Mx-1); 68 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); 69 } 70 71 ierr = DAGetCorners(daf,&i_start,0,0,&m_f,0,0);CHKERRQ(ierr); 72 ierr = DAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);CHKERRQ(ierr); 73 ierr = DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 74 75 ierr = DAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);CHKERRQ(ierr); 76 ierr = DAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);CHKERRQ(ierr); 77 ierr = DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 78 79 /* create interpolation matrix */ 80 ierr = MatCreate(((PetscObject)dac)->comm,&mat);CHKERRQ(ierr); 81 ierr = MatSetSizes(mat,m_f,m_c,mx,Mx);CHKERRQ(ierr); 82 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 83 ierr = MatSeqAIJSetPreallocation(mat,2,PETSC_NULL);CHKERRQ(ierr); 84 ierr = MatMPIAIJSetPreallocation(mat,2,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 85 86 ierr = DAGetCoordinates(daf,&vcoors);CHKERRQ(ierr); 87 if (vcoors) { 88 ierr = DAGetGhostedCoordinates(dac,&cvcoors);CHKERRQ(ierr); 89 ierr = DAVecGetArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 90 ierr = DAVecGetArray(ddc->da_coordinates,cvcoors,&ccoors);CHKERRQ(ierr); 91 } 92 /* loop over local fine grid nodes setting interpolation for those*/ 93 for (i=i_start; i<i_start+m_f; i++) { 94 /* convert to local "natural" numbering and then to PETSc global numbering */ 95 row = idx_f[dof*(i-i_start_ghost)]/dof; 96 97 i_c = (i/ratio); /* coarse grid node to left of fine grid node */ 98 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 99 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 100 101 /* 102 Only include those interpolation points that are truly 103 nonzero. Note this is very important for final grid lines 104 in x direction; since they have no right neighbor 105 */ 106 if (coors) { 107 x = (coors[i] - ccoors[i_c]); 108 /* only access the next coors point if we know there is one */ 109 /* note this is dangerous because x may not exactly equal ZERO */ 110 if (PetscAbsScalar(x) != 0.0) x = x/(ccoors[i_c+1] - ccoors[i_c]); 111 } else { 112 x = ((double)(i - i_c*ratio))/((double)ratio); 113 } 114 nc = 0; 115 /* one left and below; or we are right on it */ 116 col = dof*(i_c-i_start_ghost_c); 117 cols[nc] = idx_c[col]/dof; 118 v[nc++] = - x + 1.0; 119 /* one right? */ 120 if (i_c*ratio != i) { 121 cols[nc] = idx_c[col+dof]/dof; 122 v[nc++] = x; 123 } 124 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 125 } 126 if (vcoors) { 127 ierr = DAVecRestoreArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 128 ierr = DAVecRestoreArray(ddc->da_coordinates,cvcoors,&ccoors);CHKERRQ(ierr); 129 } 130 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 131 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 132 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 133 ierr = MatDestroy(mat);CHKERRQ(ierr); 134 ierr = PetscLogFlops(5.0*m_f);CHKERRQ(ierr); 135 PetscFunctionReturn(0); 136 } 137 138 #undef __FUNCT__ 139 #define __FUNCT__ "DAGetInterpolation_1D_Q0" 140 PetscErrorCode DAGetInterpolation_1D_Q0(DM dac,DM daf,Mat *A) 141 { 142 PetscErrorCode ierr; 143 PetscInt i,i_start,m_f,Mx,*idx_f; 144 PetscInt m_ghost,*idx_c,m_ghost_c; 145 PetscInt row,col,i_start_ghost,mx,m_c,nc,ratio; 146 PetscInt i_c,i_start_c,i_start_ghost_c,cols[2],dof; 147 PetscScalar v[2],x; 148 Mat mat; 149 DAPeriodicType pt; 150 151 PetscFunctionBegin; 152 ierr = DAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 153 ierr = DAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 154 if (pt == DA_XPERIODIC) { 155 ratio = mx/Mx; 156 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); 157 } else { 158 ratio = (mx-1)/(Mx-1); 159 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); 160 } 161 162 ierr = DAGetCorners(daf,&i_start,0,0,&m_f,0,0);CHKERRQ(ierr); 163 ierr = DAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);CHKERRQ(ierr); 164 ierr = DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 165 166 ierr = DAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);CHKERRQ(ierr); 167 ierr = DAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);CHKERRQ(ierr); 168 ierr = DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 169 170 /* create interpolation matrix */ 171 ierr = MatCreate(((PetscObject)dac)->comm,&mat);CHKERRQ(ierr); 172 ierr = MatSetSizes(mat,m_f,m_c,mx,Mx);CHKERRQ(ierr); 173 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 174 ierr = MatSeqAIJSetPreallocation(mat,2,PETSC_NULL);CHKERRQ(ierr); 175 ierr = MatMPIAIJSetPreallocation(mat,2,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 176 177 /* loop over local fine grid nodes setting interpolation for those*/ 178 for (i=i_start; i<i_start+m_f; i++) { 179 /* convert to local "natural" numbering and then to PETSc global numbering */ 180 row = idx_f[dof*(i-i_start_ghost)]/dof; 181 182 i_c = (i/ratio); /* coarse grid node to left of fine grid node */ 183 184 /* 185 Only include those interpolation points that are truly 186 nonzero. Note this is very important for final grid lines 187 in x direction; since they have no right neighbor 188 */ 189 x = ((double)(i - i_c*ratio))/((double)ratio); 190 nc = 0; 191 /* one left and below; or we are right on it */ 192 col = dof*(i_c-i_start_ghost_c); 193 cols[nc] = idx_c[col]/dof; 194 v[nc++] = - x + 1.0; 195 /* one right? */ 196 if (i_c*ratio != i) { 197 cols[nc] = idx_c[col+dof]/dof; 198 v[nc++] = x; 199 } 200 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 201 } 202 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 203 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 204 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 205 ierr = MatDestroy(mat);CHKERRQ(ierr); 206 ierr = PetscLogFlops(5.0*m_f);CHKERRQ(ierr); 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "DAGetInterpolation_2D_Q1" 212 PetscErrorCode DAGetInterpolation_2D_Q1(DM dac,DM daf,Mat *A) 213 { 214 PetscErrorCode ierr; 215 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof; 216 PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,*dnz,*onz; 217 PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj; 218 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; 219 PetscMPIInt size_c,size_f,rank_f; 220 PetscScalar v[4],x,y; 221 Mat mat; 222 DAPeriodicType pt; 223 DACoor2d **coors = 0,**ccoors; 224 Vec vcoors,cvcoors; 225 DM_DA *ddc = (DM_DA*)dac->data, *ddf = (DM_DA*)daf->data; 226 227 PetscFunctionBegin; 228 ierr = DAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 229 ierr = DAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 230 if (DAXPeriodic(pt)){ 231 ratioi = mx/Mx; 232 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); 233 } else { 234 ratioi = (mx-1)/(Mx-1); 235 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); 236 } 237 if (DAYPeriodic(pt)){ 238 ratioj = my/My; 239 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); 240 } else { 241 ratioj = (my-1)/(My-1); 242 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); 243 } 244 245 246 ierr = DAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 247 ierr = DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 248 ierr = DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 249 250 ierr = DAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 251 ierr = DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 252 ierr = DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 253 254 /* 255 Used for handling a coarse DA that lives on 1/4 the processors of the fine DA. 256 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 257 processors). It's effective length is hence 4 times its normal length, this is 258 why the col_scale is multiplied by the interpolation matrix column sizes. 259 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 260 copy of the coarse vector. A bit of a hack but you do better. 261 262 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 263 */ 264 ierr = MPI_Comm_size(((PetscObject)dac)->comm,&size_c);CHKERRQ(ierr); 265 ierr = MPI_Comm_size(((PetscObject)daf)->comm,&size_f);CHKERRQ(ierr); 266 ierr = MPI_Comm_rank(((PetscObject)daf)->comm,&rank_f);CHKERRQ(ierr); 267 col_scale = size_f/size_c; 268 col_shift = Mx*My*(rank_f/size_c); 269 270 ierr = MatPreallocateInitialize(((PetscObject)daf)->comm,m_f*n_f,col_scale*m_c*n_c,dnz,onz);CHKERRQ(ierr); 271 for (j=j_start; j<j_start+n_f; j++) { 272 for (i=i_start; i<i_start+m_f; i++) { 273 /* convert to local "natural" numbering and then to PETSc global numbering */ 274 row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 275 276 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 277 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 278 279 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 280 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 281 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 282 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 283 284 /* 285 Only include those interpolation points that are truly 286 nonzero. Note this is very important for final grid lines 287 in x and y directions; since they have no right/top neighbors 288 */ 289 nc = 0; 290 /* one left and below; or we are right on it */ 291 col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 292 cols[nc++] = col_shift + idx_c[col]/dof; 293 /* one right and below */ 294 if (i_c*ratioi != i) { 295 cols[nc++] = col_shift + idx_c[col+dof]/dof; 296 } 297 /* one left and above */ 298 if (j_c*ratioj != j) { 299 cols[nc++] = col_shift + idx_c[col+m_ghost_c*dof]/dof; 300 } 301 /* one right and above */ 302 if (j_c*ratioi != j && i_c*ratioj != i) { 303 cols[nc++] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof; 304 } 305 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 306 } 307 } 308 ierr = MatCreate(((PetscObject)daf)->comm,&mat);CHKERRQ(ierr); 309 ierr = MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);CHKERRQ(ierr); 310 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 311 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 312 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 313 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 314 315 ierr = DAGetCoordinates(daf,&vcoors);CHKERRQ(ierr); 316 if (vcoors) { 317 ierr = DAGetGhostedCoordinates(dac,&cvcoors);CHKERRQ(ierr); 318 ierr = DAVecGetArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 319 ierr = DAVecGetArray(ddc->da_coordinates,cvcoors,&ccoors);CHKERRQ(ierr); 320 } 321 322 /* loop over local fine grid nodes setting interpolation for those*/ 323 for (j=j_start; j<j_start+n_f; j++) { 324 for (i=i_start; i<i_start+m_f; i++) { 325 /* convert to local "natural" numbering and then to PETSc global numbering */ 326 row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 327 328 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 329 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 330 331 /* 332 Only include those interpolation points that are truly 333 nonzero. Note this is very important for final grid lines 334 in x and y directions; since they have no right/top neighbors 335 */ 336 if (coors) { 337 /* only access the next coors point if we know there is one */ 338 /* note this is dangerous because x may not exactly equal ZERO */ 339 x = (coors[j][i].x - ccoors[j_c][i_c].x); 340 if (PetscAbsScalar(x) != 0.0) x = x/(ccoors[j_c][i_c+1].x - ccoors[j_c][i_c].x); 341 y = (coors[j][i].y - ccoors[j_c][i_c].y); 342 if (PetscAbsScalar(y) != 0.0) y = y/(ccoors[j_c+1][i_c].y - ccoors[j_c][i_c].y); 343 } else { 344 x = ((double)(i - i_c*ratioi))/((double)ratioi); 345 y = ((double)(j - j_c*ratioj))/((double)ratioj); 346 } 347 nc = 0; 348 /* one left and below; or we are right on it */ 349 col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 350 cols[nc] = col_shift + idx_c[col]/dof; 351 v[nc++] = x*y - x - y + 1.0; 352 /* one right and below */ 353 if (i_c*ratioi != i) { 354 cols[nc] = col_shift + idx_c[col+dof]/dof; 355 v[nc++] = -x*y + x; 356 } 357 /* one left and above */ 358 if (j_c*ratioj != j) { 359 cols[nc] = col_shift + idx_c[col+m_ghost_c*dof]/dof; 360 v[nc++] = -x*y + y; 361 } 362 /* one right and above */ 363 if (j_c*ratioj != j && i_c*ratioi != i) { 364 cols[nc] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof; 365 v[nc++] = x*y; 366 } 367 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 368 } 369 } 370 if (vcoors) { 371 ierr = DAVecRestoreArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 372 ierr = DAVecRestoreArray(ddc->da_coordinates,cvcoors,&ccoors);CHKERRQ(ierr); 373 } 374 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 375 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 376 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 377 ierr = MatDestroy(mat);CHKERRQ(ierr); 378 ierr = PetscLogFlops(13.0*m_f*n_f);CHKERRQ(ierr); 379 PetscFunctionReturn(0); 380 } 381 382 /* 383 Contributed by Andrei Draganescu <aidraga@sandia.gov> 384 */ 385 #undef __FUNCT__ 386 #define __FUNCT__ "DAGetInterpolation_2D_Q0" 387 PetscErrorCode DAGetInterpolation_2D_Q0(DM dac,DM daf,Mat *A) 388 { 389 PetscErrorCode ierr; 390 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof; 391 PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,*dnz,*onz; 392 PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj; 393 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; 394 PetscMPIInt size_c,size_f,rank_f; 395 PetscScalar v[4]; 396 Mat mat; 397 DAPeriodicType pt; 398 399 PetscFunctionBegin; 400 ierr = DAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 401 ierr = DAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 402 if (DAXPeriodic(pt)) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in x"); 403 if (DAYPeriodic(pt)) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in y"); 404 ratioi = mx/Mx; 405 ratioj = my/My; 406 if (ratioi*Mx != mx) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x"); 407 if (ratioj*My != my) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y"); 408 if (ratioi != 2) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 2"); 409 if (ratioj != 2) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 2"); 410 411 ierr = DAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 412 ierr = DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 413 ierr = DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 414 415 ierr = DAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 416 ierr = DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 417 ierr = DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 418 419 /* 420 Used for handling a coarse DA that lives on 1/4 the processors of the fine DA. 421 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 422 processors). It's effective length is hence 4 times its normal length, this is 423 why the col_scale is multiplied by the interpolation matrix column sizes. 424 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 425 copy of the coarse vector. A bit of a hack but you do better. 426 427 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 428 */ 429 ierr = MPI_Comm_size(((PetscObject)dac)->comm,&size_c);CHKERRQ(ierr); 430 ierr = MPI_Comm_size(((PetscObject)daf)->comm,&size_f);CHKERRQ(ierr); 431 ierr = MPI_Comm_rank(((PetscObject)daf)->comm,&rank_f);CHKERRQ(ierr); 432 col_scale = size_f/size_c; 433 col_shift = Mx*My*(rank_f/size_c); 434 435 ierr = MatPreallocateInitialize(((PetscObject)daf)->comm,m_f*n_f,col_scale*m_c*n_c,dnz,onz);CHKERRQ(ierr); 436 for (j=j_start; j<j_start+n_f; j++) { 437 for (i=i_start; i<i_start+m_f; i++) { 438 /* convert to local "natural" numbering and then to PETSc global numbering */ 439 row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 440 441 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 442 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 443 444 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 445 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 446 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 447 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 448 449 /* 450 Only include those interpolation points that are truly 451 nonzero. Note this is very important for final grid lines 452 in x and y directions; since they have no right/top neighbors 453 */ 454 nc = 0; 455 /* one left and below; or we are right on it */ 456 col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 457 cols[nc++] = col_shift + idx_c[col]/dof; 458 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 459 } 460 } 461 ierr = MatCreate(((PetscObject)daf)->comm,&mat);CHKERRQ(ierr); 462 ierr = MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);CHKERRQ(ierr); 463 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 464 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 465 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 466 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 467 468 /* loop over local fine grid nodes setting interpolation for those*/ 469 for (j=j_start; j<j_start+n_f; j++) { 470 for (i=i_start; i<i_start+m_f; i++) { 471 /* convert to local "natural" numbering and then to PETSc global numbering */ 472 row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 473 474 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 475 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 476 nc = 0; 477 /* one left and below; or we are right on it */ 478 col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 479 cols[nc] = col_shift + idx_c[col]/dof; 480 v[nc++] = 1.0; 481 482 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 483 } 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 ierr = PetscLogFlops(13.0*m_f*n_f);CHKERRQ(ierr); 490 PetscFunctionReturn(0); 491 } 492 493 /* 494 Contributed by Jianming Yang <jianming-yang@uiowa.edu> 495 */ 496 #undef __FUNCT__ 497 #define __FUNCT__ "DAGetInterpolation_3D_Q0" 498 PetscErrorCode DAGetInterpolation_3D_Q0(DM dac,DM daf,Mat *A) 499 { 500 PetscErrorCode ierr; 501 PetscInt i,j,l,i_start,j_start,l_start,m_f,n_f,p_f,Mx,My,Mz,*idx_f,dof; 502 PetscInt m_ghost,n_ghost,p_ghost,*idx_c,m_ghost_c,n_ghost_c,p_ghost_c,nc,*dnz,*onz; 503 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; 504 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; 505 PetscMPIInt size_c,size_f,rank_f; 506 PetscScalar v[8]; 507 Mat mat; 508 DAPeriodicType pt; 509 510 PetscFunctionBegin; 511 ierr = DAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 512 ierr = DAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 513 if (DAXPeriodic(pt)) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in x"); 514 if (DAYPeriodic(pt)) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in y"); 515 if (DAZPeriodic(pt)) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in z"); 516 ratioi = mx/Mx; 517 ratioj = my/My; 518 ratiol = mz/Mz; 519 if (ratioi*Mx != mx) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x"); 520 if (ratioj*My != my) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y"); 521 if (ratiol*Mz != mz) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in z"); 522 if (ratioi != 2 && ratioi != 1) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 1 or 2"); 523 if (ratioj != 2 && ratioj != 1) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 1 or 2"); 524 if (ratiol != 2 && ratiol != 1) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in z must be 1 or 2"); 525 526 ierr = DAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 527 ierr = DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 528 ierr = DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 529 530 ierr = DAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 531 ierr = DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr); 532 ierr = DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 533 /* 534 Used for handling a coarse DA that lives on 1/4 the processors of the fine DA. 535 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 536 processors). It's effective length is hence 4 times its normal length, this is 537 why the col_scale is multiplied by the interpolation matrix column sizes. 538 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 539 copy of the coarse vector. A bit of a hack but you do better. 540 541 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 542 */ 543 ierr = MPI_Comm_size(((PetscObject)dac)->comm,&size_c);CHKERRQ(ierr); 544 ierr = MPI_Comm_size(((PetscObject)daf)->comm,&size_f);CHKERRQ(ierr); 545 ierr = MPI_Comm_rank(((PetscObject)daf)->comm,&rank_f);CHKERRQ(ierr); 546 col_scale = size_f/size_c; 547 col_shift = Mx*My*Mz*(rank_f/size_c); 548 549 ierr = MatPreallocateInitialize(((PetscObject)daf)->comm,m_f*n_f*p_f,col_scale*m_c*n_c*p_c,dnz,onz);CHKERRQ(ierr); 550 for (l=l_start; l<l_start+p_f; l++) { 551 for (j=j_start; j<j_start+n_f; j++) { 552 for (i=i_start; i<i_start+m_f; i++) { 553 /* convert to local "natural" numbering and then to PETSc global numbering */ 554 row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 555 556 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 557 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 558 l_c = (l/ratiol); 559 560 if (l_c < l_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 561 l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c); 562 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 563 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 564 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 565 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 566 567 /* 568 Only include those interpolation points that are truly 569 nonzero. Note this is very important for final grid lines 570 in x and y directions; since they have no right/top neighbors 571 */ 572 nc = 0; 573 /* one left and below; or we are right on it */ 574 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)); 575 cols[nc++] = col_shift + idx_c[col]/dof; 576 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 577 } 578 } 579 } 580 ierr = MatCreate(((PetscObject)daf)->comm,&mat);CHKERRQ(ierr); 581 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); 582 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 583 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 584 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 585 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 586 587 /* loop over local fine grid nodes setting interpolation for those*/ 588 for (l=l_start; l<l_start+p_f; l++) { 589 for (j=j_start; j<j_start+n_f; j++) { 590 for (i=i_start; i<i_start+m_f; i++) { 591 /* convert to local "natural" numbering and then to PETSc global numbering */ 592 row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 593 594 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 595 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 596 l_c = (l/ratiol); 597 nc = 0; 598 /* one left and below; or we are right on it */ 599 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)); 600 cols[nc] = col_shift + idx_c[col]/dof; 601 v[nc++] = 1.0; 602 603 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 604 } 605 } 606 } 607 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 608 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 609 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 610 ierr = MatDestroy(mat);CHKERRQ(ierr); 611 ierr = PetscLogFlops(13.0*m_f*n_f*p_f);CHKERRQ(ierr); 612 PetscFunctionReturn(0); 613 } 614 615 #undef __FUNCT__ 616 #define __FUNCT__ "DAGetInterpolation_3D_Q1" 617 PetscErrorCode DAGetInterpolation_3D_Q1(DM dac,DM daf,Mat *A) 618 { 619 PetscErrorCode ierr; 620 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof,l; 621 PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,Mz,mz; 622 PetscInt row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok; 623 PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c; 624 PetscInt l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c; 625 PetscInt l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz; 626 PetscScalar v[8],x,y,z; 627 Mat mat; 628 DAPeriodicType pt; 629 DACoor3d ***coors = 0,***ccoors; 630 Vec vcoors,cvcoors; 631 DM_DA *ddc = (DM_DA*)dac->data, *ddf = (DM_DA*)daf->data; 632 633 PetscFunctionBegin; 634 ierr = DAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 635 ierr = DAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 636 if (mx == Mx) { 637 ratioi = 1; 638 } else if (DAXPeriodic(pt)){ 639 ratioi = mx/Mx; 640 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); 641 } else { 642 ratioi = (mx-1)/(Mx-1); 643 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); 644 } 645 if (my == My) { 646 ratioj = 1; 647 } else if (DAYPeriodic(pt)){ 648 ratioj = my/My; 649 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); 650 } else { 651 ratioj = (my-1)/(My-1); 652 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); 653 } 654 if (mz == Mz) { 655 ratiok = 1; 656 } else if (DAZPeriodic(pt)){ 657 ratiok = mz/Mz; 658 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); 659 } else { 660 ratiok = (mz-1)/(Mz-1); 661 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); 662 } 663 664 ierr = DAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 665 ierr = DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 666 ierr = DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 667 668 ierr = DAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 669 ierr = DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr); 670 ierr = DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 671 672 /* create interpolation matrix, determining exact preallocation */ 673 ierr = MatPreallocateInitialize(((PetscObject)dac)->comm,m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);CHKERRQ(ierr); 674 /* loop over local fine grid nodes counting interpolating points */ 675 for (l=l_start; l<l_start+p_f; l++) { 676 for (j=j_start; j<j_start+n_f; j++) { 677 for (i=i_start; i<i_start+m_f; i++) { 678 /* convert to local "natural" numbering and then to PETSc global numbering */ 679 row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 680 i_c = (i/ratioi); 681 j_c = (j/ratioj); 682 l_c = (l/ratiok); 683 if (l_c < l_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 684 l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c); 685 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 686 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 687 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 688 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 689 690 /* 691 Only include those interpolation points that are truly 692 nonzero. Note this is very important for final grid lines 693 in x and y directions; since they have no right/top neighbors 694 */ 695 nc = 0; 696 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)); 697 cols[nc++] = idx_c[col]/dof; 698 if (i_c*ratioi != i) { 699 cols[nc++] = idx_c[col+dof]/dof; 700 } 701 if (j_c*ratioj != j) { 702 cols[nc++] = idx_c[col+m_ghost_c*dof]/dof; 703 } 704 if (l_c*ratiok != l) { 705 cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof; 706 } 707 if (j_c*ratioj != j && i_c*ratioi != i) { 708 cols[nc++] = idx_c[col+(m_ghost_c+1)*dof]/dof; 709 } 710 if (j_c*ratioj != j && l_c*ratiok != l) { 711 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof; 712 } 713 if (i_c*ratioi != i && l_c*ratiok != l) { 714 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof; 715 } 716 if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) { 717 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof; 718 } 719 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 720 } 721 } 722 } 723 ierr = MatCreate(((PetscObject)dac)->comm,&mat);CHKERRQ(ierr); 724 ierr = MatSetSizes(mat,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz);CHKERRQ(ierr); 725 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 726 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 727 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 728 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 729 730 ierr = DAGetCoordinates(daf,&vcoors);CHKERRQ(ierr); 731 if (vcoors) { 732 ierr = DAGetGhostedCoordinates(dac,&cvcoors);CHKERRQ(ierr); 733 ierr = DAVecGetArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 734 ierr = DAVecGetArray(ddc->da_coordinates,cvcoors,&ccoors);CHKERRQ(ierr); 735 } 736 737 /* loop over local fine grid nodes setting interpolation for those*/ 738 for (l=l_start; l<l_start+p_f; l++) { 739 for (j=j_start; j<j_start+n_f; j++) { 740 for (i=i_start; i<i_start+m_f; i++) { 741 /* convert to local "natural" numbering and then to PETSc global numbering */ 742 row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 743 744 i_c = (i/ratioi); 745 j_c = (j/ratioj); 746 l_c = (l/ratiok); 747 748 /* 749 Only include those interpolation points that are truly 750 nonzero. Note this is very important for final grid lines 751 in x and y directions; since they have no right/top neighbors 752 */ 753 if (coors) { 754 /* only access the next coors point if we know there is one */ 755 /* note this is dangerous because x may not exactly equal ZERO */ 756 x = (coors[l][j][i].x - ccoors[l_c][j_c][i_c].x); 757 if (PetscAbsScalar(x) != 0.0) x = x/(ccoors[l_c][j_c][i_c+1].x - ccoors[l_c][j_c][i_c].x); 758 y = (coors[l][j][i].y - ccoors[l_c][j_c][i_c].y); 759 if (PetscAbsScalar(y) != 0.0) y = y/(ccoors[l_c][j_c+1][i_c].y - ccoors[l_c][j_c][i_c].y); 760 z = (coors[l][j][i].z - ccoors[l_c][j_c][i_c].z); 761 if (PetscAbsScalar(z) != 0.0) z = z/(ccoors[l_c+1][j_c][i_c].z - ccoors[l_c][j_c][i_c].z); 762 } else { 763 x = ((double)(i - i_c*ratioi))/((double)ratioi); 764 y = ((double)(j - j_c*ratioj))/((double)ratioj); 765 z = ((double)(l - l_c*ratiok))/((double)ratiok); 766 } 767 nc = 0; 768 /* one left and below; or we are right on it */ 769 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)); 770 771 cols[nc] = idx_c[col]/dof; 772 v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.)); 773 774 if (i_c*ratioi != i) { 775 cols[nc] = idx_c[col+dof]/dof; 776 v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.)); 777 } 778 779 if (j_c*ratioj != j) { 780 cols[nc] = idx_c[col+m_ghost_c*dof]/dof; 781 v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.)); 782 } 783 784 if (l_c*ratiok != l) { 785 cols[nc] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof; 786 v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.)); 787 } 788 789 if (j_c*ratioj != j && i_c*ratioi != i) { 790 cols[nc] = idx_c[col+(m_ghost_c+1)*dof]/dof; 791 v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.)); 792 } 793 794 if (j_c*ratioj != j && l_c*ratiok != l) { 795 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof; 796 v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.)); 797 } 798 799 if (i_c*ratioi != i && l_c*ratiok != l) { 800 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof; 801 v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.)); 802 } 803 804 if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) { 805 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof; 806 v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.)); 807 } 808 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 809 } 810 } 811 } 812 if (vcoors) { 813 ierr = DAVecRestoreArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 814 ierr = DAVecRestoreArray(ddc->da_coordinates,cvcoors,&ccoors);CHKERRQ(ierr); 815 } 816 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 817 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 818 819 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 820 ierr = MatDestroy(mat);CHKERRQ(ierr); 821 ierr = PetscLogFlops(13.0*m_f*n_f);CHKERRQ(ierr); 822 PetscFunctionReturn(0); 823 } 824 825 #undef __FUNCT__ 826 #define __FUNCT__ "DAGetInterpolation" 827 /*@C 828 DAGetInterpolation - Gets an interpolation matrix that maps between 829 grids associated with two DAs. 830 831 Collective on DA 832 833 Input Parameters: 834 + dac - the coarse grid DA 835 - daf - the fine grid DA 836 837 Output Parameters: 838 + A - the interpolation matrix 839 - scale - a scaling vector used to scale the coarse grid restricted vector before applying the 840 grid function or grid Jacobian to it. 841 842 Level: intermediate 843 844 .keywords: interpolation, restriction, multigrid 845 846 .seealso: DARefine(), DAGetInjection() 847 @*/ 848 PetscErrorCode PETSCDM_DLLEXPORT DAGetInterpolation(DM dac,DM daf,Mat *A,Vec *scale) 849 { 850 PetscErrorCode ierr; 851 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 852 DAPeriodicType wrapc,wrapf; 853 DAStencilType stc,stf; 854 DM_DA *ddc = (DM_DA*)dac->data; 855 856 PetscFunctionBegin; 857 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 858 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 859 PetscValidPointer(A,3); 860 if (scale) PetscValidPointer(scale,4); 861 862 ierr = DAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&wrapc,&stc);CHKERRQ(ierr); 863 ierr = DAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&wrapf,&stf);CHKERRQ(ierr); 864 if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 865 if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DA do not match %D %D",dofc,doff);CHKERRQ(ierr); 866 if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DA do not match %D %D",sc,sf);CHKERRQ(ierr); 867 if (wrapc != wrapf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Periodic type different in two DAs");CHKERRQ(ierr); 868 if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DAs");CHKERRQ(ierr); 869 if (Mc < 2 && Mf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); 870 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"); 871 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"); 872 873 if (ddc->interptype == DA_Q1){ 874 if (dimc == 1){ 875 ierr = DAGetInterpolation_1D_Q1(dac,daf,A);CHKERRQ(ierr); 876 } else if (dimc == 2){ 877 ierr = DAGetInterpolation_2D_Q1(dac,daf,A);CHKERRQ(ierr); 878 } else if (dimc == 3){ 879 ierr = DAGetInterpolation_3D_Q1(dac,daf,A);CHKERRQ(ierr); 880 } else { 881 SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_SUP,"No support for this DA dimension %D for interpolation type %d",dimc,(int)ddc->interptype); 882 } 883 } else if (ddc->interptype == DA_Q0){ 884 if (dimc == 1){ 885 ierr = DAGetInterpolation_1D_Q0(dac,daf,A);CHKERRQ(ierr); 886 } else if (dimc == 2){ 887 ierr = DAGetInterpolation_2D_Q0(dac,daf,A);CHKERRQ(ierr); 888 } else if (dimc == 3){ 889 ierr = DAGetInterpolation_3D_Q0(dac,daf,A);CHKERRQ(ierr); 890 } else { 891 SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_SUP,"No support for this DA dimension %D for interpolation type %d",dimc,(int)ddc->interptype); 892 } 893 } 894 if (scale) { 895 ierr = DMGetInterpolationScale((DM)dac,(DM)daf,*A,scale);CHKERRQ(ierr); 896 } 897 PetscFunctionReturn(0); 898 } 899 900 #undef __FUNCT__ 901 #define __FUNCT__ "DAGetInjection_2D" 902 PetscErrorCode DAGetInjection_2D(DM dac,DM daf,VecScatter *inject) 903 { 904 PetscErrorCode ierr; 905 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof; 906 PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c; 907 PetscInt row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj; 908 PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c; 909 PetscInt *cols; 910 DAPeriodicType pt; 911 Vec vecf,vecc; 912 IS isf; 913 914 PetscFunctionBegin; 915 916 ierr = DAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 917 ierr = DAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 918 if (DAXPeriodic(pt)){ 919 ratioi = mx/Mx; 920 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); 921 } else { 922 ratioi = (mx-1)/(Mx-1); 923 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); 924 } 925 if (DAYPeriodic(pt)){ 926 ratioj = my/My; 927 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); 928 } else { 929 ratioj = (my-1)/(My-1); 930 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); 931 } 932 933 934 ierr = DAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 935 ierr = DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 936 ierr = DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 937 938 ierr = DAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 939 ierr = DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 940 ierr = DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 941 942 943 /* loop over local fine grid nodes setting interpolation for those*/ 944 nc = 0; 945 ierr = PetscMalloc(n_f*m_f*sizeof(PetscInt),&cols);CHKERRQ(ierr); 946 for (j=j_start; j<j_start+n_f; j++) { 947 for (i=i_start; i<i_start+m_f; i++) { 948 949 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 950 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 951 952 if (j_c < j_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 953 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 954 if (i_c < i_start_ghost_c) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Processor's coarse DA must lie over fine DA\n\ 955 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 956 957 if (i_c*ratioi == i && j_c*ratioj == j) { 958 /* convert to local "natural" numbering and then to PETSc global numbering */ 959 row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]; 960 cols[nc++] = row/dof; 961 } 962 } 963 } 964 965 ierr = ISCreateBlock(((PetscObject)daf)->comm,dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 966 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 967 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 968 ierr = VecScatterCreate(vecf,isf,vecc,PETSC_NULL,inject);CHKERRQ(ierr); 969 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 970 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 971 ierr = ISDestroy(isf);CHKERRQ(ierr); 972 PetscFunctionReturn(0); 973 } 974 975 #undef __FUNCT__ 976 #define __FUNCT__ "DAGetInjection" 977 /*@ 978 DAGetInjection - Gets an injection matrix that maps between 979 grids associated with two DAs. 980 981 Collective on DA 982 983 Input Parameters: 984 + dac - the coarse grid DA 985 - daf - the fine grid DA 986 987 Output Parameters: 988 . inject - the injection scatter 989 990 Level: intermediate 991 992 .keywords: interpolation, restriction, multigrid, injection 993 994 .seealso: DARefine(), DAGetInterpolation() 995 @*/ 996 PetscErrorCode PETSCDM_DLLEXPORT DAGetInjection(DM dac,DM daf,VecScatter *inject) 997 { 998 PetscErrorCode ierr; 999 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1000 DAPeriodicType wrapc,wrapf; 1001 DAStencilType stc,stf; 1002 1003 PetscFunctionBegin; 1004 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1005 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1006 PetscValidPointer(inject,3); 1007 1008 ierr = DAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&wrapc,&stc);CHKERRQ(ierr); 1009 ierr = DAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&wrapf,&stf);CHKERRQ(ierr); 1010 if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 1011 if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DA do not match %D %D",dofc,doff);CHKERRQ(ierr); 1012 if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DA do not match %D %D",sc,sf);CHKERRQ(ierr); 1013 if (wrapc != wrapf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Periodic type different in two DAs");CHKERRQ(ierr); 1014 if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DAs");CHKERRQ(ierr); 1015 if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); 1016 if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction"); 1017 if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction"); 1018 1019 if (dimc == 2){ 1020 ierr = DAGetInjection_2D(dac,daf,inject);CHKERRQ(ierr); 1021 } else { 1022 SETERRQ1(((PetscObject)daf)->comm,PETSC_ERR_SUP,"No support for this DA dimension %D",dimc); 1023 } 1024 PetscFunctionReturn(0); 1025 } 1026 1027 #undef __FUNCT__ 1028 #define __FUNCT__ "DAGetAggregates" 1029 /*@ 1030 DAGetAggregates - Gets the aggregates that map between 1031 grids associated with two DAs. 1032 1033 Collective on DA 1034 1035 Input Parameters: 1036 + dac - the coarse grid DA 1037 - daf - the fine grid DA 1038 1039 Output Parameters: 1040 . rest - the restriction matrix (transpose of the projection matrix) 1041 1042 Level: intermediate 1043 1044 .keywords: interpolation, restriction, multigrid 1045 1046 .seealso: DARefine(), DAGetInjection(), DAGetInterpolation() 1047 @*/ 1048 PetscErrorCode PETSCDM_DLLEXPORT DAGetAggregates(DM dac,DM daf,Mat *rest) 1049 { 1050 PetscErrorCode ierr; 1051 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc; 1052 PetscInt dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1053 DAPeriodicType wrapc,wrapf; 1054 DAStencilType stc,stf; 1055 /* PetscReal r_x, r_y, r_z; */ 1056 1057 PetscInt i,j,l; 1058 PetscInt i_start,j_start,l_start, m_f,n_f,p_f; 1059 PetscInt i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost; 1060 PetscInt *idx_f; 1061 PetscInt i_c,j_c,l_c; 1062 PetscInt i_start_c,j_start_c,l_start_c, m_c,n_c,p_c; 1063 PetscInt i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c; 1064 PetscInt *idx_c; 1065 PetscInt d; 1066 PetscInt a; 1067 PetscInt max_agg_size; 1068 PetscInt *fine_nodes; 1069 PetscScalar *one_vec; 1070 PetscInt fn_idx; 1071 1072 PetscFunctionBegin; 1073 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1074 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1075 PetscValidPointer(rest,3); 1076 1077 ierr = DAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&wrapc,&stc);CHKERRQ(ierr); 1078 ierr = DAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&wrapf,&stf);CHKERRQ(ierr); 1079 if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 1080 if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DA do not match %D %D",dofc,doff);CHKERRQ(ierr); 1081 if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DA do not match %D %D",sc,sf);CHKERRQ(ierr); 1082 if (wrapc != wrapf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Periodic type different in two DAs");CHKERRQ(ierr); 1083 if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DAs");CHKERRQ(ierr); 1084 1085 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); 1086 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); 1087 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); 1088 1089 if (Pc < 0) Pc = 1; 1090 if (Pf < 0) Pf = 1; 1091 if (Nc < 0) Nc = 1; 1092 if (Nf < 0) Nf = 1; 1093 1094 ierr = DAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 1095 ierr = DAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 1096 ierr = DAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 1097 1098 ierr = DAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 1099 ierr = DAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,&l_start_ghost_c,&m_ghost_c,&n_ghost_c,&p_ghost_c);CHKERRQ(ierr); 1100 ierr = DAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 1101 1102 /* 1103 Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios 1104 for dimension 1 and 2 respectively. 1105 Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1) 1106 and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate. 1107 Each specific dof on the fine grid is mapped to one dof on the coarse grid. 1108 */ 1109 1110 max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1); 1111 1112 /* create the matrix that will contain the restriction operator */ 1113 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, 1114 max_agg_size, PETSC_NULL, max_agg_size, PETSC_NULL, rest);CHKERRQ(ierr); 1115 1116 /* store nodes in the fine grid here */ 1117 ierr = PetscMalloc2(max_agg_size,PetscScalar, &one_vec,max_agg_size,PetscInt, &fine_nodes);CHKERRQ(ierr); 1118 for(i=0; i<max_agg_size; i++) one_vec[i] = 1.0; 1119 1120 /* loop over all coarse nodes */ 1121 for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) { 1122 for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) { 1123 for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) { 1124 for(d=0; d<dofc; d++) { 1125 /* convert to local "natural" numbering and then to PETSc global numbering */ 1126 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; 1127 1128 fn_idx = 0; 1129 /* Corresponding fine points are all points (i_f, j_f, l_f) such that 1130 i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc 1131 (same for other dimensions) 1132 */ 1133 for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) { 1134 for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) { 1135 for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) { 1136 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; 1137 fn_idx++; 1138 } 1139 } 1140 } 1141 /* add all these points to one aggregate */ 1142 ierr = MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);CHKERRQ(ierr); 1143 } 1144 } 1145 } 1146 } 1147 ierr = PetscFree2(one_vec,fine_nodes);CHKERRQ(ierr); 1148 ierr = MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1149 ierr = MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1150 PetscFunctionReturn(0); 1151 } 1152