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