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
CCmplxPow(CCmplx a,PetscReal n)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 }
CCmplxExp(CCmplx a)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 }
CCmplxSqrt(CCmplx a)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 }
CCmplxAdd(CCmplx a,CCmplx c)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 }
CCmplxRe(CCmplx a)46 PetscScalar CCmplxRe(CCmplx a)
47 {
48 return a.real;
49 }
CCmplxIm(CCmplx a)50 PetscScalar CCmplxIm(CCmplx a)
51 {
52 return a.imag;
53 }
54
DAApplyConformalMapping(DM da,PetscInt idx)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
DAApplyTrilinearMapping(DM da)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
DADefineXLinearField2D(DM da,Vec field)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
DADefineXLinearField3D(DM da,Vec field)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
da_test_RefineCoords1D(PetscInt mx)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
da_test_RefineCoords2D(PetscInt mx,PetscInt my)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
da_test_RefineCoords3D(PetscInt mx,PetscInt my,PetscInt mz)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
main(int argc,char ** argv)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