import import_hack
from methodsnm.mesh_1d import *
from methodsnm.visualize import *
from methodsnm.fes import *
from methodsnm.forms import *
from methodsnm.formint import *
import matplotlib.pyplot as plt
from methodsnm.poisson_1d import *
try: #include solution module (if exists)
from methodsnm.solution import *
except:
pass
Adaptive mesh refinement for the Poisson problem [ Paula John + Marvin Langer ]ΒΆ
Implement an adaptive mesh refinement strategy for the Poisson problem.
All functions used for mesh refinement are located in poisson_1d.py.
PrerequisitesΒΆ
To solve the problem below, the second derivative was needed. For the P1_Segment_FE
and P2_Segment_FE
we implemented the derivative by hand and for the rest we implemented a recursive function to calculate the numerical derivative for any order n. See fe_1d.py.
Problem: $- \alpha u'' + \gamma u = f$ with $\gamma = 1$, $\alpha \ll 1$ and $f = \chi_{[0.4,0.6]}$ (i.e. a bump in the middle of the domain), $\Omega = (0,1)$.ΒΆ
Computing the numerical solution for $\alpha \in \{ 10^{-i} \}_{i=1,2,3,4,5}$.ΒΆ
ne = 30
mesh = Mesh1D((0,1),ne)
fes = Lagrange_Segments_Space(mesh, order=3)
f = GlobalFunction(lambda x: 1 if 0.4<=x and x<= 0.6 else 0, mesh = mesh)
for i in range(1,6):
alpha = 10**(-1*i)
gamma = 1
uh = solve_poi_1D(mesh,fes, alpha, gamma, f)
DrawFunction1D(uh, label = "i = %s"% i )
plt.legend()
plt.show()
DrawFunction1D(f, label = "f" )
plt.legend()
plt.show()
Posteriori error estimation for each element $T \in \mathcal{T}_h$:ΒΆ
$$ \eta_T = h_T^2 \| - \alpha u_h'' + \gamma u_h - f \|_{L^2(T)}^2 + \sum_{F \subset \partial T} h_T \| \llbracket u_h' \rrbracket \|_{L^2(F)}^2 $$
Here $h_T$ denotes the diameter of the element $T$ and $\llbracket u' \rrbracket$ denotes the jump of the derivative of $u$ across the face $F$ (vertex in 1D case).
ne = 30
for i in range(1,6):
alpha = 10**(-1*i)
gamma = 1
mesh = Mesh1D((0,1),ne)
fes = P1_Segments_Space(mesh)
f = GlobalFunction(lambda x: 1 if 0.4<=x and x<= 0.6 else 0, mesh = mesh)
uh = solve_poi_1D(mesh,fes, alpha, gamma, f)
mesh_error = compute_mesh_error(uh, alpha, gamma, f)
piecew = lambda x: mesh_error[int(x)]
xs = np.linspace(0,len(mesh_error)-1,1000)
ys = [piecew(x) for x in xs]
plt.plot(xs,ys,label = "i = %s"% i)
plt.legend()
plt.show()
"DΓΆrfler" marking: selecting the elements with largest $\eta_T$ that add up to $20 \%$ of the total error indicator $\sum_T \eta_T$ΒΆ
example: $i=5$
ne =30
i = 5
alpha = 10**(-1*i)
gamma = 1
f = GlobalFunction(lambda x: 1 if 0.4<=x and x<= 0.6 else 0, mesh = mesh)
mesh = Mesh1D((0,1),ne)
fes = Lagrange_Segments_Space(mesh, order=3)
uh = solve_poi_1D(mesh,fes, alpha, gamma, f)
#what to refine
mesh_error = compute_mesh_error(uh, alpha, gamma, f)
refine = compute_refine(mesh_error)
xs = np.arange(0,len(mesh_error),1)
points = mesh.points[:-1]
widths = np.array([(mesh.points[p+1]-mesh.points[p])*0.9 for p in range(len(points))])
plt.bar(points[refine],mesh_error[refine],align="edge",width = widths[refine], color = "orange",alpha=0.8, label = "refine")
plt.bar(points[np.invert(refine)],mesh_error[np.invert(refine)],align="edge",width = widths[np.invert(refine)], color = "blue",alpha=0.6, label = "ok")
plt.legend()
plt.show()
Refining the mesh by adding an extra point in the middle of the elements we want to refine
mesh = refine_mesh(mesh, refine)
fes = Lagrange_Segments_Space(mesh, order=3)
uh = solve_poi_1D(mesh,fes, alpha, gamma, f)
mesh_error = compute_mesh_error(uh, alpha, gamma, f)
refine = compute_refine(mesh_error)
xs = np.arange(0,len(mesh_error),1)
points = mesh.points[:-1]
widths = np.array([(mesh.points[p+1]-mesh.points[p])*0.9 for p in range(len(points))])
plt.bar(points[refine],mesh_error[refine],align="edge",width = widths[refine], color = "orange",alpha=0.8, label = "refine")
plt.bar(points[np.invert(refine)],mesh_error[np.invert(refine)],align="edge",width = widths[np.invert(refine)], color = "blue",alpha=0.6, label = "ok")
plt.legend()
plt.show()
ne =30
i = 5
tol = 0.01
alpha = 10**(-1*i)
gamma = 1
f = GlobalFunction(lambda x: 1 if 0.4<=x and x<= 0.6 else 0, mesh = mesh)
mesh = Mesh1D((0,1),ne)
fes = Lagrange_Segments_Space(mesh, order=3)
uh = solve_poi_1D(mesh,fes, alpha, gamma, f)
mesh_error = compute_mesh_error(uh, alpha, gamma, f)
n = 0
while sum(mesh_error)> tol:
refine = compute_refine(mesh_error)
xs = np.arange(0,len(mesh_error),1)
DrawFunction1D(uh, label = "refinements: %s"%n, show_mesh=True)
plt.show()
points = mesh.points[:-1]
widths = np.array([(mesh.points[p+1]-mesh.points[p])*0.9 for p in range(len(points))])
plt.bar(points[refine],mesh_error[refine],align="edge",width = widths[refine], color = "orange", alpha=0.8, label = "refine")
plt.bar(points[np.invert(refine)],mesh_error[np.invert(refine)],align="edge",width = widths[np.invert(refine)], color = "blue",alpha=0.6, label = "ok")
plt.legend()
plt.yscale("log")
plt.show()
mesh = refine_mesh(mesh, refine)
f = GlobalFunction(lambda x: 1 if 0.4<=x and x<= 0.6 else 0, mesh = mesh)
fes = Lagrange_Segments_Space(mesh, order=3)
uh = solve_poi_1D(mesh,fes, alpha, gamma, f)
n+=1
mesh_error = compute_mesh_error(uh, alpha, gamma, f)
DrawFunction1D(uh, label = "refinements: %s"%n, show_mesh=True)
plt.show()
DrawFunction1D(uh, label = "refinements: %s"%n, show_mesh=True)
plt.yscale("log")
plt.show()
Ausblick 2DΒΆ
sinnvolle Wahl von $h_T$
Ableitungen implementieren
Randintegrale (1D)
- durchgehen der Edges
- Edge muss benachbarte finite Elemente kennen (analog zu faces2edges in StructuredRectangleMesh)
- Orientierung beachten (gleiches Problem wie bei P3_Triangle_Space)
- compute_refine funktioniert gleich
- refine mesh?