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