Element Matrizen 2d Fall

Inhalt

9.2.4. Element Matrizen 2d Fall#

Die Wahl der Basisfunktionen für finite Elemente Räume ist ein grosses Gebiet in der Mathematik. In NGSolve ist für die Approximation des \(H^1\) eine \(hp\) finite Elemente Methode implementiert. Wir betrachten die FE-Basis Funktionen auf dem Einheitsdreieck und Einheitsquadrat.

from ngsolve import *
from ngsolve.webgui import Draw
import numpy as np
from pandas import DataFrame

def highlight_ortho(s):
    is_ortho = np.abs(s) < 1e-13
    return ['background-color: yellow' if v else '' for v in is_ortho]

Wir betrachten die in NGSolve implementierten Basisfunktionen auf verschiedenen Elemente:

# Einheits Dreieck (0,0),(1,0),(0,1)
def Trig():
    from netgen.meshing import Mesh, MeshPoint, Element2D, Pnt
    mesh = Mesh(2)
    mesh.AddRegion("",2)
    pnts = [mesh.Add(MeshPoint(Pnt(0,0,0))), mesh.Add(MeshPoint(Pnt(1,0,0))), mesh.Add(MeshPoint(Pnt(0,1,0)))]
    mesh.Add(Element2D(1, pnts))
    return mesh

# Einheits Dreieck (-1,-1),(1,-1),(0,1)
def Trig2():
    from netgen.meshing import Mesh, MeshPoint, Element2D, Pnt
    mesh = Mesh(2)
    mesh.AddRegion("",2)
    pnts = [mesh.Add(MeshPoint(Pnt(-1,-1,0))), mesh.Add(MeshPoint(Pnt(1,-1,0))), mesh.Add(MeshPoint(Pnt(0,1,0)))]
    mesh.Add(Element2D(1, pnts))
    return mesh

# Dreieck (1,1),(1.5,-1),(2,1.2)
def Trig3():
    from netgen.meshing import Mesh, MeshPoint, Element2D, Pnt
    mesh = Mesh(2)
    mesh.AddRegion("",2)
    pnts = [mesh.Add(MeshPoint(Pnt(1,1,0))),
            mesh.Add(MeshPoint(Pnt(1.5,-1,0))),
            mesh.Add(MeshPoint(Pnt(2,1.2,0)))]
    mesh.Add(Element2D(1, pnts))
    return mesh

# Einheits Quadrat
def Quad():
    from netgen.meshing import Mesh, MeshPoint, Element2D, Pnt
    mesh = Mesh(2)
    mesh.AddRegion("",2)
    pnts = [mesh.Add(MeshPoint(Pnt(0,0,0))), mesh.Add(MeshPoint(Pnt(1,0,0))), mesh.Add(MeshPoint(Pnt(1,1,0))), mesh.Add(MeshPoint(Pnt(0,1,0)))]
    mesh.Add(Element2D(1, pnts))
    return mesh

Einheitsdreieck#

mesh = Mesh(Trig())
Draw(mesh);

Lineare Ordnung#

V = H1(mesh,order=1)
gfu = GridFunction(V)

print('Anzahl Basis Funktionen: ',V.ndof)

Anzahl Basis Funktionen:  3

Die Basisfunktionen auf dem Einheitsdreick sind gegeben durch

\[\{1-x-y, x, y\}\]

gfu.vec.FV()[:] = 0
gfu.vec.FV()[1] = 1
Draw(gfu,mesh);

gfu1 = GridFunction(V)
gfu2 = GridFunction(V)
ortho = []
for i in range(0,V.ndof):
    gfu1.vec.FV()[:] = 0
    gfu1.vec.FV()[i] = 1
    oi = []
    for j in range(0,V.ndof):
        gfu2.vec.FV()[:] = 0
        gfu2.vec.FV()[j] = 1
        oi.append(Integrate(grad(gfu1)*grad(gfu2),mesh))
    ortho.append(oi)
ortho = DataFrame(ortho)
ortho.style.\
    apply(highlight_ortho).\
    set_table_attributes('style="font-size: 12px"')

	0	1	2
0	1.000000	-0.500000	-0.500000
1	-0.500000	0.500000	0.000000
2	-0.500000	0.000000	0.500000

gfu1 = GridFunction(V)
gfu2 = GridFunction(V)
ortho = []
for i in range(0,V.ndof):
    gfu1.vec.FV()[:] = 0
    gfu1.vec.FV()[i] = 1
    oi = []
    for j in range(0,V.ndof):
        gfu2.vec.FV()[:] = 0
        gfu2.vec.FV()[j] = 1
        oi.append(Integrate(gfu1*gfu2,mesh))
    ortho.append(oi)
ortho = DataFrame(ortho)
ortho.style.\
    apply(highlight_ortho).\
    set_table_attributes('style="font-size: 12px"')

	0	1	2
0	0.083333	0.041667	0.041667
1	0.041667	0.083333	0.041667
2	0.041667	0.041667	0.083333

High Order FEM#

V = H1(mesh,order=2)
gfu = GridFunction(V)

print('Anzahl Basis Funktionen: ',V.ndof)

Anzahl Basis Funktionen:  6

Hier wird es bezüglich Basisfunktionen schon komplexer. Die ersten drei Basisfunktionen sind die linearen. Danach kommen die drei quadratische Funktionen:

gfu.vec.FV()[:] = 0
gfu.vec.FV()[4] = 1
s=Draw(gfu,mesh)

gfu.vec.FV()[:] = 0
gfu.vec.FV()[3] = 1
s.Redraw()

Gewisse Basisfunktionen sind zueinander orthogonal bezüglich dem Skalarprodukt \((\nabla u,\nabla v)\):

gfu1 = GridFunction(V)
gfu2 = GridFunction(V)
ortho = []
for i in range(0,V.ndof):
    gfu1.vec.FV()[:] = 0
    gfu1.vec.FV()[i] = 1
    oi = []
    for j in range(0,V.ndof):
        gfu2.vec.FV()[:] = 0
        gfu2.vec.FV()[j] = 1
        oi.append(Integrate(grad(gfu1)*grad(gfu2),mesh))
    ortho.append(oi)
ortho = DataFrame(ortho)
ortho.style.\
    apply(highlight_ortho).\
    set_table_attributes('style="font-size: 10px"')

	0	1	2	3	4	5
0	1.000000	-0.500000	-0.500000	-0.083333	-0.083333	0.166667
1	-0.500000	0.500000	0.000000	-0.000000	0.083333	-0.083333
2	-0.500000	0.000000	0.500000	0.083333	0.000000	-0.083333
3	-0.083333	-0.000000	0.083333	0.041667	0.000000	-0.020833
4	-0.083333	0.083333	0.000000	0.000000	0.041667	-0.020833
5	0.166667	-0.083333	-0.083333	-0.020833	-0.020833	0.041667

Konstruktion von Basisfunktionen#

from sympy import symbols, solve
from sympy.abc import x,y

Wir berechnen Basisfunktionen mit dem Ansatz der Lagrange Polynome aus dem 1d. Dazu definieren wir ein allgemeines Polynom mit einer bestimmten maximalen Ordnung.

\[p(x,y) = \sum_{i,j}^{i+j\le n} c_k x^i y^j\]

def p(x,y,c,n):
    z = 0
    k = 0
    for i in range(n+1):
        for j in range(n+1):
            if i+j<=n:
                z += c[k]*x**i*y**j
                k += 1
    return z

Für den Ansatz erster Ordnung fordern wir analog zu den Lagrange Polynome im 1d, dass die Basisfunktionen jeweilen an einer Ecke den Funktionswert 1 annehmen und an den andern 0.

# Koordinaten
pi = [[0,0],[0,1],[1,0]]
# Funktionswerte an den Koordinaten
yi = np.eye(len(pi))
# Polynom Koeffizienten
c = symbols('c:3')

Wir berechnen die zugehörigen Polynom Koeffizienten

sol=[solve(np.array([p(*pii,c,1) for pii in pi])-yi[:,i]) for i in range(yi.shape[0])]

und erhalten damit einen möglichen Ansatz für eine Basis mit der vorgeschriebenen Ordnung:

for i in range(yi.shape[0]):
    display(p(x,y,[sol[i][ci] for ci in c],1))

\[\displaystyle - 1.0 x - 1.0 y + 1.0\]

\[\displaystyle 1.0 y\]

\[\displaystyle 1.0 x\]

In dem Fall entspricht der Ansatz auch den in NGSolve implementierten.

Analog betrachten wir nun den Ansatz für Polynoe zweiter Ordnung. Dazu führen wir zusätzlich Knoten jeweilen in der Mitte jeder Kanten ein. Aus dem Pascalschen Dreieck sehen wir, dass das Polynom für beliebige Ordnung \(p\)

\[\frac{(p+1)\,(p+2)}{2}\]

Koeffizienten hat. Für die Ordnung \(p=2\) erhalten wir 6 Koeffizienten:

order=2
(order+1)*(order+2)/2

6.0

# Knoten
pi = [[0,0],[0,1],[1,0],[0.5,0],[1/2,1/2],[0,.5]]
# Funktionswerte
yi = np.eye(len(pi))
# Koeffizienten
c = symbols('c:'+str(int((order+1)*(order+2)/2)))

Berechnen wir wiederum die zugehörigen Polynom Koeffizienten, so erhalten wir den folgenden Ansatz:

sol=[solve(np.array([p(*pii,c,2) for pii in pi])-yi[:,i]) for i in range(yi.shape[0])]
phis=[p(x,y,[sol[i][ci] for ci in c],2) for i in range(yi.shape[0])]
for phii in phis:
    display(phii)

\[\displaystyle 2.0 x^{2} + 4.0 x y - 3.0 x + 2.0 y^{2} - 3.0 y + 1.0\]

\[\displaystyle 2.0 y^{2} - 1.0 y\]

\[\displaystyle 2.0 x^{2} - 1.0 x\]

\[\displaystyle - 4.0 x^{2} - 4.0 x y + 4.0 x\]

\[\displaystyle 4.0 x y\]

\[\displaystyle - 4.0 x y - 4.0 y^{2} + 4.0 y\]

Wie zu sehen ist, erhalten wir jetzt im Gegensatz zu den in NGSolve implementierten Basisfunktionen sechs Polynome zweiter Ordnung. Es stellt sich natürlich die Frage, ob dieser Ansatz geeignet ist. Per Konstruktion sind die Polynome linearunabhängig. Die in NGSolve gewählten high Order Basisfunktionen führen zu einer exponentiellen Konvergenz Ordnung unter Anwendung von \(hp\) finite Elemente Methoden.

Einheitsquadrat#

mesh = Mesh(Quad())
Draw(mesh);

Lineare Ordnung#

V = H1(mesh,order=1)
gfu = GridFunction(V)

print('Anzahl Basis Funktionen: ',V.ndof)

Anzahl Basis Funktionen:  4

gfu.vec.FV()[:] = 0
gfu.vec.FV()[0] = 1
s2=Draw(gfu,mesh);

gfu.vec.FV()[:] = 0
gfu.vec.FV()[1] = 1
s2.Redraw()

Quadratische Ordnung#

V = H1(mesh,order=2)
gfu = GridFunction(V)

print('Anzahl Basis Funktionen: ',V.ndof)

Anzahl Basis Funktionen:  9

Hier wird es bezüglich Basisfunktionen schon komplexer. Die ersten drei Basisfunktionen sind die linearen. Danach kommen die drei quadratische Funktionen:

gfu.vec.FV()[:] = 0
gfu.vec.FV()[4] = 1
s=Draw(gfu,mesh);

gfu.vec.FV()[:] = 0
gfu.vec.FV()[8] = 1
s.Redraw()

Die Basisfunktionen sind orthogonal bezüglich dem Skalarprodukt \((\nabla u,\nabla v)\)

gfu1 = GridFunction(V)
gfu2 = GridFunction(V)

ortho = []
for i in range(0,V.ndof):
    gfu1.vec.FV()[:] = 0
    gfu1.vec.FV()[i] = 1
    oi = []
    for j in range(0,V.ndof):
        gfu2.vec.FV()[:] = 0
        gfu2.vec.FV()[j] = 1

        oi.append(Integrate(grad(gfu1)*grad(gfu2),mesh))
    ortho.append(oi)
ortho = DataFrame(ortho)
ortho.style.\
    apply(highlight_ortho).\
    set_table_attributes('style="font-size: 10px"')

	0	1	2	3	4	5	6	7	8
0	0.666667	-0.166667	-0.333333	-0.166667	-0.041667	-0.041667	0.041667	0.041667	-0.000000
1	-0.166667	0.666667	-0.166667	-0.333333	-0.041667	0.041667	-0.041667	0.041667	0.000000
2	-0.333333	-0.166667	0.666667	-0.166667	0.041667	0.041667	-0.041667	-0.041667	0.000000
3	-0.166667	-0.333333	-0.166667	0.666667	0.041667	-0.041667	0.041667	-0.041667	0.000000
4	-0.041667	-0.041667	0.041667	0.041667	0.036111	-0.000000	0.000000	0.005556	-0.003472
5	-0.041667	0.041667	0.041667	-0.041667	-0.000000	0.036111	0.005556	0.000000	-0.003472
6	0.041667	-0.041667	-0.041667	0.041667	0.000000	0.005556	0.036111	0.000000	-0.003472
7	0.041667	0.041667	-0.041667	-0.041667	0.005556	0.000000	0.000000	0.036111	-0.003472
8	-0.000000	0.000000	0.000000	0.000000	-0.003472	-0.003472	-0.003472	-0.003472	0.001389