1 #define PETSCDM_DLL 2 3 /* 4 Code for interpolating between grids represented by DMDAs 5 */ 6 7 #include "private/daimpl.h" /*I "petscdm.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__ "DMGetInterpolation_DA_1D_Q1" 47 PetscErrorCode DMGetInterpolation_DA_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 DMDAPeriodicType pt; 57 Vec vcoors,cvcoors; 58 DM_DA *ddc = (DM_DA*)dac->data, *ddf = (DM_DA*)daf->data; 59 60 PetscFunctionBegin; 61 ierr = DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 62 ierr = DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 63 if (pt == DMDA_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 = DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);CHKERRQ(ierr); 72 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);CHKERRQ(ierr); 73 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 74 75 ierr = DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);CHKERRQ(ierr); 76 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);CHKERRQ(ierr); 77 ierr = DMDAGetGlobalIndices(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 = DMDAGetCoordinates(daf,&vcoors);CHKERRQ(ierr); 87 if (vcoors) { 88 ierr = DMDAGetGhostedCoordinates(dac,&cvcoors);CHKERRQ(ierr); 89 ierr = DMDAVecGetArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 90 ierr = DMDAVecGetArray(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 DMDA must lie over fine DMDA\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 = DMDAVecRestoreArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 128 ierr = DMDAVecRestoreArray(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__ "DMGetInterpolation_DA_1D_Q0" 140 PetscErrorCode DMGetInterpolation_DA_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 DMDAPeriodicType pt; 150 151 PetscFunctionBegin; 152 ierr = DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 153 ierr = DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 154 if (pt == DMDA_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 = DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);CHKERRQ(ierr); 163 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);CHKERRQ(ierr); 164 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 165 166 ierr = DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);CHKERRQ(ierr); 167 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);CHKERRQ(ierr); 168 ierr = DMDAGetGlobalIndices(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__ "DMGetInterpolation_DA_2D_Q1" 212 PetscErrorCode DMGetInterpolation_DA_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 DMDAPeriodicType pt; 223 DMDACoor2d **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 = DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 229 ierr = DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 230 if (DMDAXPeriodic(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 (DMDAYPeriodic(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 = DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 247 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 248 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 249 250 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 251 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 252 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 253 254 /* 255 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 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 DMDA must lie over fine DMDA\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 DMDA must lie over fine DMDA\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 = DMDAGetCoordinates(daf,&vcoors);CHKERRQ(ierr); 316 if (vcoors) { 317 ierr = DMDAGetGhostedCoordinates(dac,&cvcoors);CHKERRQ(ierr); 318 ierr = DMDAVecGetArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 319 ierr = DMDAVecGetArray(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 = DMDAVecRestoreArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 372 ierr = DMDAVecRestoreArray(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__ "DMGetInterpolation_DA_2D_Q0" 387 PetscErrorCode DMGetInterpolation_DA_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 DMDAPeriodicType pt; 398 399 PetscFunctionBegin; 400 ierr = DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 401 ierr = DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 402 if (DMDAXPeriodic(pt)) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in x"); 403 if (DMDAYPeriodic(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 = DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 412 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 413 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 414 415 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 416 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 417 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 418 419 /* 420 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 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 DMDA must lie over fine DMDA\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 DMDA must lie over fine DMDA\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__ "DMGetInterpolation_DA_3D_Q0" 498 PetscErrorCode DMGetInterpolation_DA_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 DMDAPeriodicType pt; 509 510 PetscFunctionBegin; 511 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 512 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 513 if (DMDAXPeriodic(pt)) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in x"); 514 if (DMDAYPeriodic(pt)) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in y"); 515 if (DMDAZPeriodic(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 = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 527 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 528 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 529 530 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 531 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); 532 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 533 /* 534 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 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 DMDA must lie over fine DMDA\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 DMDA must lie over fine DMDA\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 DMDA must lie over fine DMDA\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__ "DMGetInterpolation_DA_3D_Q1" 617 PetscErrorCode DMGetInterpolation_DA_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 DMDAPeriodicType pt; 629 DMDACoor3d ***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 = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 635 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 636 if (mx == Mx) { 637 ratioi = 1; 638 } else if (DMDAXPeriodic(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 (DMDAYPeriodic(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 (DMDAZPeriodic(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 = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 665 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 666 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 667 668 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 669 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); 670 ierr = DMDAGetGlobalIndices(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 DMDA must lie over fine DMDA\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 DMDA must lie over fine DMDA\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 DMDA must lie over fine DMDA\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 = DMDAGetCoordinates(daf,&vcoors);CHKERRQ(ierr); 731 if (vcoors) { 732 ierr = DMDAGetGhostedCoordinates(dac,&cvcoors);CHKERRQ(ierr); 733 ierr = DMDAVecGetArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 734 ierr = DMDAVecGetArray(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 = DMDAVecRestoreArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr); 814 ierr = DMDAVecRestoreArray(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__ "DMGetInterpolation_DA" 827 PetscErrorCode PETSCDM_DLLEXPORT DMGetInterpolation_DA(DM dac,DM daf,Mat *A,Vec *scale) 828 { 829 PetscErrorCode ierr; 830 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 831 DMDAPeriodicType wrapc,wrapf; 832 DMDAStencilType stc,stf; 833 DM_DA *ddc = (DM_DA*)dac->data; 834 835 PetscFunctionBegin; 836 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 837 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 838 PetscValidPointer(A,3); 839 if (scale) PetscValidPointer(scale,4); 840 841 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&wrapc,&stc);CHKERRQ(ierr); 842 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&wrapf,&stf);CHKERRQ(ierr); 843 if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 844 if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);CHKERRQ(ierr); 845 if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);CHKERRQ(ierr); 846 if (wrapc != wrapf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Periodic type different in two DMDAs");CHKERRQ(ierr); 847 if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");CHKERRQ(ierr); 848 if (Mc < 2 && Mf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); 849 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"); 850 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"); 851 852 if (ddc->interptype == DMDA_Q1){ 853 if (dimc == 1){ 854 ierr = DMGetInterpolation_DA_1D_Q1(dac,daf,A);CHKERRQ(ierr); 855 } else if (dimc == 2){ 856 ierr = DMGetInterpolation_DA_2D_Q1(dac,daf,A);CHKERRQ(ierr); 857 } else if (dimc == 3){ 858 ierr = DMGetInterpolation_DA_3D_Q1(dac,daf,A);CHKERRQ(ierr); 859 } else { 860 SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype); 861 } 862 } else if (ddc->interptype == DMDA_Q0){ 863 if (dimc == 1){ 864 ierr = DMGetInterpolation_DA_1D_Q0(dac,daf,A);CHKERRQ(ierr); 865 } else if (dimc == 2){ 866 ierr = DMGetInterpolation_DA_2D_Q0(dac,daf,A);CHKERRQ(ierr); 867 } else if (dimc == 3){ 868 ierr = DMGetInterpolation_DA_3D_Q0(dac,daf,A);CHKERRQ(ierr); 869 } else { 870 SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype); 871 } 872 } 873 if (scale) { 874 ierr = DMGetInterpolationScale((DM)dac,(DM)daf,*A,scale);CHKERRQ(ierr); 875 } 876 PetscFunctionReturn(0); 877 } 878 879 #undef __FUNCT__ 880 #define __FUNCT__ "DMGetInjection_DA_2D" 881 PetscErrorCode DMGetInjection_DA_2D(DM dac,DM daf,VecScatter *inject) 882 { 883 PetscErrorCode ierr; 884 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof; 885 PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c; 886 PetscInt row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj; 887 PetscInt i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c; 888 PetscInt *cols; 889 DMDAPeriodicType pt; 890 Vec vecf,vecc; 891 IS isf; 892 893 PetscFunctionBegin; 894 ierr = DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 895 ierr = DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 896 if (DMDAXPeriodic(pt)){ 897 ratioi = mx/Mx; 898 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); 899 } else { 900 ratioi = (mx-1)/(Mx-1); 901 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); 902 } 903 if (DMDAYPeriodic(pt)){ 904 ratioj = my/My; 905 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); 906 } else { 907 ratioj = (my-1)/(My-1); 908 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); 909 } 910 911 ierr = DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 912 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 913 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 914 915 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 916 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 917 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 918 919 920 /* loop over local fine grid nodes setting interpolation for those*/ 921 nc = 0; 922 ierr = PetscMalloc(n_f*m_f*sizeof(PetscInt),&cols);CHKERRQ(ierr); 923 for (j=j_start_c; j<j_start_c+n_c; j++) { 924 for (i=i_start_c; i<i_start_c+m_c; i++) { 925 PetscInt i_f = i*ratioi,j_f = j*ratioj; 926 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\ 927 j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost); 928 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\ 929 i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 930 row = idx_f[dof*(m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))]; 931 cols[nc++] = row/dof; 932 } 933 } 934 935 ierr = ISCreateBlock(((PetscObject)daf)->comm,dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 936 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 937 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 938 ierr = VecScatterCreate(vecf,isf,vecc,PETSC_NULL,inject);CHKERRQ(ierr); 939 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 940 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 941 ierr = ISDestroy(isf);CHKERRQ(ierr); 942 PetscFunctionReturn(0); 943 } 944 945 #undef __FUNCT__ 946 #define __FUNCT__ "DMGetInjection_DA_3D" 947 PetscErrorCode DMGetInjection_DA_3D(DM dac,DM daf,VecScatter *inject) 948 { 949 PetscErrorCode ierr; 950 PetscInt i,j,k,i_start,j_start,k_start,m_f,n_f,p_f,Mx,My,Mz; 951 PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c; 952 PetscInt i_start_ghost,j_start_ghost,k_start_ghost; 953 PetscInt mx,my,mz,ratioi,ratioj,ratiok; 954 PetscInt i_start_c,j_start_c,k_start_c; 955 PetscInt m_c,n_c,p_c; 956 PetscInt i_start_ghost_c,j_start_ghost_c,k_start_ghost_c; 957 PetscInt row,nc,dof; 958 PetscInt *idx_c,*idx_f; 959 PetscInt *cols; 960 DMDAPeriodicType pt; 961 Vec vecf,vecc; 962 IS isf; 963 964 PetscFunctionBegin; 965 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&pt,0);CHKERRQ(ierr); 966 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);CHKERRQ(ierr); 967 968 if (DMDAXPeriodic(pt)){ 969 ratioi = mx/Mx; 970 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); 971 } else { 972 ratioi = (mx-1)/(Mx-1); 973 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); 974 } 975 if (DMDAYPeriodic(pt)){ 976 ratioj = my/My; 977 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); 978 } else { 979 ratioj = (my-1)/(My-1); 980 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); 981 } 982 if (DMDAZPeriodic(pt)){ 983 ratiok = mz/Mz; 984 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); 985 } else { 986 ratiok = (mz-1)/(Mz-1); 987 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); 988 } 989 990 ierr = DMDAGetCorners(daf,&i_start,&j_start,&k_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 991 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&k_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 992 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 993 994 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&k_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 995 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); 996 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 997 998 999 /* loop over local fine grid nodes setting interpolation for those*/ 1000 nc = 0; 1001 ierr = PetscMalloc(n_f*m_f*p_f*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1002 for (k=k_start_c; k<k_start_c+p_c; k++) { 1003 for (j=j_start_c; j<j_start_c+n_c; j++) { 1004 for (i=i_start_c; i<i_start_c+m_c; i++) { 1005 PetscInt i_f = i*ratioi,j_f = j*ratioj,k_f = k*ratiok; 1006 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 " 1007 "k_c %D k_f %D fine ghost range [%D,%D]",k,k_f,k_start_ghost,k_start_ghost+p_ghost); 1008 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 " 1009 "j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost); 1010 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 " 1011 "i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1012 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))]; 1013 cols[nc++] = row/dof; 1014 } 1015 } 1016 } 1017 1018 ierr = ISCreateBlock(((PetscObject)daf)->comm,dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1019 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1020 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1021 ierr = VecScatterCreate(vecf,isf,vecc,PETSC_NULL,inject);CHKERRQ(ierr); 1022 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1023 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1024 ierr = ISDestroy(isf);CHKERRQ(ierr); 1025 PetscFunctionReturn(0); 1026 } 1027 1028 #undef __FUNCT__ 1029 #define __FUNCT__ "DMGetInjection_DA" 1030 PetscErrorCode PETSCDM_DLLEXPORT DMGetInjection_DA(DM dac,DM daf,VecScatter *inject) 1031 { 1032 PetscErrorCode ierr; 1033 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1034 DMDAPeriodicType wrapc,wrapf; 1035 DMDAStencilType stc,stf; 1036 1037 PetscFunctionBegin; 1038 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1039 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1040 PetscValidPointer(inject,3); 1041 1042 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&wrapc,&stc);CHKERRQ(ierr); 1043 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&wrapf,&stf);CHKERRQ(ierr); 1044 if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 1045 if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);CHKERRQ(ierr); 1046 if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);CHKERRQ(ierr); 1047 if (wrapc != wrapf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Periodic type different in two DMDAs");CHKERRQ(ierr); 1048 if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");CHKERRQ(ierr); 1049 if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); 1050 if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction"); 1051 if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction"); 1052 1053 if (dimc == 2){ 1054 ierr = DMGetInjection_DA_2D(dac,daf,inject);CHKERRQ(ierr); 1055 } else if (dimc == 3) { 1056 ierr = DMGetInjection_DA_3D(dac,daf,inject);CHKERRQ(ierr); 1057 } else { 1058 SETERRQ1(((PetscObject)daf)->comm,PETSC_ERR_SUP,"No support for this DMDA dimension %D",dimc); 1059 } 1060 PetscFunctionReturn(0); 1061 } 1062 1063 #undef __FUNCT__ 1064 #define __FUNCT__ "DMGetAggregates_DA" 1065 PetscErrorCode PETSCDM_DLLEXPORT DMGetAggregates_DA(DM dac,DM daf,Mat *rest) 1066 { 1067 PetscErrorCode ierr; 1068 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc; 1069 PetscInt dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1070 DMDAPeriodicType wrapc,wrapf; 1071 DMDAStencilType stc,stf; 1072 PetscInt i,j,l; 1073 PetscInt i_start,j_start,l_start, m_f,n_f,p_f; 1074 PetscInt i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost; 1075 PetscInt *idx_f; 1076 PetscInt i_c,j_c,l_c; 1077 PetscInt i_start_c,j_start_c,l_start_c, m_c,n_c,p_c; 1078 PetscInt i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c; 1079 PetscInt *idx_c; 1080 PetscInt d; 1081 PetscInt a; 1082 PetscInt max_agg_size; 1083 PetscInt *fine_nodes; 1084 PetscScalar *one_vec; 1085 PetscInt fn_idx; 1086 1087 PetscFunctionBegin; 1088 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1089 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1090 PetscValidPointer(rest,3); 1091 1092 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&wrapc,&stc);CHKERRQ(ierr); 1093 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&wrapf,&stf);CHKERRQ(ierr); 1094 if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 1095 if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);CHKERRQ(ierr); 1096 if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);CHKERRQ(ierr); 1097 if (wrapc != wrapf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Periodic type different in two DMDAs");CHKERRQ(ierr); 1098 if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");CHKERRQ(ierr); 1099 1100 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); 1101 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); 1102 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); 1103 1104 if (Pc < 0) Pc = 1; 1105 if (Pf < 0) Pf = 1; 1106 if (Nc < 0) Nc = 1; 1107 if (Nf < 0) Nf = 1; 1108 1109 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 1110 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 1111 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 1112 1113 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 1114 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); 1115 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 1116 1117 /* 1118 Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios 1119 for dimension 1 and 2 respectively. 1120 Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1) 1121 and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate. 1122 Each specific dof on the fine grid is mapped to one dof on the coarse grid. 1123 */ 1124 1125 max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1); 1126 1127 /* create the matrix that will contain the restriction operator */ 1128 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, 1129 max_agg_size, PETSC_NULL, max_agg_size, PETSC_NULL, rest);CHKERRQ(ierr); 1130 1131 /* store nodes in the fine grid here */ 1132 ierr = PetscMalloc2(max_agg_size,PetscScalar, &one_vec,max_agg_size,PetscInt, &fine_nodes);CHKERRQ(ierr); 1133 for(i=0; i<max_agg_size; i++) one_vec[i] = 1.0; 1134 1135 /* loop over all coarse nodes */ 1136 for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) { 1137 for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) { 1138 for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) { 1139 for(d=0; d<dofc; d++) { 1140 /* convert to local "natural" numbering and then to PETSc global numbering */ 1141 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; 1142 1143 fn_idx = 0; 1144 /* Corresponding fine points are all points (i_f, j_f, l_f) such that 1145 i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc 1146 (same for other dimensions) 1147 */ 1148 for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) { 1149 for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) { 1150 for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) { 1151 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; 1152 fn_idx++; 1153 } 1154 } 1155 } 1156 /* add all these points to one aggregate */ 1157 ierr = MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);CHKERRQ(ierr); 1158 } 1159 } 1160 } 1161 } 1162 ierr = PetscFree2(one_vec,fine_nodes);CHKERRQ(ierr); 1163 ierr = MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1164 ierr = MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1165 PetscFunctionReturn(0); 1166 } 1167