-
Notifications
You must be signed in to change notification settings - Fork 2
/
interptest.py
100 lines (82 loc) · 2.97 KB
/
interptest.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
# -*- coding: utf-8; -*-
"""Investigate the interpolation of different element types on the unit square.
The degree can be specified:
python -m demo.interptest 1
python -m demo.interptest 2
python -m demo.interptest 3
Note this demo exercises also `quad_to_tri` and the quad-plotting functionality
of `extrafeathers`.
"""
import sys
from unpythonic import flatten1
import numpy as np
import matplotlib.pyplot as plt
import dolfin
from extrafeathers import plotmagic
# --------------------------------------------------------------------------------
# Set up the mesh
N = 2
arg = sys.argv[1] if len(sys.argv) > 1 else "1"
degree = int(arg)
# data = [+1.0, -1.0, -1.0, +1.0]
data = [0.0, 0.0, 0.0, 1.0]
# data = [0.0, 1.0, 1.0, 1.0]
# data = [0.0, 9.0, 4.0, 1.0]
# def bilerp(data, x):
# """
# `data`: length-4 vector, layout:
# 2-3
# | |
# 0-1
# `x`: length-2 vector, point in unit square
# """
# return (data[0] * (1 - x[0]) * (1 - x[1]) +
# data[1] * x[0] * (1 - x[1]) +
# data[2] * (1 - x[0]) * x[1] +
# data[3] * x[0] * x[1])
bilerp = dolfin.Expression("""(d0 * (1.0 - x[0]) * (1.0 - x[1]) +
d1 * x[0] * (1.0 - x[1]) +
d2 * (1.0 - x[0]) * x[1] +
d3 * x[0] * x[1])""",
d0=data[0], d1=data[1], d2=data[2], d3=data[3],
degree=2)
cases = [(f"{N}×{N} grid, P{degree} (right diagonals)", "P",
dolfin.UnitSquareMesh(N, N)),
(f"{N}×{N} grid, P{degree} (crossed diagonals)", "P",
dolfin.UnitSquareMesh(N, N, "crossed")),
(f"{N}×{N} grid, Q{degree} (plotted as tri)", "Q",
dolfin.UnitSquareMesh.create(N, N, dolfin.CellType.Type.quadrilateral))]
if dolfin.MPI.comm_world.rank == 0:
fig, ax = plt.subplots(2, 2, constrained_layout=True, figsize=(8, 8))
ia = iter(flatten1(ax))
else:
ia = iter(range(4))
for ((title, family, mesh), a) in zip(cases, ia):
if dolfin.MPI.comm_world.rank == 0:
plt.sca(a)
V = dolfin.FunctionSpace(mesh, family, degree)
u = dolfin.interpolate(bilerp, V)
theplot = plotmagic.mpiplot(u)
if dolfin.MPI.comm_world.rank == 0:
plt.colorbar(theplot)
plt.axis("equal")
plotmagic.mpiplot_mesh(V, show_partitioning=True)
if dolfin.MPI.comm_world.rank == 0:
plt.legend(loc="upper right") # show the labels for the mesh parts
plt.title(title)
# exact
if dolfin.MPI.comm_world.rank == 0:
plt.sca(next(ia))
xx = np.linspace(0, 1, 101)
X, Y = np.meshgrid(xx, xx, indexing='xy')
zs = []
for x, y in zip(np.ravel(X), np.ravel(Y)):
zs.append(bilerp(x, y))
Z = np.reshape(zs, np.shape(X))
theplot = plt.contourf(X, Y, Z, levels=32)
plt.colorbar(theplot)
plt.axis("equal")
plt.title("Exact bilinear")
if dolfin.MPI.comm_world.rank == 0:
plt.suptitle("Interpolation using different elements")
plt.show()