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