Integrate Results of 3D FEM on mesh with arbitrary integration limits

2 次查看（过去 30 天）
Ryan Woodall
Ryan Woodall 2018-5-17
编辑： Ryan Woodall 2018-5-17
I have results of a 3D finite element analysis stored as a list of points [x,y,z], and corresponding values V. I would like to perform 3D integration on small sub-regions of the mesh, using a piece-wise linear interpolation of this data.
I would like to be able to specify a cube in which to integrate the interpolated solution. If a portion of the cube lies outside the mesh, I would like the interpolator to return 0, instead of NaN.
I have tried the below solution, but it throws a warning (see below) and returns 0 when the integration region hangs outside the mesh. It does work on internal regions, but it is very slow (~15 seconds per integration, n_nodes ~40,000), and n_pts is ~1000. This will eventually be wrapped in lsqnonlin() for parameter estimation, so it needs to run relatively quickly for multiple iterations.
function v = RemoveNans(Si,X,Y,Z)
v = Si(X,Y,Z);
v(isnan(v)) = 0;
end
nodes = csvread('nodes.csv'); %3D location of points, irregularly spaced, size = [n_nodes,3]
vals = csvread('values.csv'); %solution to FEM at above nodes, size = [n_nodes,1]
int_pts = csvread('int_points.csv'); %list of points to integrate over, size = [n_pts,3], n_pts << n_nodes
R = 1; %size of cube for integration
Si = scatteredInterpolant(nodes, vals, 'linear', 'none')
f = @(X,Y,Z)RemoveNans(Si,X,Y,Z);
integrated_values = zeros(size(int_pts));
for i=1:length(integrated_values)
integrated_values(i) = integral3(f,...
int_pts(i,1)-R/2,int_pts(i,1)+R/2,...
int_pts(i,2)-R/2,int_pts(i,2)+R/2,...
int_pts(i,3)-R/2,int_pts(i,3)+R/2);
end
"Warning: Reached the maximum number of function evaluations (10000). The result passes the global error test."
Any help on getting this to work properly, or a "better" way to do this would be very helpful.

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by