1 2 static char help[] = "Checks the functionality of DMGetInterpolation() on deformed grids.\n\n"; 3 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 7 typedef struct _n_CCmplx CCmplx; 8 struct _n_CCmplx { 9 PetscReal real; 10 PetscReal imag; 11 }; 12 13 CCmplx CCmplxPow(CCmplx a,PetscReal n) 14 { 15 CCmplx b; 16 PetscReal r,theta; 17 r = PetscSqrtReal(a.real*a.real + a.imag*a.imag); 18 theta = PetscAtan2Real(a.imag,a.real); 19 b.real = PetscPowReal(r,n) * PetscCosReal(n*theta); 20 b.imag = PetscPowReal(r,n) * PetscSinReal(n*theta); 21 return b; 22 } 23 CCmplx CCmplxExp(CCmplx a) 24 { 25 CCmplx b; 26 b.real = PetscExpReal(a.real) * PetscCosReal(a.imag); 27 b.imag = PetscExpReal(a.real) * PetscSinReal(a.imag); 28 return b; 29 } 30 CCmplx CCmplxSqrt(CCmplx a) 31 { 32 CCmplx b; 33 PetscReal r,theta; 34 r = PetscSqrtReal(a.real*a.real + a.imag*a.imag); 35 theta = PetscAtan2Real(a.imag,a.real); 36 b.real = PetscSqrtReal(r) * PetscCosReal(0.5*theta); 37 b.imag = PetscSqrtReal(r) * PetscSinReal(0.5*theta); 38 return b; 39 } 40 CCmplx CCmplxAdd(CCmplx a,CCmplx c) 41 { 42 CCmplx b; 43 b.real = a.real +c.real; 44 b.imag = a.imag +c.imag; 45 return b; 46 } 47 PetscScalar CCmplxRe(CCmplx a) 48 { 49 return (PetscScalar)a.real; 50 } 51 PetscScalar CCmplxIm(CCmplx a) 52 { 53 return (PetscScalar)a.imag; 54 } 55 56 PetscErrorCode DAApplyConformalMapping(DM da,PetscInt idx) 57 { 58 PetscInt i,n; 59 PetscInt sx,nx,sy,ny,sz,nz,dim; 60 Vec Gcoords; 61 PetscScalar *XX; 62 PetscScalar xx,yy,zz; 63 DM cda; 64 65 PetscFunctionBeginUser; 66 if (idx==0) { 67 PetscFunctionReturn(0); 68 } else if (idx==1) { /* dam break */ 69 CHKERRQ(DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 70 } else if (idx==2) { /* stagnation in a corner */ 71 CHKERRQ(DMDASetUniformCoordinates(da, -1.0,1.0, 0.0,1.0, -1.0,1.0)); 72 } else if (idx==3) { /* nautilis */ 73 CHKERRQ(DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 74 } else if (idx==4) { 75 CHKERRQ(DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 76 } 77 78 CHKERRQ(DMGetCoordinateDM(da,&cda)); 79 CHKERRQ(DMGetCoordinates(da,&Gcoords)); 80 81 CHKERRQ(VecGetArray(Gcoords,&XX)); 82 CHKERRQ(DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz)); 83 CHKERRQ(DMDAGetInfo(da, &dim, 0,0,0, 0,0,0, 0,0,0,0,0,0)); 84 CHKERRQ(VecGetLocalSize(Gcoords,&n)); 85 n = n / dim; 86 87 for (i=0; i<n; i++) { 88 if ((dim==3) && (idx!=2)) { 89 PetscScalar Ni[8]; 90 PetscScalar xi = XX[dim*i]; 91 PetscScalar eta = XX[dim*i+1]; 92 PetscScalar zeta = XX[dim*i+2]; 93 PetscScalar xn[] = {-1.0,1.0,-1.0,1.0, -1.0,1.0,-1.0,1.0 }; 94 PetscScalar yn[] = {-1.0,-1.0,1.0,1.0, -1.0,-1.0,1.0,1.0 }; 95 PetscScalar zn[] = {-0.1,-4.0,-0.2,-1.0, 0.1,4.0,0.2,1.0 }; 96 PetscInt p; 97 98 Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta); 99 Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta); 100 Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta); 101 Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta); 102 103 Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta); 104 Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta); 105 Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta); 106 Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta); 107 108 xx = yy = zz = 0.0; 109 for (p=0; p<8; p++) { 110 xx += Ni[p]*xn[p]; 111 yy += Ni[p]*yn[p]; 112 zz += Ni[p]*zn[p]; 113 } 114 XX[dim*i] = xx; 115 XX[dim*i+1] = yy; 116 XX[dim*i+2] = zz; 117 } 118 119 if (idx==1) { 120 CCmplx zeta,t1,t2; 121 122 xx = XX[dim*i] - 0.8; 123 yy = XX[dim*i+1] + 1.5; 124 125 zeta.real = PetscRealPart(xx); 126 zeta.imag = PetscRealPart(yy); 127 128 t1 = CCmplxPow(zeta,-1.0); 129 t2 = CCmplxAdd(zeta,t1); 130 131 XX[dim*i] = CCmplxRe(t2); 132 XX[dim*i+1] = CCmplxIm(t2); 133 } else if (idx==2) { 134 CCmplx zeta,t1; 135 136 xx = XX[dim*i]; 137 yy = XX[dim*i+1]; 138 zeta.real = PetscRealPart(xx); 139 zeta.imag = PetscRealPart(yy); 140 141 t1 = CCmplxSqrt(zeta); 142 XX[dim*i] = CCmplxRe(t1); 143 XX[dim*i+1] = CCmplxIm(t1); 144 } else if (idx==3) { 145 CCmplx zeta,t1,t2; 146 147 xx = XX[dim*i] - 0.8; 148 yy = XX[dim*i+1] + 1.5; 149 150 zeta.real = PetscRealPart(xx); 151 zeta.imag = PetscRealPart(yy); 152 t1 = CCmplxPow(zeta,-1.0); 153 t2 = CCmplxAdd(zeta,t1); 154 XX[dim*i] = CCmplxRe(t2); 155 XX[dim*i+1] = CCmplxIm(t2); 156 157 xx = XX[dim*i]; 158 yy = XX[dim*i+1]; 159 zeta.real = PetscRealPart(xx); 160 zeta.imag = PetscRealPart(yy); 161 t1 = CCmplxExp(zeta); 162 XX[dim*i] = CCmplxRe(t1); 163 XX[dim*i+1] = CCmplxIm(t1); 164 165 xx = XX[dim*i] + 0.4; 166 yy = XX[dim*i+1]; 167 zeta.real = PetscRealPart(xx); 168 zeta.imag = PetscRealPart(yy); 169 t1 = CCmplxPow(zeta,2.0); 170 XX[dim*i] = CCmplxRe(t1); 171 XX[dim*i+1] = CCmplxIm(t1); 172 } else if (idx==4) { 173 PetscScalar Ni[4]; 174 PetscScalar xi = XX[dim*i]; 175 PetscScalar eta = XX[dim*i+1]; 176 PetscScalar xn[] = {0.0,2.0,0.2,3.5}; 177 PetscScalar yn[] = {-1.3,0.0,2.0,4.0}; 178 PetscInt p; 179 180 Ni[0] = 0.25*(1.0-xi)*(1.0-eta); 181 Ni[1] = 0.25*(1.0+xi)*(1.0-eta); 182 Ni[2] = 0.25*(1.0-xi)*(1.0+eta); 183 Ni[3] = 0.25*(1.0+xi)*(1.0+eta); 184 185 xx = yy = 0.0; 186 for (p=0; p<4; p++) { 187 xx += Ni[p]*xn[p]; 188 yy += Ni[p]*yn[p]; 189 } 190 XX[dim*i] = xx; 191 XX[dim*i+1] = yy; 192 } 193 } 194 CHKERRQ(VecRestoreArray(Gcoords,&XX)); 195 PetscFunctionReturn(0); 196 } 197 198 PetscErrorCode DAApplyTrilinearMapping(DM da) 199 { 200 PetscInt i,j,k; 201 PetscInt sx,nx,sy,ny,sz,nz; 202 Vec Gcoords; 203 DMDACoor3d ***XX; 204 PetscScalar xx,yy,zz; 205 DM cda; 206 207 PetscFunctionBeginUser; 208 CHKERRQ(DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 209 CHKERRQ(DMGetCoordinateDM(da,&cda)); 210 CHKERRQ(DMGetCoordinates(da,&Gcoords)); 211 212 CHKERRQ(DMDAVecGetArrayRead(cda,Gcoords,&XX)); 213 CHKERRQ(DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz)); 214 215 for (i=sx; i<sx+nx; i++) { 216 for (j=sy; j<sy+ny; j++) { 217 for (k=sz; k<sz+nz; k++) { 218 PetscScalar Ni[8]; 219 PetscScalar xi = XX[k][j][i].x; 220 PetscScalar eta = XX[k][j][i].y; 221 PetscScalar zeta = XX[k][j][i].z; 222 PetscScalar xn[] = {0.0,2.0,0.2,3.5, 0.0,2.1,0.23,3.125 }; 223 PetscScalar yn[] = {-1.3,0.0,2.0,4.0, -1.45,-0.1,2.24,3.79 }; 224 PetscScalar zn[] = {0.0,0.3,-0.1,0.123, 0.956,1.32,1.12,0.798 }; 225 PetscInt p; 226 227 Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta); 228 Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta); 229 Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta); 230 Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta); 231 232 Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta); 233 Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta); 234 Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta); 235 Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta); 236 237 xx = yy = zz = 0.0; 238 for (p=0; p<8; p++) { 239 xx += Ni[p]*xn[p]; 240 yy += Ni[p]*yn[p]; 241 zz += Ni[p]*zn[p]; 242 } 243 XX[k][j][i].x = xx; 244 XX[k][j][i].y = yy; 245 XX[k][j][i].z = zz; 246 } 247 } 248 } 249 CHKERRQ(DMDAVecRestoreArrayRead(cda,Gcoords,&XX)); 250 PetscFunctionReturn(0); 251 } 252 253 PetscErrorCode DADefineXLinearField2D(DM da,Vec field) 254 { 255 PetscInt i,j; 256 PetscInt sx,nx,sy,ny; 257 Vec Gcoords; 258 DMDACoor2d **XX; 259 PetscScalar **FF; 260 DM cda; 261 262 PetscFunctionBeginUser; 263 CHKERRQ(DMGetCoordinateDM(da,&cda)); 264 CHKERRQ(DMGetCoordinates(da,&Gcoords)); 265 266 CHKERRQ(DMDAVecGetArrayRead(cda,Gcoords,&XX)); 267 CHKERRQ(DMDAVecGetArray(da,field,&FF)); 268 269 CHKERRQ(DMDAGetCorners(da,&sx,&sy,0,&nx,&ny,0)); 270 271 for (i=sx; i<sx+nx; i++) { 272 for (j=sy; j<sy+ny; j++) { 273 FF[j][i] = 10.0 + 3.0 * XX[j][i].x + 5.5 * XX[j][i].y + 8.003 * XX[j][i].x * XX[j][i].y; 274 } 275 } 276 277 CHKERRQ(DMDAVecRestoreArray(da,field,&FF)); 278 CHKERRQ(DMDAVecRestoreArrayRead(cda,Gcoords,&XX)); 279 PetscFunctionReturn(0); 280 } 281 282 PetscErrorCode DADefineXLinearField3D(DM da,Vec field) 283 { 284 PetscInt i,j,k; 285 PetscInt sx,nx,sy,ny,sz,nz; 286 Vec Gcoords; 287 DMDACoor3d ***XX; 288 PetscScalar ***FF; 289 DM cda; 290 291 PetscFunctionBeginUser; 292 CHKERRQ(DMGetCoordinateDM(da,&cda)); 293 CHKERRQ(DMGetCoordinates(da,&Gcoords)); 294 295 CHKERRQ(DMDAVecGetArrayRead(cda,Gcoords,&XX)); 296 CHKERRQ(DMDAVecGetArray(da,field,&FF)); 297 298 CHKERRQ(DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz)); 299 300 for (k=sz; k<sz+nz; k++) { 301 for (j=sy; j<sy+ny; j++) { 302 for (i=sx; i<sx+nx; i++) { 303 FF[k][j][i] = 10.0 304 + 4.05 * XX[k][j][i].x 305 + 5.50 * XX[k][j][i].y 306 + 1.33 * XX[k][j][i].z 307 + 2.03 * XX[k][j][i].x * XX[k][j][i].y 308 + 0.03 * XX[k][j][i].x * XX[k][j][i].z 309 + 0.83 * XX[k][j][i].y * XX[k][j][i].z 310 + 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z; 311 } 312 } 313 } 314 315 CHKERRQ(DMDAVecRestoreArray(da,field,&FF)); 316 CHKERRQ(DMDAVecRestoreArrayRead(cda,Gcoords,&XX)); 317 PetscFunctionReturn(0); 318 } 319 320 PetscErrorCode da_test_RefineCoords1D(PetscInt mx) 321 { 322 DM dac,daf; 323 PetscViewer vv; 324 Vec ac,af; 325 PetscInt Mx; 326 Mat II,INTERP; 327 Vec scale; 328 PetscBool output = PETSC_FALSE; 329 330 PetscFunctionBeginUser; 331 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,mx+1,1, /* 1 dof */1, /* stencil = 1 */NULL,&dac)); 332 CHKERRQ(DMSetFromOptions(dac)); 333 CHKERRQ(DMSetUp(dac)); 334 335 CHKERRQ(DMRefine(dac,MPI_COMM_NULL,&daf)); 336 CHKERRQ(DMDAGetInfo(daf,0,&Mx,0,0,0,0,0,0,0,0,0,0,0)); 337 Mx--; 338 339 CHKERRQ(DMDASetUniformCoordinates(dac, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE)); 340 CHKERRQ(DMDASetUniformCoordinates(daf, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE)); 341 342 { 343 DM cdaf,cdac; 344 Vec coordsc,coordsf; 345 346 CHKERRQ(DMGetCoordinateDM(dac,&cdac)); 347 CHKERRQ(DMGetCoordinateDM(daf,&cdaf)); 348 349 CHKERRQ(DMGetCoordinates(dac,&coordsc)); 350 CHKERRQ(DMGetCoordinates(daf,&coordsf)); 351 352 CHKERRQ(DMCreateInterpolation(cdac,cdaf,&II,&scale)); 353 CHKERRQ(MatInterpolate(II,coordsc,coordsf)); 354 CHKERRQ(MatDestroy(&II)); 355 CHKERRQ(VecDestroy(&scale)); 356 } 357 358 CHKERRQ(DMCreateInterpolation(dac,daf,&INTERP,NULL)); 359 360 CHKERRQ(DMCreateGlobalVector(dac,&ac)); 361 CHKERRQ(VecSet(ac,66.99)); 362 363 CHKERRQ(DMCreateGlobalVector(daf,&af)); 364 365 CHKERRQ(MatMult(INTERP,ac, af)); 366 367 { 368 Vec afexact; 369 PetscReal nrm; 370 PetscInt N; 371 372 CHKERRQ(DMCreateGlobalVector(daf,&afexact)); 373 CHKERRQ(VecSet(afexact,66.99)); 374 CHKERRQ(VecAXPY(afexact,-1.0,af)); /* af <= af - afinterp */ 375 CHKERRQ(VecNorm(afexact,NORM_2,&nrm)); 376 CHKERRQ(VecGetSize(afexact,&N)); 377 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%D=>%D, interp err = %1.4e\n",mx,Mx,(double)(nrm/PetscSqrtReal((PetscReal)N)))); 378 CHKERRQ(VecDestroy(&afexact)); 379 } 380 381 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL)); 382 if (output) { 383 CHKERRQ(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv)); 384 CHKERRQ(VecView(ac, vv)); 385 CHKERRQ(PetscViewerDestroy(&vv)); 386 387 CHKERRQ(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv)); 388 CHKERRQ(VecView(af, vv)); 389 CHKERRQ(PetscViewerDestroy(&vv)); 390 } 391 392 CHKERRQ(MatDestroy(&INTERP)); 393 CHKERRQ(DMDestroy(&dac)); 394 CHKERRQ(DMDestroy(&daf)); 395 CHKERRQ(VecDestroy(&ac)); 396 CHKERRQ(VecDestroy(&af)); 397 PetscFunctionReturn(0); 398 } 399 400 PetscErrorCode da_test_RefineCoords2D(PetscInt mx,PetscInt my) 401 { 402 DM dac,daf; 403 PetscViewer vv; 404 Vec ac,af; 405 PetscInt map_id,Mx,My; 406 Mat II,INTERP; 407 Vec scale; 408 PetscBool output = PETSC_FALSE; 409 410 PetscFunctionBeginUser; 411 CHKERRQ(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,mx+1,my+1,PETSC_DECIDE, PETSC_DECIDE,1, /* 1 dof */1, /* stencil = 1 */NULL, NULL,&dac)); 412 CHKERRQ(DMSetFromOptions(dac)); 413 CHKERRQ(DMSetUp(dac)); 414 415 CHKERRQ(DMRefine(dac,MPI_COMM_NULL,&daf)); 416 CHKERRQ(DMDAGetInfo(daf,0,&Mx,&My,0,0,0,0,0,0,0,0,0,0)); 417 Mx--; My--; 418 419 CHKERRQ(DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE)); 420 CHKERRQ(DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE)); 421 422 /* apply conformal mappings */ 423 map_id = 0; 424 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL)); 425 if (map_id >= 1) { 426 CHKERRQ(DAApplyConformalMapping(dac,map_id)); 427 } 428 429 { 430 DM cdaf,cdac; 431 Vec coordsc,coordsf; 432 433 CHKERRQ(DMGetCoordinateDM(dac,&cdac)); 434 CHKERRQ(DMGetCoordinateDM(daf,&cdaf)); 435 436 CHKERRQ(DMGetCoordinates(dac,&coordsc)); 437 CHKERRQ(DMGetCoordinates(daf,&coordsf)); 438 439 CHKERRQ(DMCreateInterpolation(cdac,cdaf,&II,&scale)); 440 CHKERRQ(MatInterpolate(II,coordsc,coordsf)); 441 CHKERRQ(MatDestroy(&II)); 442 CHKERRQ(VecDestroy(&scale)); 443 } 444 445 CHKERRQ(DMCreateInterpolation(dac,daf,&INTERP,NULL)); 446 447 CHKERRQ(DMCreateGlobalVector(dac,&ac)); 448 CHKERRQ(DADefineXLinearField2D(dac,ac)); 449 450 CHKERRQ(DMCreateGlobalVector(daf,&af)); 451 CHKERRQ(MatMult(INTERP,ac, af)); 452 453 { 454 Vec afexact; 455 PetscReal nrm; 456 PetscInt N; 457 458 CHKERRQ(DMCreateGlobalVector(daf,&afexact)); 459 CHKERRQ(VecZeroEntries(afexact)); 460 CHKERRQ(DADefineXLinearField2D(daf,afexact)); 461 CHKERRQ(VecAXPY(afexact,-1.0,af)); /* af <= af - afinterp */ 462 CHKERRQ(VecNorm(afexact,NORM_2,&nrm)); 463 CHKERRQ(VecGetSize(afexact,&N)); 464 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[%D x %D]=>[%D x %D], interp err = %1.4e\n",mx,my,Mx,My,(double)(nrm/PetscSqrtReal((PetscReal)N)))); 465 CHKERRQ(VecDestroy(&afexact)); 466 } 467 468 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL)); 469 if (output) { 470 CHKERRQ(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv)); 471 CHKERRQ(VecView(ac, vv)); 472 CHKERRQ(PetscViewerDestroy(&vv)); 473 474 CHKERRQ(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv)); 475 CHKERRQ(VecView(af, vv)); 476 CHKERRQ(PetscViewerDestroy(&vv)); 477 } 478 479 CHKERRQ(MatDestroy(&INTERP)); 480 CHKERRQ(DMDestroy(&dac)); 481 CHKERRQ(DMDestroy(&daf)); 482 CHKERRQ(VecDestroy(&ac)); 483 CHKERRQ(VecDestroy(&af)); 484 PetscFunctionReturn(0); 485 } 486 487 PetscErrorCode da_test_RefineCoords3D(PetscInt mx,PetscInt my,PetscInt mz) 488 { 489 PetscErrorCode ierr; 490 DM dac,daf; 491 PetscViewer vv; 492 Vec ac,af; 493 PetscInt map_id,Mx,My,Mz; 494 Mat II,INTERP; 495 Vec scale; 496 PetscBool output = PETSC_FALSE; 497 498 PetscFunctionBeginUser; 499 ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx+1, my+1,mz+1,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,1, /* 1 dof */ 500 1, /* stencil = 1 */NULL,NULL,NULL,&dac);CHKERRQ(ierr); 501 CHKERRQ(DMSetFromOptions(dac)); 502 CHKERRQ(DMSetUp(dac)); 503 504 CHKERRQ(DMRefine(dac,MPI_COMM_NULL,&daf)); 505 CHKERRQ(DMDAGetInfo(daf,0,&Mx,&My,&Mz,0,0,0,0,0,0,0,0,0)); 506 Mx--; My--; Mz--; 507 508 CHKERRQ(DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 509 CHKERRQ(DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 510 511 /* apply trilinear mappings */ 512 /*CHKERRQ(DAApplyTrilinearMapping(dac));*/ 513 /* apply conformal mappings */ 514 map_id = 0; 515 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL)); 516 if (map_id >= 1) { 517 CHKERRQ(DAApplyConformalMapping(dac,map_id)); 518 } 519 520 { 521 DM cdaf,cdac; 522 Vec coordsc,coordsf; 523 524 CHKERRQ(DMGetCoordinateDM(dac,&cdac)); 525 CHKERRQ(DMGetCoordinateDM(daf,&cdaf)); 526 527 CHKERRQ(DMGetCoordinates(dac,&coordsc)); 528 CHKERRQ(DMGetCoordinates(daf,&coordsf)); 529 530 CHKERRQ(DMCreateInterpolation(cdac,cdaf,&II,&scale)); 531 CHKERRQ(MatInterpolate(II,coordsc,coordsf)); 532 CHKERRQ(MatDestroy(&II)); 533 CHKERRQ(VecDestroy(&scale)); 534 } 535 536 CHKERRQ(DMCreateInterpolation(dac,daf,&INTERP,NULL)); 537 538 CHKERRQ(DMCreateGlobalVector(dac,&ac)); 539 CHKERRQ(VecZeroEntries(ac)); 540 CHKERRQ(DADefineXLinearField3D(dac,ac)); 541 542 CHKERRQ(DMCreateGlobalVector(daf,&af)); 543 CHKERRQ(VecZeroEntries(af)); 544 545 CHKERRQ(MatMult(INTERP,ac, af)); 546 547 { 548 Vec afexact; 549 PetscReal nrm; 550 PetscInt N; 551 552 CHKERRQ(DMCreateGlobalVector(daf,&afexact)); 553 CHKERRQ(VecZeroEntries(afexact)); 554 CHKERRQ(DADefineXLinearField3D(daf,afexact)); 555 CHKERRQ(VecAXPY(afexact,-1.0,af)); /* af <= af - afinterp */ 556 CHKERRQ(VecNorm(afexact,NORM_2,&nrm)); 557 CHKERRQ(VecGetSize(afexact,&N)); 558 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[%D x %D x %D]=>[%D x %D x %D], interp err = %1.4e\n",mx,my,mz,Mx,My,Mz,(double)(nrm/PetscSqrtReal((PetscReal)N)))); 559 CHKERRQ(VecDestroy(&afexact)); 560 } 561 562 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL)); 563 if (output) { 564 CHKERRQ(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv)); 565 CHKERRQ(VecView(ac, vv)); 566 CHKERRQ(PetscViewerDestroy(&vv)); 567 568 CHKERRQ(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv)); 569 CHKERRQ(VecView(af, vv)); 570 CHKERRQ(PetscViewerDestroy(&vv)); 571 } 572 573 CHKERRQ(MatDestroy(&INTERP)); 574 CHKERRQ(DMDestroy(&dac)); 575 CHKERRQ(DMDestroy(&daf)); 576 CHKERRQ(VecDestroy(&ac)); 577 CHKERRQ(VecDestroy(&af)); 578 PetscFunctionReturn(0); 579 } 580 581 int main(int argc,char **argv) 582 { 583 PetscErrorCode ierr; 584 PetscInt mx = 2,my = 2,mz = 2,l,nl,dim; 585 586 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 587 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mx", &mx, 0)); 588 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-my", &my, 0)); 589 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mz", &mz, 0)); 590 nl = 1; 591 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-nl", &nl, 0)); 592 dim = 2; 593 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-dim", &dim, 0)); 594 595 for (l=0; l<nl; l++) { 596 if (dim==1) { 597 CHKERRQ(da_test_RefineCoords1D(mx)); 598 } else if (dim==2) { 599 CHKERRQ(da_test_RefineCoords2D(mx,my)); 600 } else if (dim==3) { 601 CHKERRQ(da_test_RefineCoords3D(mx,my,mz)); 602 } 603 mx = mx * 2; 604 my = my * 2; 605 mz = mz * 2; 606 } 607 ierr = PetscFinalize(); 608 return ierr; 609 } 610 611 /*TEST 612 613 test: 614 suffix: 1d 615 args: -mx 10 -nl 6 -dim 1 616 617 test: 618 suffix: 2d 619 output_file: output/ex36_2d.out 620 args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}} 621 622 test: 623 suffix: 2dp1 624 nsize: 8 625 args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4 626 timeoutfactor: 2 627 628 test: 629 suffix: 2dp2 630 nsize: 8 631 args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1 632 timeoutfactor: 2 633 634 test: 635 suffix: 3d 636 args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3 637 638 test: 639 suffix: 3dp1 640 nsize: 32 641 args: -mx 5 -my 5 -mz 5 -nl 3 -dim 3 -cmap 1 -da_refine_x 1 -da_refine_y 3 -da_refine_z 4 642 643 TEST*/ 644