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. 10 */ 11 #define NEWVERSION 0 12 13 #include <private/daimpl.h> /*I "petscdmda.h" I*/ 14 #include <petscpcmg.h> 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "DMGetInterpolationScale" 18 /*@ 19 DMGetInterpolationScale - 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: DMGetInterpolation() 33 34 @*/ 35 PetscErrorCode DMGetInterpolationScale(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__ "DMGetInterpolation_DA_1D_Q1" 53 PetscErrorCode DMGetInterpolation_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,0,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__ "DMGetInterpolation_DA_1D_Q0" 179 PetscErrorCode DMGetInterpolation_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__ "DMGetInterpolation_DA_2D_Q1" 251 PetscErrorCode DMGetInterpolation_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__ "DMGetInterpolation_DA_2D_Q0" 477 PetscErrorCode DMGetInterpolation_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__ "DMGetInterpolation_DA_3D_Q0" 588 PetscErrorCode DMGetInterpolation_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__ "DMGetInterpolation_DA_3D_Q1" 707 PetscErrorCode DMGetInterpolation_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__ "DMGetInterpolation_DA" 988 PetscErrorCode DMGetInterpolation_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 = DMGetInterpolation_DA_1D_Q1(dac,daf,A);CHKERRQ(ierr); 1016 } else if (dimc == 2){ 1017 ierr = DMGetInterpolation_DA_2D_Q1(dac,daf,A);CHKERRQ(ierr); 1018 } else if (dimc == 3){ 1019 ierr = DMGetInterpolation_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 = DMGetInterpolation_DA_1D_Q0(dac,daf,A);CHKERRQ(ierr); 1024 } else if (dimc == 2){ 1025 ierr = DMGetInterpolation_DA_2D_Q0(dac,daf,A);CHKERRQ(ierr); 1026 } else if (dimc == 3){ 1027 ierr = DMGetInterpolation_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 = DMGetInterpolationScale((DM)dac,(DM)daf,*A,scale);CHKERRQ(ierr); 1032 } 1033 PetscFunctionReturn(0); 1034 } 1035 1036 #undef __FUNCT__ 1037 #define __FUNCT__ "DMGetInjection_DA_1D" 1038 PetscErrorCode DMGetInjection_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__ "DMGetInjection_DA_2D" 1097 PetscErrorCode DMGetInjection_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__ "DMGetInjection_DA_3D" 1163 PetscErrorCode DMGetInjection_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__ "DMGetInjection_DA" 1246 PetscErrorCode DMGetInjection_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 = DMGetInjection_DA_1D(dac,daf,inject);CHKERRQ(ierr); 1271 } else if (dimc == 2) { 1272 ierr = DMGetInjection_DA_2D(dac,daf,inject);CHKERRQ(ierr); 1273 } else if (dimc == 3) { 1274 ierr = DMGetInjection_DA_3D(dac,daf,inject);CHKERRQ(ierr); 1275 } 1276 PetscFunctionReturn(0); 1277 } 1278 1279 #undef __FUNCT__ 1280 #define __FUNCT__ "DMGetAggregates_DA" 1281 PetscErrorCode DMGetAggregates_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 = 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, 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