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