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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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++) 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; 271 } 272 273 PetscCall(DMDAVecRestoreArray(da, field, &FF)); 274 PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX)); 275 PetscFunctionReturn(PETSC_SUCCESS); 276 } 277 278 PetscErrorCode DADefineXLinearField3D(DM da, Vec field) 279 { 280 PetscInt i, j, k; 281 PetscInt sx, nx, sy, ny, sz, nz; 282 Vec Gcoords; 283 DMDACoor3d ***XX; 284 PetscScalar ***FF; 285 DM cda; 286 287 PetscFunctionBeginUser; 288 PetscCall(DMGetCoordinateDM(da, &cda)); 289 PetscCall(DMGetCoordinates(da, &Gcoords)); 290 291 PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX)); 292 PetscCall(DMDAVecGetArray(da, field, &FF)); 293 294 PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz)); 295 296 for (k = sz; k < sz + nz; k++) { 297 for (j = sy; j < sy + ny; j++) { 298 for (i = sx; i < sx + nx; i++) { 299 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 + 300 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z; 301 } 302 } 303 } 304 305 PetscCall(DMDAVecRestoreArray(da, field, &FF)); 306 PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX)); 307 PetscFunctionReturn(PETSC_SUCCESS); 308 } 309 310 PetscErrorCode da_test_RefineCoords1D(PetscInt mx) 311 { 312 DM dac, daf; 313 PetscViewer vv; 314 Vec ac, af; 315 PetscInt Mx; 316 Mat II, INTERP; 317 Vec scale; 318 PetscBool output = PETSC_FALSE; 319 320 PetscFunctionBeginUser; 321 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, mx + 1, 1, /* 1 dof */ 1, /* stencil = 1 */ NULL, &dac)); 322 PetscCall(DMSetFromOptions(dac)); 323 PetscCall(DMSetUp(dac)); 324 325 PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf)); 326 PetscCall(DMDAGetInfo(daf, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 327 Mx--; 328 329 PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE)); 330 PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE)); 331 332 { 333 DM cdaf, cdac; 334 Vec coordsc, coordsf; 335 336 PetscCall(DMGetCoordinateDM(dac, &cdac)); 337 PetscCall(DMGetCoordinateDM(daf, &cdaf)); 338 339 PetscCall(DMGetCoordinates(dac, &coordsc)); 340 PetscCall(DMGetCoordinates(daf, &coordsf)); 341 342 PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale)); 343 PetscCall(MatInterpolate(II, coordsc, coordsf)); 344 PetscCall(MatDestroy(&II)); 345 PetscCall(VecDestroy(&scale)); 346 } 347 348 PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL)); 349 350 PetscCall(DMCreateGlobalVector(dac, &ac)); 351 PetscCall(VecSet(ac, 66.99)); 352 353 PetscCall(DMCreateGlobalVector(daf, &af)); 354 355 PetscCall(MatMult(INTERP, ac, af)); 356 357 { 358 Vec afexact; 359 PetscReal nrm; 360 PetscInt N; 361 362 PetscCall(DMCreateGlobalVector(daf, &afexact)); 363 PetscCall(VecSet(afexact, 66.99)); 364 PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */ 365 PetscCall(VecNorm(afexact, NORM_2, &nrm)); 366 PetscCall(VecGetSize(afexact, &N)); 367 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "=>%" PetscInt_FMT ", interp err = %1.4e\n", mx, Mx, (double)(nrm / PetscSqrtReal((PetscReal)N)))); 368 PetscCall(VecDestroy(&afexact)); 369 } 370 371 PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL)); 372 if (output) { 373 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv)); 374 PetscCall(VecView(ac, vv)); 375 PetscCall(PetscViewerDestroy(&vv)); 376 377 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv)); 378 PetscCall(VecView(af, vv)); 379 PetscCall(PetscViewerDestroy(&vv)); 380 } 381 382 PetscCall(MatDestroy(&INTERP)); 383 PetscCall(DMDestroy(&dac)); 384 PetscCall(DMDestroy(&daf)); 385 PetscCall(VecDestroy(&ac)); 386 PetscCall(VecDestroy(&af)); 387 PetscFunctionReturn(PETSC_SUCCESS); 388 } 389 390 PetscErrorCode da_test_RefineCoords2D(PetscInt mx, PetscInt my) 391 { 392 DM dac, daf; 393 PetscViewer vv; 394 Vec ac, af; 395 PetscInt map_id, Mx, My; 396 Mat II, INTERP; 397 Vec scale; 398 PetscBool output = PETSC_FALSE; 399 400 PetscFunctionBeginUser; 401 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)); 402 PetscCall(DMSetFromOptions(dac)); 403 PetscCall(DMSetUp(dac)); 404 405 PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf)); 406 PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 407 Mx--; 408 My--; 409 410 PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE)); 411 PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE)); 412 413 /* apply conformal mappings */ 414 map_id = 0; 415 PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL)); 416 if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id)); 417 418 { 419 DM cdaf, cdac; 420 Vec coordsc, coordsf; 421 422 PetscCall(DMGetCoordinateDM(dac, &cdac)); 423 PetscCall(DMGetCoordinateDM(daf, &cdaf)); 424 425 PetscCall(DMGetCoordinates(dac, &coordsc)); 426 PetscCall(DMGetCoordinates(daf, &coordsf)); 427 428 PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale)); 429 PetscCall(MatInterpolate(II, coordsc, coordsf)); 430 PetscCall(MatDestroy(&II)); 431 PetscCall(VecDestroy(&scale)); 432 } 433 434 PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL)); 435 436 PetscCall(DMCreateGlobalVector(dac, &ac)); 437 PetscCall(DADefineXLinearField2D(dac, ac)); 438 439 PetscCall(DMCreateGlobalVector(daf, &af)); 440 PetscCall(MatMult(INTERP, ac, af)); 441 442 { 443 Vec afexact; 444 PetscReal nrm; 445 PetscInt N; 446 447 PetscCall(DMCreateGlobalVector(daf, &afexact)); 448 PetscCall(VecZeroEntries(afexact)); 449 PetscCall(DADefineXLinearField2D(daf, afexact)); 450 PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */ 451 PetscCall(VecNorm(afexact, NORM_2, &nrm)); 452 PetscCall(VecGetSize(afexact, &N)); 453 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)))); 454 PetscCall(VecDestroy(&afexact)); 455 } 456 457 PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL)); 458 if (output) { 459 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv)); 460 PetscCall(VecView(ac, vv)); 461 PetscCall(PetscViewerDestroy(&vv)); 462 463 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv)); 464 PetscCall(VecView(af, vv)); 465 PetscCall(PetscViewerDestroy(&vv)); 466 } 467 468 PetscCall(MatDestroy(&INTERP)); 469 PetscCall(DMDestroy(&dac)); 470 PetscCall(DMDestroy(&daf)); 471 PetscCall(VecDestroy(&ac)); 472 PetscCall(VecDestroy(&af)); 473 PetscFunctionReturn(PETSC_SUCCESS); 474 } 475 476 PetscErrorCode da_test_RefineCoords3D(PetscInt mx, PetscInt my, PetscInt mz) 477 { 478 DM dac, daf; 479 PetscViewer vv; 480 Vec ac, af; 481 PetscInt map_id, Mx, My, Mz; 482 Mat II, INTERP; 483 Vec scale; 484 PetscBool output = PETSC_FALSE; 485 486 PetscFunctionBeginUser; 487 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 */ 488 1, /* stencil = 1 */ NULL, NULL, NULL, &dac)); 489 PetscCall(DMSetFromOptions(dac)); 490 PetscCall(DMSetUp(dac)); 491 492 PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf)); 493 PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, &Mz, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 494 Mx--; 495 My--; 496 Mz--; 497 498 PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0)); 499 PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0)); 500 501 /* apply trilinear mappings */ 502 /*PetscCall(DAApplyTrilinearMapping(dac));*/ 503 /* apply conformal mappings */ 504 map_id = 0; 505 PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL)); 506 if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id)); 507 508 { 509 DM cdaf, cdac; 510 Vec coordsc, coordsf; 511 512 PetscCall(DMGetCoordinateDM(dac, &cdac)); 513 PetscCall(DMGetCoordinateDM(daf, &cdaf)); 514 515 PetscCall(DMGetCoordinates(dac, &coordsc)); 516 PetscCall(DMGetCoordinates(daf, &coordsf)); 517 518 PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale)); 519 PetscCall(MatInterpolate(II, coordsc, coordsf)); 520 PetscCall(MatDestroy(&II)); 521 PetscCall(VecDestroy(&scale)); 522 } 523 524 PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL)); 525 526 PetscCall(DMCreateGlobalVector(dac, &ac)); 527 PetscCall(VecZeroEntries(ac)); 528 PetscCall(DADefineXLinearField3D(dac, ac)); 529 530 PetscCall(DMCreateGlobalVector(daf, &af)); 531 PetscCall(VecZeroEntries(af)); 532 533 PetscCall(MatMult(INTERP, ac, af)); 534 535 { 536 Vec afexact; 537 PetscReal nrm; 538 PetscInt N; 539 540 PetscCall(DMCreateGlobalVector(daf, &afexact)); 541 PetscCall(VecZeroEntries(afexact)); 542 PetscCall(DADefineXLinearField3D(daf, afexact)); 543 PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */ 544 PetscCall(VecNorm(afexact, NORM_2, &nrm)); 545 PetscCall(VecGetSize(afexact, &N)); 546 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)))); 547 PetscCall(VecDestroy(&afexact)); 548 } 549 550 PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL)); 551 if (output) { 552 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv)); 553 PetscCall(VecView(ac, vv)); 554 PetscCall(PetscViewerDestroy(&vv)); 555 556 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv)); 557 PetscCall(VecView(af, vv)); 558 PetscCall(PetscViewerDestroy(&vv)); 559 } 560 561 PetscCall(MatDestroy(&INTERP)); 562 PetscCall(DMDestroy(&dac)); 563 PetscCall(DMDestroy(&daf)); 564 PetscCall(VecDestroy(&ac)); 565 PetscCall(VecDestroy(&af)); 566 PetscFunctionReturn(PETSC_SUCCESS); 567 } 568 569 int main(int argc, char **argv) 570 { 571 PetscInt mx = 2, my = 2, mz = 2, l, nl, dim; 572 573 PetscFunctionBeginUser; 574 PetscCall(PetscInitialize(&argc, &argv, 0, help)); 575 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0)); 576 PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, 0)); 577 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0)); 578 nl = 1; 579 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &nl, 0)); 580 dim = 2; 581 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, 0)); 582 583 for (l = 0; l < nl; l++) { 584 if (dim == 1) PetscCall(da_test_RefineCoords1D(mx)); 585 else if (dim == 2) PetscCall(da_test_RefineCoords2D(mx, my)); 586 else if (dim == 3) PetscCall(da_test_RefineCoords3D(mx, my, mz)); 587 mx = mx * 2; 588 my = my * 2; 589 mz = mz * 2; 590 } 591 PetscCall(PetscFinalize()); 592 return 0; 593 } 594 595 /*TEST 596 597 test: 598 suffix: 1d 599 args: -mx 10 -nl 6 -dim 1 600 601 test: 602 suffix: 2d 603 output_file: output/ex36_2d.out 604 args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}} 605 606 test: 607 suffix: 2dp1 608 nsize: 8 609 args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4 610 timeoutfactor: 2 611 612 test: 613 suffix: 2dp2 614 nsize: 8 615 args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1 616 timeoutfactor: 2 617 618 test: 619 suffix: 3d 620 args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3 621 622 test: 623 suffix: 3dp1 624 nsize: 32 625 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 626 627 TEST*/ 628