-
Notifications
You must be signed in to change notification settings - Fork 1
/
gaussQuad2.py
34 lines (31 loc) · 1.09 KB
/
gaussQuad2.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
from gaussNodes import *
import numpy as np
def gaussQuad2(f,x,y,m):
def jac(x,y,s,t):
J = np.zeros((2,2))
J[0,0] = -(1.0 - t)*x[0] + (1.0 - t)*x[1] \
+ (1.0 + t)*x[2] - (1.0 + t)*x[3]
J[0,1] = -(1.0 - t)*y[0] + (1.0 - t)*y[1] \
+ (1.0 + t)*y[2] - (1.0 + t)*y[3]
J[1,0] = -(1.0 - s)*x[0] - (1.0 + s)*x[1] \
+ (1.0 + s)*x[2] + (1.0 - s)*x[3]
J[1,1] = -(1.0 - s)*y[0] - (1.0 + s)*y[1] \
+ (1.0 + s)*y[2] + (1.0 - s)*y[3]
return (J[0,0]*J[1,1] - J[0,1]*J[1,0])/16.0
def map(x,y,s,t):
N = np.zeros(4)
N[0] = (1.0 - s)*(1.0 - t)/4.0
N[1] = (1.0 + s)*(1.0 - t)/4.0
N[2] = (1.0 + s)*(1.0 + t)/4.0
N[3] = (1.0 - s)*(1.0 + t)/4.0
xCoord = np.dot(N,x)
yCoord = np.dot(N,y)
return xCoord,yCoord
s,A = gaussNodes(m)
sum = 0.0
for i in range(m):
for j in range(m):
xCoord,yCoord = map(x,y,s[i],s[j])
sum = sum + A[i]*A[j]*jac(x,y,s[i],s[j]) \
*f(xCoord,yCoord)
return sum