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