-
Notifications
You must be signed in to change notification settings - Fork 0
/
L2projection2D.py
84 lines (58 loc) · 1.96 KB
/
L2projection2D.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
# This is the main file.
# this code can be used to find the L2 projection of a given function in 1D.
# Code by: Maged Shaban
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from meshgen import *
from poltall import *
p,t = meshgen()
def PolyArea(x,y):
return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
# This function to evaluate the Mass Matrix
def massmat(p,t):
np1 = p.shape[1] # number of nodes
nt = t.shape[1] # number of elements
M = np.zeros((np1,np1)) # allocate the Mass Matrix
for k in range(nt):
loc2glob = t[0:3,k]-1
x = p[0,loc2glob]
y = p[1,loc2glob]
area = PolyArea(x,y)
MK = np.array(([2,1,1],[1,2,1],[1,1,2]))/12*area
#assmble the local matrix to the global
M[np.ix_(loc2glob,loc2glob)] += MK
return M
# The given function
def func(x,y):
return y + x
#np.sin(np.sqrt(x**2+ y**2))
# This laod vector function evaluation.
def LoadVec2D(p,t):
np1 = p.shape[1]
nt = t.shape[1]
b = np.zeros((np1,1))
for k in range(nt):
loc2glob = t[0:3,k]-1
x = p[0,loc2glob]
y = p[1,loc2glob]
area = PolyArea(x,y)
bk =np.array(([func(x[0],y[0])],
[func(x[1],y[1])],
[func(x[2],y[2])])) /3*area
b[loc2glob] += bk
return b
M = massmat(p,t)
b = LoadVec2D(p,t)
pf = np.linalg.inv(M)@b
#print('The Mass Matrix(M) = \n',M)
print('\nThe load vector (b)= \n',b)
print('\nThe L2 Projection (pf) of the given function f(x) is: \n',pf)
print('\n=======================================\n')
print('The total number of nodes = ',p.shape[1])
print('The total number of elements = ', t.shape[1])
print('\n=======================================\n')
#delunay_tri(M)
poltall(x = p[0,:],y = p[1,:], zt = pf)