InΒ [1]:
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
source module for methodsNM imported.

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}$.ΒΆ

InΒ [2]:
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()
Warning: Using NumPy Gauss rules!
Warning: Using numerical differentiation for deriv evaluation in <class 'methodsnm.fe_1d.Lagrange_Segment_FE'> object.
No description has been provided for this image
No description has been provided for this image

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).

InΒ [3]:
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()
No description has been provided for this image

"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$

InΒ [4]:
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()
No description has been provided for this image

Refining the mesh by adding an extra point in the middle of the elements we want to refine

InΒ [5]:
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()
No description has been provided for this image
InΒ [6]:
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()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

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)

&quot;&quot;

  • compute_refine funktioniert gleich
  • refine mesh?

&quot;&quot;