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