I found this discussion in my previous search which shows how to use Python for numerical integration of a function on points in space.
Looking for Python package for numerical integration over a tessellated domain
My first question is will the code given in the above link work for an existing STL mesh? I have all the points and the surface normals for the mesh. I tried running this code and it doesn't work for me by passing the vertices and the indices of the triangles.
I also found this link in my google searches for numerical integration for FEM and found a small code snippet for integration of a function over a triangle in space as shown below:
def dJ( u, v, p1, p2, p3 ):
x1, y1 = p1
x2, y2 = p2
x3, y3 = p3
dxdu = ( (1-v)*x2 + v*x3 - x1 )
dxdv = ( u*x3 - u*x2 )
dydu = ( (1-v)*y2 + v*y3 - y1 )
dydv = ( u*y3 - u*y2 )
return np.abs( dxdu*dydv - dxdv*dydu )
def tridblquad( integrand, p1, p2, p3 ):
#http://connor-johnson.com/2014/03/09/integrals-over-arbitrary-triangular-regions-for-fem/
'''
Perform double quadtrature integration on triangular domain.
Input: function to integrate, points of triangle as tuples.
Output: integral and estimated absolute error as a tuple.
'''
x1, y1 = p1 ; x2, y2 = p2 ; x3, y3 = p3
# transformation to the unit square
g = lambda u, v, c1, c2, c3: (1-u)*c1 + u*( (1-v)*c2 + v*c3 )
# transformation for the integrand,
# including the Jacobian scaling factor
def h( u, v ):
x = g( u, v, x1, x2, x3 )
y = g( u, v, y1, y2, y3 )
I = integrand( x, y )
I *= dJ( u, v, p1, p2, p3 )
return I
# perfrom the double integration using quadrature in the transformed space
integral, error = scipy.integrate.dblquad( h, 0, 1, lambda x: 0, lambda x: 1, epsrel=1e-6, epsabs=0 )
return integral, error
However, this code works only for 2 dimensional triangles. Can anyone suggest a better more general way for performing numerical integration of a function over a tessellated mesh?
Thanks.
