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) cols[nc++] = col_shift + idx_c[col+dof]/dof; 331 /* one left and above */ 332 if (j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+m_ghost_c*dof]/dof; 333 /* one right and above */ 334 if (i_c*ratioi != i && j_c*ratioj != j) cols[nc++] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof; 335 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 336 } 337 } 338 ierr = MatCreate(((PetscObject)daf)->comm,&mat);CHKERRQ(ierr); 339 ierr = MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);CHKERRQ(ierr); 340 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 341 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 342 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 343 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 344 345 /* loop over local fine grid nodes setting interpolation for those*/ 346 if (!NEWVERSION) { 347 348 for (j=j_start; j<j_start+n_f; j++) { 349 for (i=i_start; i<i_start+m_f; i++) { 350 /* convert to local "natural" numbering and then to PETSc global numbering */ 351 row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 352 353 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 354 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 355 356 /* 357 Only include those interpolation points that are truly 358 nonzero. Note this is very important for final grid lines 359 in x and y directions; since they have no right/top neighbors 360 */ 361 x = ((double)(i - i_c*ratioi))/((double)ratioi); 362 y = ((double)(j - j_c*ratioj))/((double)ratioj); 363 364 nc = 0; 365 /* one left and below; or we are right on it */ 366 col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 367 cols[nc] = col_shift + idx_c[col]/dof; 368 v[nc++] = x*y - x - y + 1.0; 369 /* one right and below */ 370 if (i_c*ratioi != i) { 371 cols[nc] = col_shift + idx_c[col+dof]/dof; 372 v[nc++] = -x*y + x; 373 } 374 /* one left and above */ 375 if (j_c*ratioj != j) { 376 cols[nc] = col_shift + idx_c[col+m_ghost_c*dof]/dof; 377 v[nc++] = -x*y + y; 378 } 379 /* one right and above */ 380 if (j_c*ratioj != j && i_c*ratioi != i) { 381 cols[nc] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof; 382 v[nc++] = x*y; 383 } 384 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 385 } 386 } 387 388 } else { 389 PetscScalar Ni[4]; 390 PetscScalar *xi,*eta; 391 PetscInt li,nxi,lj,neta; 392 393 /* compute local coordinate arrays */ 394 nxi = ratioi + 1; 395 neta = ratioj + 1; 396 ierr = PetscMalloc(sizeof(PetscScalar)*nxi,&xi);CHKERRQ(ierr); 397 ierr = PetscMalloc(sizeof(PetscScalar)*neta,&eta);CHKERRQ(ierr); 398 for (li=0; li<nxi; li++) { 399 xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1)); 400 } 401 for (lj=0; lj<neta; lj++) { 402 eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1)); 403 } 404 405 /* loop over local fine grid nodes setting interpolation for those*/ 406 for (j=j_start; j<j_start+n_f; j++) { 407 for (i=i_start; i<i_start+m_f; i++) { 408 /* convert to local "natural" numbering and then to PETSc global numbering */ 409 row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 410 411 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 412 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 413 414 /* remainders */ 415 li = i - ratioi * (i/ratioi); 416 if (i==mx-1) li = nxi-1; 417 lj = j - ratioj * (j/ratioj); 418 if (j==my-1) lj = neta-1; 419 420 /* corners */ 421 col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 422 cols[0] = col_shift + idx_c[col]/dof; /* left, below */ 423 Ni[0] = 1.0; 424 if ((li==0) || (li==nxi-1)) { 425 if ((lj==0) || (lj==neta-1)) { 426 ierr = MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr); 427 continue; 428 } 429 } 430 431 /* edges + interior */ 432 /* remainders */ 433 if (i==mx-1) i_c--; 434 if (j==my-1) j_c--; 435 436 col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 437 cols[0] = col_shift + idx_c[col]/dof; /* left, below */ 438 cols[1] = col_shift + idx_c[col+dof]/dof; /* right, below */ 439 cols[2] = col_shift + idx_c[col+m_ghost_c*dof]/dof; /* left, above */ 440 cols[3] = col_shift + idx_c[col+(m_ghost_c+1)*dof]/dof; /* right, above */ 441 442 Ni[0] = 0.25*(1.0-xi[li])*(1.0-eta[lj]); 443 Ni[1] = 0.25*(1.0+xi[li])*(1.0-eta[lj]); 444 Ni[2] = 0.25*(1.0-xi[li])*(1.0+eta[lj]); 445 Ni[3] = 0.25*(1.0+xi[li])*(1.0+eta[lj]); 446 447 nc = 0; 448 if (PetscAbsScalar(Ni[0])<1.0e-32) cols[0]=-1; 449 if (PetscAbsScalar(Ni[1])<1.0e-32) cols[1]=-1; 450 if (PetscAbsScalar(Ni[2])<1.0e-32) cols[2]=-1; 451 if (PetscAbsScalar(Ni[3])<1.0e-32) cols[3]=-1; 452 453 ierr = MatSetValues(mat,1,&row,4,cols,Ni,INSERT_VALUES);CHKERRQ(ierr); 454 } 455 } 456 ierr = PetscFree(xi);CHKERRQ(ierr); 457 ierr = PetscFree(eta);CHKERRQ(ierr); 458 } 459 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 460 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 461 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 462 ierr = MatDestroy(&mat);CHKERRQ(ierr); 463 PetscFunctionReturn(0); 464 } 465 466 /* 467 Contributed by Andrei Draganescu <aidraga@sandia.gov> 468 */ 469 #undef __FUNCT__ 470 #define __FUNCT__ "DMCreateInterpolation_DA_2D_Q0" 471 PetscErrorCode DMCreateInterpolation_DA_2D_Q0(DM dac,DM daf,Mat *A) 472 { 473 PetscErrorCode ierr; 474 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof; 475 PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,*dnz,*onz; 476 PetscInt row,col,i_start_ghost,j_start_ghost,cols[4],mx,m_c,my,nc,ratioi,ratioj; 477 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; 478 PetscMPIInt size_c,size_f,rank_f; 479 PetscScalar v[4]; 480 Mat mat; 481 DMDABoundaryType bx,by; 482 483 PetscFunctionBegin; 484 ierr = DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);CHKERRQ(ierr); 485 ierr = DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 486 if (bx == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in x"); 487 if (by == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in y"); 488 ratioi = mx/Mx; 489 ratioj = my/My; 490 if (ratioi*Mx != mx) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x"); 491 if (ratioj*My != my) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y"); 492 if (ratioi != 2) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 2"); 493 if (ratioj != 2) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 2"); 494 495 ierr = DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 496 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 497 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 498 499 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 500 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 501 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 502 503 /* 504 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 505 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 506 processors). It's effective length is hence 4 times its normal length, this is 507 why the col_scale is multiplied by the interpolation matrix column sizes. 508 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 509 copy of the coarse vector. A bit of a hack but you do better. 510 511 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 512 */ 513 ierr = MPI_Comm_size(((PetscObject)dac)->comm,&size_c);CHKERRQ(ierr); 514 ierr = MPI_Comm_size(((PetscObject)daf)->comm,&size_f);CHKERRQ(ierr); 515 ierr = MPI_Comm_rank(((PetscObject)daf)->comm,&rank_f);CHKERRQ(ierr); 516 col_scale = size_f/size_c; 517 col_shift = Mx*My*(rank_f/size_c); 518 519 ierr = MatPreallocateInitialize(((PetscObject)daf)->comm,m_f*n_f,col_scale*m_c*n_c,dnz,onz);CHKERRQ(ierr); 520 for (j=j_start; j<j_start+n_f; j++) { 521 for (i=i_start; i<i_start+m_f; i++) { 522 /* convert to local "natural" numbering and then to PETSc global numbering */ 523 row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 524 525 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 526 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 527 528 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\ 529 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 530 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\ 531 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 532 533 /* 534 Only include those interpolation points that are truly 535 nonzero. Note this is very important for final grid lines 536 in x and y directions; since they have no right/top neighbors 537 */ 538 nc = 0; 539 /* one left and below; or we are right on it */ 540 col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 541 cols[nc++] = col_shift + idx_c[col]/dof; 542 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 543 } 544 } 545 ierr = MatCreate(((PetscObject)daf)->comm,&mat);CHKERRQ(ierr); 546 ierr = MatSetSizes(mat,m_f*n_f,col_scale*m_c*n_c,mx*my,col_scale*Mx*My);CHKERRQ(ierr); 547 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 548 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 549 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 550 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 551 552 /* loop over local fine grid nodes setting interpolation for those*/ 553 for (j=j_start; j<j_start+n_f; j++) { 554 for (i=i_start; i<i_start+m_f; i++) { 555 /* convert to local "natural" numbering and then to PETSc global numbering */ 556 row = idx_f[dof*(m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 557 558 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 559 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 560 nc = 0; 561 /* one left and below; or we are right on it */ 562 col = dof*(m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c)); 563 cols[nc] = col_shift + idx_c[col]/dof; 564 v[nc++] = 1.0; 565 566 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 567 } 568 } 569 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 570 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 571 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 572 ierr = MatDestroy(&mat);CHKERRQ(ierr); 573 ierr = PetscLogFlops(13.0*m_f*n_f);CHKERRQ(ierr); 574 PetscFunctionReturn(0); 575 } 576 577 /* 578 Contributed by Jianming Yang <jianming-yang@uiowa.edu> 579 */ 580 #undef __FUNCT__ 581 #define __FUNCT__ "DMCreateInterpolation_DA_3D_Q0" 582 PetscErrorCode DMCreateInterpolation_DA_3D_Q0(DM dac,DM daf,Mat *A) 583 { 584 PetscErrorCode ierr; 585 PetscInt i,j,l,i_start,j_start,l_start,m_f,n_f,p_f,Mx,My,Mz,*idx_f,dof; 586 PetscInt m_ghost,n_ghost,p_ghost,*idx_c,m_ghost_c,n_ghost_c,p_ghost_c,nc,*dnz,*onz; 587 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; 588 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; 589 PetscMPIInt size_c,size_f,rank_f; 590 PetscScalar v[8]; 591 Mat mat; 592 DMDABoundaryType bx,by,bz; 593 594 PetscFunctionBegin; 595 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 596 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 597 if (bx == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in x"); 598 if (by == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in y"); 599 if (bz == DMDA_BOUNDARY_PERIODIC) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Cannot handle periodic grid in z"); 600 ratioi = mx/Mx; 601 ratioj = my/My; 602 ratiol = mz/Mz; 603 if (ratioi*Mx != mx) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in x"); 604 if (ratioj*My != my) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in y"); 605 if (ratiol*Mz != mz) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Fine grid points must be multiple of coarse grid points in z"); 606 if (ratioi != 2 && ratioi != 1) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in x must be 1 or 2"); 607 if (ratioj != 2 && ratioj != 1) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in y must be 1 or 2"); 608 if (ratiol != 2 && ratiol != 1) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_WRONG,"Coarsening factor in z must be 1 or 2"); 609 610 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 611 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 612 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 613 614 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 615 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); 616 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 617 /* 618 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 619 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 620 processors). It's effective length is hence 4 times its normal length, this is 621 why the col_scale is multiplied by the interpolation matrix column sizes. 622 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 623 copy of the coarse vector. A bit of a hack but you do better. 624 625 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 626 */ 627 ierr = MPI_Comm_size(((PetscObject)dac)->comm,&size_c);CHKERRQ(ierr); 628 ierr = MPI_Comm_size(((PetscObject)daf)->comm,&size_f);CHKERRQ(ierr); 629 ierr = MPI_Comm_rank(((PetscObject)daf)->comm,&rank_f);CHKERRQ(ierr); 630 col_scale = size_f/size_c; 631 col_shift = Mx*My*Mz*(rank_f/size_c); 632 633 ierr = MatPreallocateInitialize(((PetscObject)daf)->comm,m_f*n_f*p_f,col_scale*m_c*n_c*p_c,dnz,onz);CHKERRQ(ierr); 634 for (l=l_start; l<l_start+p_f; l++) { 635 for (j=j_start; j<j_start+n_f; j++) { 636 for (i=i_start; i<i_start+m_f; i++) { 637 /* convert to local "natural" numbering and then to PETSc global numbering */ 638 row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 639 640 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 641 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 642 l_c = (l/ratiol); 643 644 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\ 645 l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c); 646 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\ 647 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 648 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\ 649 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 650 651 /* 652 Only include those interpolation points that are truly 653 nonzero. Note this is very important for final grid lines 654 in x and y directions; since they have no right/top neighbors 655 */ 656 nc = 0; 657 /* one left and below; or we are right on it */ 658 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)); 659 cols[nc++] = col_shift + idx_c[col]/dof; 660 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 661 } 662 } 663 } 664 ierr = MatCreate(((PetscObject)daf)->comm,&mat);CHKERRQ(ierr); 665 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); 666 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 667 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 668 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 669 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 670 671 /* loop over local fine grid nodes setting interpolation for those*/ 672 for (l=l_start; l<l_start+p_f; l++) { 673 for (j=j_start; j<j_start+n_f; j++) { 674 for (i=i_start; i<i_start+m_f; i++) { 675 /* convert to local "natural" numbering and then to PETSc global numbering */ 676 row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 677 678 i_c = (i/ratioi); /* coarse grid node to left of fine grid node */ 679 j_c = (j/ratioj); /* coarse grid node below fine grid node */ 680 l_c = (l/ratiol); 681 nc = 0; 682 /* one left and below; or we are right on it */ 683 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)); 684 cols[nc] = col_shift + idx_c[col]/dof; 685 v[nc++] = 1.0; 686 687 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 688 } 689 } 690 } 691 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 692 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 693 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 694 ierr = MatDestroy(&mat);CHKERRQ(ierr); 695 ierr = PetscLogFlops(13.0*m_f*n_f*p_f);CHKERRQ(ierr); 696 PetscFunctionReturn(0); 697 } 698 699 #undef __FUNCT__ 700 #define __FUNCT__ "DMCreateInterpolation_DA_3D_Q1" 701 PetscErrorCode DMCreateInterpolation_DA_3D_Q1(DM dac,DM daf,Mat *A) 702 { 703 PetscErrorCode ierr; 704 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof,l; 705 PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,Mz,mz; 706 PetscInt row,col,i_start_ghost,j_start_ghost,cols[8],mx,m_c,my,nc,ratioi,ratioj,ratiok; 707 PetscInt i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c; 708 PetscInt l_start,p_f,l_start_ghost,p_ghost,l_start_c,p_c; 709 PetscInt l_start_ghost_c,p_ghost_c,l_c,*dnz,*onz; 710 PetscScalar v[8],x,y,z; 711 Mat mat; 712 DMDABoundaryType bx,by,bz; 713 714 PetscFunctionBegin; 715 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 716 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 717 if (mx == Mx) { 718 ratioi = 1; 719 } else if (bx == DMDA_BOUNDARY_PERIODIC) { 720 ratioi = mx/Mx; 721 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); 722 } else { 723 ratioi = (mx-1)/(Mx-1); 724 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); 725 } 726 if (my == My) { 727 ratioj = 1; 728 } else if (by == DMDA_BOUNDARY_PERIODIC) { 729 ratioj = my/My; 730 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); 731 } else { 732 ratioj = (my-1)/(My-1); 733 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); 734 } 735 if (mz == Mz) { 736 ratiok = 1; 737 } else if (bz == DMDA_BOUNDARY_PERIODIC) { 738 ratiok = mz/Mz; 739 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); 740 } else { 741 ratiok = (mz-1)/(Mz-1); 742 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); 743 } 744 745 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 746 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 747 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 748 749 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 750 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); 751 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 752 753 /* create interpolation matrix, determining exact preallocation */ 754 ierr = MatPreallocateInitialize(((PetscObject)dac)->comm,m_f*n_f*p_f,m_c*n_c*p_c,dnz,onz);CHKERRQ(ierr); 755 /* loop over local fine grid nodes counting interpolating points */ 756 for (l=l_start; l<l_start+p_f; l++) { 757 for (j=j_start; j<j_start+n_f; j++) { 758 for (i=i_start; i<i_start+m_f; i++) { 759 /* convert to local "natural" numbering and then to PETSc global numbering */ 760 row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 761 i_c = (i/ratioi); 762 j_c = (j/ratioj); 763 l_c = (l/ratiok); 764 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\ 765 l_start %D l_c %D l_start_ghost_c %D",l_start,l_c,l_start_ghost_c); 766 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\ 767 j_start %D j_c %D j_start_ghost_c %D",j_start,j_c,j_start_ghost_c); 768 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\ 769 i_start %D i_c %D i_start_ghost_c %D",i_start,i_c,i_start_ghost_c); 770 771 /* 772 Only include those interpolation points that are truly 773 nonzero. Note this is very important for final grid lines 774 in x and y directions; since they have no right/top neighbors 775 */ 776 nc = 0; 777 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)); 778 cols[nc++] = idx_c[col]/dof; 779 if (i_c*ratioi != i) { 780 cols[nc++] = idx_c[col+dof]/dof; 781 } 782 if (j_c*ratioj != j) { 783 cols[nc++] = idx_c[col+m_ghost_c*dof]/dof; 784 } 785 if (l_c*ratiok != l) { 786 cols[nc++] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof; 787 } 788 if (j_c*ratioj != j && i_c*ratioi != i) { 789 cols[nc++] = idx_c[col+(m_ghost_c+1)*dof]/dof; 790 } 791 if (j_c*ratioj != j && l_c*ratiok != l) { 792 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof; 793 } 794 if (i_c*ratioi != i && l_c*ratiok != l) { 795 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof; 796 } 797 if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) { 798 cols[nc++] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof; 799 } 800 ierr = MatPreallocateSet(row,nc,cols,dnz,onz);CHKERRQ(ierr); 801 } 802 } 803 } 804 ierr = MatCreate(((PetscObject)dac)->comm,&mat);CHKERRQ(ierr); 805 ierr = MatSetSizes(mat,m_f*n_f*p_f,m_c*n_c*p_c,mx*my*mz,Mx*My*Mz);CHKERRQ(ierr); 806 ierr = MatSetType(mat,MATAIJ);CHKERRQ(ierr); 807 ierr = MatSeqAIJSetPreallocation(mat,0,dnz);CHKERRQ(ierr); 808 ierr = MatMPIAIJSetPreallocation(mat,0,dnz,0,onz);CHKERRQ(ierr); 809 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 810 811 /* loop over local fine grid nodes setting interpolation for those*/ 812 if (!NEWVERSION) { 813 814 for (l=l_start; l<l_start+p_f; l++) { 815 for (j=j_start; j<j_start+n_f; j++) { 816 for (i=i_start; i<i_start+m_f; i++) { 817 /* convert to local "natural" numbering and then to PETSc global numbering */ 818 row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 819 820 i_c = (i/ratioi); 821 j_c = (j/ratioj); 822 l_c = (l/ratiok); 823 824 /* 825 Only include those interpolation points that are truly 826 nonzero. Note this is very important for final grid lines 827 in x and y directions; since they have no right/top neighbors 828 */ 829 x = ((double)(i - i_c*ratioi))/((double)ratioi); 830 y = ((double)(j - j_c*ratioj))/((double)ratioj); 831 z = ((double)(l - l_c*ratiok))/((double)ratiok); 832 833 nc = 0; 834 /* one left and below; or we are right on it */ 835 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)); 836 837 cols[nc] = idx_c[col]/dof; 838 v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.)); 839 840 if (i_c*ratioi != i) { 841 cols[nc] = idx_c[col+dof]/dof; 842 v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.)); 843 } 844 845 if (j_c*ratioj != j) { 846 cols[nc] = idx_c[col+m_ghost_c*dof]/dof; 847 v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.)); 848 } 849 850 if (l_c*ratiok != l) { 851 cols[nc] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof; 852 v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.)); 853 } 854 855 if (j_c*ratioj != j && i_c*ratioi != i) { 856 cols[nc] = idx_c[col+(m_ghost_c+1)*dof]/dof; 857 v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.)); 858 } 859 860 if (j_c*ratioj != j && l_c*ratiok != l) { 861 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof; 862 v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.)); 863 } 864 865 if (i_c*ratioi != i && l_c*ratiok != l) { 866 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof; 867 v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.)); 868 } 869 870 if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) { 871 cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof; 872 v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.)); 873 } 874 ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr); 875 } 876 } 877 } 878 879 } else { 880 PetscScalar *xi,*eta,*zeta; 881 PetscInt li,nxi,lj,neta,lk,nzeta,n; 882 PetscScalar Ni[8]; 883 884 /* compute local coordinate arrays */ 885 nxi = ratioi + 1; 886 neta = ratioj + 1; 887 nzeta = ratiok + 1; 888 ierr = PetscMalloc(sizeof(PetscScalar)*nxi,&xi);CHKERRQ(ierr); 889 ierr = PetscMalloc(sizeof(PetscScalar)*neta,&eta);CHKERRQ(ierr); 890 ierr = PetscMalloc(sizeof(PetscScalar)*nzeta,&zeta);CHKERRQ(ierr); 891 for (li=0; li<nxi; li++) xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1)); 892 for (lj=0; lj<neta; lj++) eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1)); 893 for (lk=0; lk<nzeta; lk++) zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1)); 894 895 for (l=l_start; l<l_start+p_f; l++) { 896 for (j=j_start; j<j_start+n_f; j++) { 897 for (i=i_start; i<i_start+m_f; i++) { 898 /* convert to local "natural" numbering and then to PETSc global numbering */ 899 row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof; 900 901 i_c = (i/ratioi); 902 j_c = (j/ratioj); 903 l_c = (l/ratiok); 904 905 /* remainders */ 906 li = i - ratioi * (i/ratioi); 907 if (i==mx-1) li = nxi-1; 908 lj = j - ratioj * (j/ratioj); 909 if (j==my-1) lj = neta-1; 910 lk = l - ratiok * (l/ratiok); 911 if (l==mz-1) lk = nzeta-1; 912 913 /* corners */ 914 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)); 915 cols[0] = idx_c[col]/dof; 916 Ni[0] = 1.0; 917 if ((li==0) || (li==nxi-1)) { 918 if ((lj==0) || (lj==neta-1)) { 919 if ((lk==0) || (lk==nzeta-1)) { 920 ierr = MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr); 921 continue; 922 } 923 } 924 } 925 926 /* edges + interior */ 927 /* remainders */ 928 if (i==mx-1) i_c--; 929 if (j==my-1) j_c--; 930 if (l==mz-1) l_c--; 931 932 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)); 933 cols[0] = idx_c[col]/dof; /* one left and below; or we are right on it */ 934 cols[1] = idx_c[col+dof]/dof; /* one right and below */ 935 cols[2] = idx_c[col+m_ghost_c*dof]/dof; /* one left and above */ 936 cols[3] = idx_c[col+(m_ghost_c+1)*dof]/dof; /* one right and above */ 937 938 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 */ 939 cols[5] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof; /* one right and below, and front */ 940 cols[6] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof; /* one left and above and front*/ 941 cols[7] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof; /* one right and above and front */ 942 943 Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]); 944 Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]); 945 Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]); 946 Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]); 947 948 Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]); 949 Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]); 950 Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]); 951 Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]); 952 953 for (n=0; n<8; n++) { 954 if (PetscAbsScalar(Ni[n])<1.0e-32) cols[n]=-1; 955 } 956 ierr = MatSetValues(mat,1,&row,8,cols,Ni,INSERT_VALUES);CHKERRQ(ierr); 957 958 } 959 } 960 } 961 ierr = PetscFree(xi);CHKERRQ(ierr); 962 ierr = PetscFree(eta);CHKERRQ(ierr); 963 ierr = PetscFree(zeta);CHKERRQ(ierr); 964 } 965 966 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 967 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 968 969 ierr = MatCreateMAIJ(mat,dof,A);CHKERRQ(ierr); 970 ierr = MatDestroy(&mat);CHKERRQ(ierr); 971 PetscFunctionReturn(0); 972 } 973 974 #undef __FUNCT__ 975 #define __FUNCT__ "DMCreateInterpolation_DA" 976 PetscErrorCode DMCreateInterpolation_DA(DM dac,DM daf,Mat *A,Vec *scale) 977 { 978 PetscErrorCode ierr; 979 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 980 DMDABoundaryType bxc,byc,bzc,bxf,byf,bzf; 981 DMDAStencilType stc,stf; 982 DM_DA *ddc = (DM_DA*)dac->data; 983 984 PetscFunctionBegin; 985 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 986 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 987 PetscValidPointer(A,3); 988 if (scale) PetscValidPointer(scale,4); 989 990 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr); 991 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr); 992 if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 993 if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);CHKERRQ(ierr); 994 if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);CHKERRQ(ierr); 995 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");CHKERRQ(ierr); 996 if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");CHKERRQ(ierr); 997 if (Mc < 2 && Mf > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); 998 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"); 999 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"); 1000 1001 if (ddc->interptype == DMDA_Q1) { 1002 if (dimc == 1) { 1003 ierr = DMCreateInterpolation_DA_1D_Q1(dac,daf,A);CHKERRQ(ierr); 1004 } else if (dimc == 2) { 1005 ierr = DMCreateInterpolation_DA_2D_Q1(dac,daf,A);CHKERRQ(ierr); 1006 } else if (dimc == 3) { 1007 ierr = DMCreateInterpolation_DA_3D_Q1(dac,daf,A);CHKERRQ(ierr); 1008 } else SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype); 1009 } else if (ddc->interptype == DMDA_Q0) { 1010 if (dimc == 1) { 1011 ierr = DMCreateInterpolation_DA_1D_Q0(dac,daf,A);CHKERRQ(ierr); 1012 } else if (dimc == 2) { 1013 ierr = DMCreateInterpolation_DA_2D_Q0(dac,daf,A);CHKERRQ(ierr); 1014 } else if (dimc == 3) { 1015 ierr = DMCreateInterpolation_DA_3D_Q0(dac,daf,A);CHKERRQ(ierr); 1016 } else SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_SUP,"No support for this DMDA dimension %D for interpolation type %d",dimc,(int)ddc->interptype); 1017 } 1018 if (scale) { 1019 ierr = DMCreateInterpolationScale((DM)dac,(DM)daf,*A,scale);CHKERRQ(ierr); 1020 } 1021 PetscFunctionReturn(0); 1022 } 1023 1024 #undef __FUNCT__ 1025 #define __FUNCT__ "DMCreateInjection_DA_1D" 1026 PetscErrorCode DMCreateInjection_DA_1D(DM dac,DM daf,VecScatter *inject) 1027 { 1028 PetscErrorCode ierr; 1029 PetscInt i,i_start,m_f,Mx,*idx_f,dof; 1030 PetscInt m_ghost,*idx_c,m_ghost_c; 1031 PetscInt row,i_start_ghost,mx,m_c,nc,ratioi; 1032 PetscInt i_start_c,i_start_ghost_c; 1033 PetscInt *cols; 1034 DMDABoundaryType bx; 1035 Vec vecf,vecc; 1036 IS isf; 1037 1038 PetscFunctionBegin; 1039 ierr = DMDAGetInfo(dac,0,&Mx,0,0,0,0,0,0,0,&bx,0,0,0);CHKERRQ(ierr); 1040 ierr = DMDAGetInfo(daf,0,&mx,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 1041 if (bx == DMDA_BOUNDARY_PERIODIC) { 1042 ratioi = mx/Mx; 1043 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); 1044 } else { 1045 ratioi = (mx-1)/(Mx-1); 1046 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); 1047 } 1048 1049 ierr = DMDAGetCorners(daf,&i_start,0,0,&m_f,0,0);CHKERRQ(ierr); 1050 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,0,0,&m_ghost,0,0);CHKERRQ(ierr); 1051 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 1052 1053 ierr = DMDAGetCorners(dac,&i_start_c,0,0,&m_c,0,0);CHKERRQ(ierr); 1054 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,0,0,&m_ghost_c,0,0);CHKERRQ(ierr); 1055 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 1056 1057 1058 /* loop over local fine grid nodes setting interpolation for those*/ 1059 nc = 0; 1060 ierr = PetscMalloc(m_f*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1061 1062 1063 for (i=i_start_c; i<i_start_c+m_c; i++) { 1064 PetscInt i_f = i*ratioi; 1065 1066 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\ni_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1067 1068 row = idx_f[dof*(i_f-i_start_ghost)]; 1069 cols[nc++] = row/dof; 1070 } 1071 1072 1073 ierr = ISCreateBlock(((PetscObject)daf)->comm,dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1074 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1075 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1076 ierr = VecScatterCreate(vecf,isf,vecc,PETSC_NULL,inject);CHKERRQ(ierr); 1077 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1078 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1079 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1080 PetscFunctionReturn(0); 1081 } 1082 1083 #undef __FUNCT__ 1084 #define __FUNCT__ "DMCreateInjection_DA_2D" 1085 PetscErrorCode DMCreateInjection_DA_2D(DM dac,DM daf,VecScatter *inject) 1086 { 1087 PetscErrorCode ierr; 1088 PetscInt i,j,i_start,j_start,m_f,n_f,Mx,My,*idx_f,dof; 1089 PetscInt m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c; 1090 PetscInt row,i_start_ghost,j_start_ghost,mx,m_c,my,nc,ratioi,ratioj; 1091 PetscInt i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c; 1092 PetscInt *cols; 1093 DMDABoundaryType bx,by; 1094 Vec vecf,vecc; 1095 IS isf; 1096 1097 PetscFunctionBegin; 1098 ierr = DMDAGetInfo(dac,0,&Mx,&My,0,0,0,0,0,0,&bx,&by,0,0);CHKERRQ(ierr); 1099 ierr = DMDAGetInfo(daf,0,&mx,&my,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 1100 if (bx == DMDA_BOUNDARY_PERIODIC) { 1101 ratioi = mx/Mx; 1102 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); 1103 } else { 1104 ratioi = (mx-1)/(Mx-1); 1105 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); 1106 } 1107 if (by == DMDA_BOUNDARY_PERIODIC) { 1108 ratioj = my/My; 1109 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); 1110 } else { 1111 ratioj = (my-1)/(My-1); 1112 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); 1113 } 1114 1115 ierr = DMDAGetCorners(daf,&i_start,&j_start,0,&m_f,&n_f,0);CHKERRQ(ierr); 1116 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);CHKERRQ(ierr); 1117 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 1118 1119 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,0,&m_c,&n_c,0);CHKERRQ(ierr); 1120 ierr = DMDAGetGhostCorners(dac,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);CHKERRQ(ierr); 1121 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 1122 1123 1124 /* loop over local fine grid nodes setting interpolation for those*/ 1125 nc = 0; 1126 ierr = PetscMalloc(n_f*m_f*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1127 for (j=j_start_c; j<j_start_c+n_c; j++) { 1128 for (i=i_start_c; i<i_start_c+m_c; i++) { 1129 PetscInt i_f = i*ratioi,j_f = j*ratioj; 1130 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\ 1131 j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost); 1132 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\ 1133 i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1134 row = idx_f[dof*(m_ghost*(j_f-j_start_ghost) + (i_f-i_start_ghost))]; 1135 cols[nc++] = row/dof; 1136 } 1137 } 1138 1139 ierr = ISCreateBlock(((PetscObject)daf)->comm,dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1140 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1141 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1142 ierr = VecScatterCreate(vecf,isf,vecc,PETSC_NULL,inject);CHKERRQ(ierr); 1143 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1144 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1145 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1146 PetscFunctionReturn(0); 1147 } 1148 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "DMCreateInjection_DA_3D" 1151 PetscErrorCode DMCreateInjection_DA_3D(DM dac,DM daf,VecScatter *inject) 1152 { 1153 PetscErrorCode ierr; 1154 PetscInt i,j,k,i_start,j_start,k_start,m_f,n_f,p_f,Mx,My,Mz; 1155 PetscInt m_ghost,n_ghost,p_ghost,m_ghost_c,n_ghost_c,p_ghost_c; 1156 PetscInt i_start_ghost,j_start_ghost,k_start_ghost; 1157 PetscInt mx,my,mz,ratioi,ratioj,ratiok; 1158 PetscInt i_start_c,j_start_c,k_start_c; 1159 PetscInt m_c,n_c,p_c; 1160 PetscInt i_start_ghost_c,j_start_ghost_c,k_start_ghost_c; 1161 PetscInt row,nc,dof; 1162 PetscInt *idx_c,*idx_f; 1163 PetscInt *cols; 1164 DMDABoundaryType bx,by,bz; 1165 Vec vecf,vecc; 1166 IS isf; 1167 1168 PetscFunctionBegin; 1169 ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&bx,&by,&bz,0);CHKERRQ(ierr); 1170 ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 1171 1172 if (bx == DMDA_BOUNDARY_PERIODIC) { 1173 ratioi = mx/Mx; 1174 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); 1175 } else { 1176 ratioi = (mx-1)/(Mx-1); 1177 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); 1178 } 1179 if (by == DMDA_BOUNDARY_PERIODIC) { 1180 ratioj = my/My; 1181 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); 1182 } else { 1183 ratioj = (my-1)/(My-1); 1184 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); 1185 } 1186 if (bz == DMDA_BOUNDARY_PERIODIC) { 1187 ratiok = mz/Mz; 1188 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); 1189 } else { 1190 ratiok = (mz-1)/(Mz-1); 1191 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); 1192 } 1193 1194 ierr = DMDAGetCorners(daf,&i_start,&j_start,&k_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 1195 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&k_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 1196 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 1197 1198 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&k_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 1199 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); 1200 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 1201 1202 1203 /* loop over local fine grid nodes setting interpolation for those*/ 1204 nc = 0; 1205 ierr = PetscMalloc(n_f*m_f*p_f*sizeof(PetscInt),&cols);CHKERRQ(ierr); 1206 for (k=k_start_c; k<k_start_c+p_c; k++) { 1207 for (j=j_start_c; j<j_start_c+n_c; j++) { 1208 for (i=i_start_c; i<i_start_c+m_c; i++) { 1209 PetscInt i_f = i*ratioi,j_f = j*ratioj,k_f = k*ratiok; 1210 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 " 1211 "k_c %D k_f %D fine ghost range [%D,%D]",k,k_f,k_start_ghost,k_start_ghost+p_ghost); 1212 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 " 1213 "j_c %D j_f %D fine ghost range [%D,%D]",j,j_f,j_start_ghost,j_start_ghost+n_ghost); 1214 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 " 1215 "i_c %D i_f %D fine ghost range [%D,%D]",i,i_f,i_start_ghost,i_start_ghost+m_ghost); 1216 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))]; 1217 cols[nc++] = row/dof; 1218 } 1219 } 1220 } 1221 1222 ierr = ISCreateBlock(((PetscObject)daf)->comm,dof,nc,cols,PETSC_OWN_POINTER,&isf);CHKERRQ(ierr); 1223 ierr = DMGetGlobalVector(dac,&vecc);CHKERRQ(ierr); 1224 ierr = DMGetGlobalVector(daf,&vecf);CHKERRQ(ierr); 1225 ierr = VecScatterCreate(vecf,isf,vecc,PETSC_NULL,inject);CHKERRQ(ierr); 1226 ierr = DMRestoreGlobalVector(dac,&vecc);CHKERRQ(ierr); 1227 ierr = DMRestoreGlobalVector(daf,&vecf);CHKERRQ(ierr); 1228 ierr = ISDestroy(&isf);CHKERRQ(ierr); 1229 PetscFunctionReturn(0); 1230 } 1231 1232 #undef __FUNCT__ 1233 #define __FUNCT__ "DMCreateInjection_DA" 1234 PetscErrorCode DMCreateInjection_DA(DM dac,DM daf,VecScatter *inject) 1235 { 1236 PetscErrorCode ierr; 1237 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc,dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1238 DMDABoundaryType bxc,byc,bzc,bxf,byf,bzf; 1239 DMDAStencilType stc,stf; 1240 1241 PetscFunctionBegin; 1242 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1243 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1244 PetscValidPointer(inject,3); 1245 1246 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr); 1247 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr); 1248 if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 1249 if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);CHKERRQ(ierr); 1250 if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);CHKERRQ(ierr); 1251 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");CHKERRQ(ierr); 1252 if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");CHKERRQ(ierr); 1253 if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); 1254 if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction"); 1255 if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction"); 1256 1257 if (dimc == 1) { 1258 ierr = DMCreateInjection_DA_1D(dac,daf,inject);CHKERRQ(ierr); 1259 } else if (dimc == 2) { 1260 ierr = DMCreateInjection_DA_2D(dac,daf,inject);CHKERRQ(ierr); 1261 } else if (dimc == 3) { 1262 ierr = DMCreateInjection_DA_3D(dac,daf,inject);CHKERRQ(ierr); 1263 } 1264 PetscFunctionReturn(0); 1265 } 1266 1267 #undef __FUNCT__ 1268 #define __FUNCT__ "DMCreateAggregates_DA" 1269 PetscErrorCode DMCreateAggregates_DA(DM dac,DM daf,Mat *rest) 1270 { 1271 PetscErrorCode ierr; 1272 PetscInt dimc,Mc,Nc,Pc,mc,nc,pc,dofc,sc; 1273 PetscInt dimf,Mf,Nf,Pf,mf,nf,pf,doff,sf; 1274 DMDABoundaryType bxc,byc,bzc,bxf,byf,bzf; 1275 DMDAStencilType stc,stf; 1276 PetscInt i,j,l; 1277 PetscInt i_start,j_start,l_start, m_f,n_f,p_f; 1278 PetscInt i_start_ghost,j_start_ghost,l_start_ghost,m_ghost,n_ghost,p_ghost; 1279 PetscInt *idx_f; 1280 PetscInt i_c,j_c,l_c; 1281 PetscInt i_start_c,j_start_c,l_start_c, m_c,n_c,p_c; 1282 PetscInt i_start_ghost_c,j_start_ghost_c,l_start_ghost_c,m_ghost_c,n_ghost_c,p_ghost_c; 1283 PetscInt *idx_c; 1284 PetscInt d; 1285 PetscInt a; 1286 PetscInt max_agg_size; 1287 PetscInt *fine_nodes; 1288 PetscScalar *one_vec; 1289 PetscInt fn_idx; 1290 1291 PetscFunctionBegin; 1292 PetscValidHeaderSpecific(dac,DM_CLASSID,1); 1293 PetscValidHeaderSpecific(daf,DM_CLASSID,2); 1294 PetscValidPointer(rest,3); 1295 1296 ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&bxc,&byc,&bzc,&stc);CHKERRQ(ierr); 1297 ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&bxf,&byf,&bzf,&stf);CHKERRQ(ierr); 1298 if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);CHKERRQ(ierr); 1299 if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);CHKERRQ(ierr); 1300 if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);CHKERRQ(ierr); 1301 if (bxc != bxf || byc != byf || bzc != bzf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Boundary type different in two DMDAs");CHKERRQ(ierr); 1302 if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");CHKERRQ(ierr); 1303 1304 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); 1305 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); 1306 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); 1307 1308 if (Pc < 0) Pc = 1; 1309 if (Pf < 0) Pf = 1; 1310 if (Nc < 0) Nc = 1; 1311 if (Nf < 0) Nf = 1; 1312 1313 ierr = DMDAGetCorners(daf,&i_start,&j_start,&l_start,&m_f,&n_f,&p_f);CHKERRQ(ierr); 1314 ierr = DMDAGetGhostCorners(daf,&i_start_ghost,&j_start_ghost,&l_start_ghost,&m_ghost,&n_ghost,&p_ghost);CHKERRQ(ierr); 1315 ierr = DMDAGetGlobalIndices(daf,PETSC_NULL,&idx_f);CHKERRQ(ierr); 1316 1317 ierr = DMDAGetCorners(dac,&i_start_c,&j_start_c,&l_start_c,&m_c,&n_c,&p_c);CHKERRQ(ierr); 1318 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); 1319 ierr = DMDAGetGlobalIndices(dac,PETSC_NULL,&idx_c);CHKERRQ(ierr); 1320 1321 /* 1322 Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios 1323 for dimension 1 and 2 respectively. 1324 Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1) 1325 and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate. 1326 Each specific dof on the fine grid is mapped to one dof on the coarse grid. 1327 */ 1328 1329 max_agg_size = (Mf/Mc+1)*(Nf/Nc+1)*(Pf/Pc+1); 1330 1331 /* create the matrix that will contain the restriction operator */ 1332 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, 1333 max_agg_size, PETSC_NULL, max_agg_size, PETSC_NULL, rest);CHKERRQ(ierr); 1334 1335 /* store nodes in the fine grid here */ 1336 ierr = PetscMalloc2(max_agg_size,PetscScalar, &one_vec,max_agg_size,PetscInt, &fine_nodes);CHKERRQ(ierr); 1337 for (i=0; i<max_agg_size; i++) one_vec[i] = 1.0; 1338 1339 /* loop over all coarse nodes */ 1340 for (l_c=l_start_c; l_c<l_start_c+p_c; l_c++) { 1341 for (j_c=j_start_c; j_c<j_start_c+n_c; j_c++) { 1342 for (i_c=i_start_c; i_c<i_start_c+m_c; i_c++) { 1343 for (d=0; d<dofc; d++) { 1344 /* convert to local "natural" numbering and then to PETSc global numbering */ 1345 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; 1346 1347 fn_idx = 0; 1348 /* Corresponding fine points are all points (i_f, j_f, l_f) such that 1349 i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc 1350 (same for other dimensions) 1351 */ 1352 for (l=l_c*Pf/Pc; l<PetscMin((l_c+1)*Pf/Pc,Pf); l++) { 1353 for (j=j_c*Nf/Nc; j<PetscMin((j_c+1)*Nf/Nc,Nf); j++) { 1354 for (i=i_c*Mf/Mc; i<PetscMin((i_c+1)*Mf/Mc,Mf); i++) { 1355 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; 1356 fn_idx++; 1357 } 1358 } 1359 } 1360 /* add all these points to one aggregate */ 1361 ierr = MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES);CHKERRQ(ierr); 1362 } 1363 } 1364 } 1365 } 1366 ierr = PetscFree2(one_vec,fine_nodes);CHKERRQ(ierr); 1367 ierr = MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1368 ierr = MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1369 PetscFunctionReturn(0); 1370 } 1371