-
Notifications
You must be signed in to change notification settings - Fork 0
/
integrate.py
93 lines (84 loc) · 3.71 KB
/
integrate.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
85
86
87
88
89
90
91
92
93
import numpy as np
from scipy.special import roots_jacobi
def integral(v,t,g,q):
# =============================================================================
# # computes the integral of a function g on the mesh [v,t]
# v is an (nv by 2) matrix of the mesh vertices
# t is an (nt by 3) matrix of the mesh triangles
# output: int_{Omega} f(x,y) dxdy
# # uses Stroud conical quadrature at Gauss-Jacobi nodes
# see Ainsworth, Andriamaro, Davydov, SISC 2011
# test: g = lambda x,y: 0.*x + 1.0
# q = 5
# val = integral(v,t,g,q)
# check that val = area of the region
# by Shelvean Kapita, 2022
# =============================================================================
a = area(v,t)
[t1,w1] = roots_jacobi(q,1,0) # Gauss Jacobi nodes/weights on [-1,1]
[t2,w2] = roots_jacobi(q,0,0)
wx = 0.5*w1 # Gauss Jacobi weights on [0,1]
wy = 0.5*w2
x,y = stroudnodes(v,t,q) # Stroud nodes on the mesh
f = g(x,y) # evaluate g on the Stroud nodes
qs = q*q # number of points per triangle
w = np.outer(wy,wx).reshape(1,q,q)
ntri = int(f.shape[0]/qs) # number of triangles
nn = np.arange(f.size)
idx = np.floor(nn/qs).astype(int) # find triangle containing Stroud point
a1 = a[idx] # find area of triangle containing Stroud point
a1 = a1.reshape(a1.size,1)
ze = (a1*f).reshape(ntri,q,q)
v = np.sum(np.tensordot(ze,w,axes=([1,2],[1,2]))) # Stroud quadrature
return v
def stroudnodes(v,t,q):
# =============================================================================
# tensorized computation of Stroud nodes on the mesh
# this maps the Gauss-Jacobi points on the square [0,1]^2 to the physical triangles
# computes Stroud nodes on all triangles at once as a tensor multiplication
# v is an (nv by 2) matrix of the vertex coordinates
# t is an (nt by 3) matrix of the mesh triangles
# q is the number of Gauss-Jacobi points
# q^2 is the number of points per triangle
# output: x, y coordinates of the Stroud points on the mesh.
# =============================================================================
vt = v[t[:,:],:]
[t11,w11] = roots_jacobi(q,1,0) # Gauss Jacobi nodes on [-1,1]
[t22,w22] = roots_jacobi(q,0,0)
t1 = 0.5 + 0.5*t11 # transform Gauss-Jacobi nodes to [0,1]
t2 = 0.5 + 0.5*t22
tx,ty = np.meshgrid(t1,t2)
u = tx # barycentric coordinates at the Gauss-Jacobi nodes on square
v = ty*(1-u)
w = 1-u-v
z = np.multiply.outer(vt[:,0,:],u)+np.multiply.outer(vt[:,1,:],v)+\
np.multiply.outer(vt[:,2,:],w)
x = z[:,0,:].flatten() # corresponding points on the mesh
y = z[:,1,:].flatten()
return x.reshape(x.size,1),y.reshape(y.size,1)
def triarea(v1, v2, v3):
# =============================================================================
# # computes the signed area of a triangle with vertices v1, v2, v3
# # cc oriented triangles have positive signed area
# # clockwise oriented triangles return negative signed area
# =============================================================================
v1 = v1.reshape(1,2)
v2 = v2.reshape(1,2)
v3 = v3.reshape(1,2)
x = v1[0,0]; y = v1[0,1]
a = v2[0,0]; b = v2[0,1]
c = v3[0,0]; d = v3[0,1]
A = ((a-x)*(d-y)-(c-x)*(b-y))/2
return A
def area(v,t):
# =============================================================================
# # returns vector of areas of the triangulation
# =============================================================================
n = t.shape[0]
A = np.zeros(shape=(n,1)).flatten()
for i in np.arange(n):
v1 = v[t[i,0],:]
v2 = v[t[i,1],:]
v3 = v[t[i,2],:]
A[i] = np.abs(triarea(v1,v2,v3))
return A