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 PetscCall(DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 70 } else if (idx==2) { /* stagnation in a corner */ 71 PetscCall(DMDASetUniformCoordinates(da, -1.0,1.0, 0.0,1.0, -1.0,1.0)); 72 } else if (idx==3) { /* nautilis */ 73 PetscCall(DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 74 } else if (idx==4) { 75 PetscCall(DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 76 } 77 78 PetscCall(DMGetCoordinateDM(da,&cda)); 79 PetscCall(DMGetCoordinates(da,&Gcoords)); 80 81 PetscCall(VecGetArray(Gcoords,&XX)); 82 PetscCall(DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz)); 83 PetscCall(DMDAGetInfo(da, &dim, 0,0,0, 0,0,0, 0,0,0,0,0,0)); 84 PetscCall(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 PetscCall(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 PetscCall(DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 209 PetscCall(DMGetCoordinateDM(da,&cda)); 210 PetscCall(DMGetCoordinates(da,&Gcoords)); 211 212 PetscCall(DMDAVecGetArrayRead(cda,Gcoords,&XX)); 213 PetscCall(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 PetscCall(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 PetscCall(DMGetCoordinateDM(da,&cda)); 264 PetscCall(DMGetCoordinates(da,&Gcoords)); 265 266 PetscCall(DMDAVecGetArrayRead(cda,Gcoords,&XX)); 267 PetscCall(DMDAVecGetArray(da,field,&FF)); 268 269 PetscCall(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 PetscCall(DMDAVecRestoreArray(da,field,&FF)); 278 PetscCall(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 PetscCall(DMGetCoordinateDM(da,&cda)); 293 PetscCall(DMGetCoordinates(da,&Gcoords)); 294 295 PetscCall(DMDAVecGetArrayRead(cda,Gcoords,&XX)); 296 PetscCall(DMDAVecGetArray(da,field,&FF)); 297 298 PetscCall(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 PetscCall(DMDAVecRestoreArray(da,field,&FF)); 316 PetscCall(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 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,mx+1,1, /* 1 dof */1, /* stencil = 1 */NULL,&dac)); 332 PetscCall(DMSetFromOptions(dac)); 333 PetscCall(DMSetUp(dac)); 334 335 PetscCall(DMRefine(dac,MPI_COMM_NULL,&daf)); 336 PetscCall(DMDAGetInfo(daf,0,&Mx,0,0,0,0,0,0,0,0,0,0,0)); 337 Mx--; 338 339 PetscCall(DMDASetUniformCoordinates(dac, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE)); 340 PetscCall(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 PetscCall(DMGetCoordinateDM(dac,&cdac)); 347 PetscCall(DMGetCoordinateDM(daf,&cdaf)); 348 349 PetscCall(DMGetCoordinates(dac,&coordsc)); 350 PetscCall(DMGetCoordinates(daf,&coordsf)); 351 352 PetscCall(DMCreateInterpolation(cdac,cdaf,&II,&scale)); 353 PetscCall(MatInterpolate(II,coordsc,coordsf)); 354 PetscCall(MatDestroy(&II)); 355 PetscCall(VecDestroy(&scale)); 356 } 357 358 PetscCall(DMCreateInterpolation(dac,daf,&INTERP,NULL)); 359 360 PetscCall(DMCreateGlobalVector(dac,&ac)); 361 PetscCall(VecSet(ac,66.99)); 362 363 PetscCall(DMCreateGlobalVector(daf,&af)); 364 365 PetscCall(MatMult(INTERP,ac, af)); 366 367 { 368 Vec afexact; 369 PetscReal nrm; 370 PetscInt N; 371 372 PetscCall(DMCreateGlobalVector(daf,&afexact)); 373 PetscCall(VecSet(afexact,66.99)); 374 PetscCall(VecAXPY(afexact,-1.0,af)); /* af <= af - afinterp */ 375 PetscCall(VecNorm(afexact,NORM_2,&nrm)); 376 PetscCall(VecGetSize(afexact,&N)); 377 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%D=>%D, interp err = %1.4e\n",mx,Mx,(double)(nrm/PetscSqrtReal((PetscReal)N)))); 378 PetscCall(VecDestroy(&afexact)); 379 } 380 381 PetscCall(PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL)); 382 if (output) { 383 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv)); 384 PetscCall(VecView(ac, vv)); 385 PetscCall(PetscViewerDestroy(&vv)); 386 387 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv)); 388 PetscCall(VecView(af, vv)); 389 PetscCall(PetscViewerDestroy(&vv)); 390 } 391 392 PetscCall(MatDestroy(&INTERP)); 393 PetscCall(DMDestroy(&dac)); 394 PetscCall(DMDestroy(&daf)); 395 PetscCall(VecDestroy(&ac)); 396 PetscCall(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 PetscCall(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 PetscCall(DMSetFromOptions(dac)); 413 PetscCall(DMSetUp(dac)); 414 415 PetscCall(DMRefine(dac,MPI_COMM_NULL,&daf)); 416 PetscCall(DMDAGetInfo(daf,0,&Mx,&My,0,0,0,0,0,0,0,0,0,0)); 417 Mx--; My--; 418 419 PetscCall(DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE)); 420 PetscCall(DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE)); 421 422 /* apply conformal mappings */ 423 map_id = 0; 424 PetscCall(PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL)); 425 if (map_id >= 1) { 426 PetscCall(DAApplyConformalMapping(dac,map_id)); 427 } 428 429 { 430 DM cdaf,cdac; 431 Vec coordsc,coordsf; 432 433 PetscCall(DMGetCoordinateDM(dac,&cdac)); 434 PetscCall(DMGetCoordinateDM(daf,&cdaf)); 435 436 PetscCall(DMGetCoordinates(dac,&coordsc)); 437 PetscCall(DMGetCoordinates(daf,&coordsf)); 438 439 PetscCall(DMCreateInterpolation(cdac,cdaf,&II,&scale)); 440 PetscCall(MatInterpolate(II,coordsc,coordsf)); 441 PetscCall(MatDestroy(&II)); 442 PetscCall(VecDestroy(&scale)); 443 } 444 445 PetscCall(DMCreateInterpolation(dac,daf,&INTERP,NULL)); 446 447 PetscCall(DMCreateGlobalVector(dac,&ac)); 448 PetscCall(DADefineXLinearField2D(dac,ac)); 449 450 PetscCall(DMCreateGlobalVector(daf,&af)); 451 PetscCall(MatMult(INTERP,ac, af)); 452 453 { 454 Vec afexact; 455 PetscReal nrm; 456 PetscInt N; 457 458 PetscCall(DMCreateGlobalVector(daf,&afexact)); 459 PetscCall(VecZeroEntries(afexact)); 460 PetscCall(DADefineXLinearField2D(daf,afexact)); 461 PetscCall(VecAXPY(afexact,-1.0,af)); /* af <= af - afinterp */ 462 PetscCall(VecNorm(afexact,NORM_2,&nrm)); 463 PetscCall(VecGetSize(afexact,&N)); 464 PetscCall(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 PetscCall(VecDestroy(&afexact)); 466 } 467 468 PetscCall(PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL)); 469 if (output) { 470 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv)); 471 PetscCall(VecView(ac, vv)); 472 PetscCall(PetscViewerDestroy(&vv)); 473 474 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv)); 475 PetscCall(VecView(af, vv)); 476 PetscCall(PetscViewerDestroy(&vv)); 477 } 478 479 PetscCall(MatDestroy(&INTERP)); 480 PetscCall(DMDestroy(&dac)); 481 PetscCall(DMDestroy(&daf)); 482 PetscCall(VecDestroy(&ac)); 483 PetscCall(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);PetscCall(ierr); 501 PetscCall(DMSetFromOptions(dac)); 502 PetscCall(DMSetUp(dac)); 503 504 PetscCall(DMRefine(dac,MPI_COMM_NULL,&daf)); 505 PetscCall(DMDAGetInfo(daf,0,&Mx,&My,&Mz,0,0,0,0,0,0,0,0,0)); 506 Mx--; My--; Mz--; 507 508 PetscCall(DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 509 PetscCall(DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0)); 510 511 /* apply trilinear mappings */ 512 /*PetscCall(DAApplyTrilinearMapping(dac));*/ 513 /* apply conformal mappings */ 514 map_id = 0; 515 PetscCall(PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL)); 516 if (map_id >= 1) { 517 PetscCall(DAApplyConformalMapping(dac,map_id)); 518 } 519 520 { 521 DM cdaf,cdac; 522 Vec coordsc,coordsf; 523 524 PetscCall(DMGetCoordinateDM(dac,&cdac)); 525 PetscCall(DMGetCoordinateDM(daf,&cdaf)); 526 527 PetscCall(DMGetCoordinates(dac,&coordsc)); 528 PetscCall(DMGetCoordinates(daf,&coordsf)); 529 530 PetscCall(DMCreateInterpolation(cdac,cdaf,&II,&scale)); 531 PetscCall(MatInterpolate(II,coordsc,coordsf)); 532 PetscCall(MatDestroy(&II)); 533 PetscCall(VecDestroy(&scale)); 534 } 535 536 PetscCall(DMCreateInterpolation(dac,daf,&INTERP,NULL)); 537 538 PetscCall(DMCreateGlobalVector(dac,&ac)); 539 PetscCall(VecZeroEntries(ac)); 540 PetscCall(DADefineXLinearField3D(dac,ac)); 541 542 PetscCall(DMCreateGlobalVector(daf,&af)); 543 PetscCall(VecZeroEntries(af)); 544 545 PetscCall(MatMult(INTERP,ac, af)); 546 547 { 548 Vec afexact; 549 PetscReal nrm; 550 PetscInt N; 551 552 PetscCall(DMCreateGlobalVector(daf,&afexact)); 553 PetscCall(VecZeroEntries(afexact)); 554 PetscCall(DADefineXLinearField3D(daf,afexact)); 555 PetscCall(VecAXPY(afexact,-1.0,af)); /* af <= af - afinterp */ 556 PetscCall(VecNorm(afexact,NORM_2,&nrm)); 557 PetscCall(VecGetSize(afexact,&N)); 558 PetscCall(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 PetscCall(VecDestroy(&afexact)); 560 } 561 562 PetscCall(PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL)); 563 if (output) { 564 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv)); 565 PetscCall(VecView(ac, vv)); 566 PetscCall(PetscViewerDestroy(&vv)); 567 568 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv)); 569 PetscCall(VecView(af, vv)); 570 PetscCall(PetscViewerDestroy(&vv)); 571 } 572 573 PetscCall(MatDestroy(&INTERP)); 574 PetscCall(DMDestroy(&dac)); 575 PetscCall(DMDestroy(&daf)); 576 PetscCall(VecDestroy(&ac)); 577 PetscCall(VecDestroy(&af)); 578 PetscFunctionReturn(0); 579 } 580 581 int main(int argc,char **argv) 582 { 583 PetscInt mx = 2,my = 2,mz = 2,l,nl,dim; 584 585 PetscCall(PetscInitialize(&argc,&argv,0,help)); 586 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mx", &mx, 0)); 587 PetscCall(PetscOptionsGetInt(NULL,NULL,"-my", &my, 0)); 588 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mz", &mz, 0)); 589 nl = 1; 590 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nl", &nl, 0)); 591 dim = 2; 592 PetscCall(PetscOptionsGetInt(NULL,NULL,"-dim", &dim, 0)); 593 594 for (l=0; l<nl; l++) { 595 if (dim==1) { 596 PetscCall(da_test_RefineCoords1D(mx)); 597 } else if (dim==2) { 598 PetscCall(da_test_RefineCoords2D(mx,my)); 599 } else if (dim==3) { 600 PetscCall(da_test_RefineCoords3D(mx,my,mz)); 601 } 602 mx = mx * 2; 603 my = my * 2; 604 mz = mz * 2; 605 } 606 PetscCall(PetscFinalize()); 607 return 0; 608 } 609 610 /*TEST 611 612 test: 613 suffix: 1d 614 args: -mx 10 -nl 6 -dim 1 615 616 test: 617 suffix: 2d 618 output_file: output/ex36_2d.out 619 args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}} 620 621 test: 622 suffix: 2dp1 623 nsize: 8 624 args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4 625 timeoutfactor: 2 626 627 test: 628 suffix: 2dp2 629 nsize: 8 630 args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1 631 timeoutfactor: 2 632 633 test: 634 suffix: 3d 635 args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3 636 637 test: 638 suffix: 3dp1 639 nsize: 32 640 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 641 642 TEST*/ 643