✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
基于Matlab的非线性能量的快速 FEM 评估
⛄ 部分代码
clear variables
close all
clc
add_paths;
%=========================================================================
all.levels = 1:5; % mesh refinement levels
all.p = 1:1; % polynomial degrees
all.methods_on = [1 1 0 0 0]; % options
%=========================================================================
all.RUN = kron(all.methods_on,ones(numel(all.levels),1));
all.RUN(5:end,4) = 0; % exclude some evaluations for option 4
all.RUN(5:end,5) = 0; % exclude some evaluations for option 5
%=========================================================================
elasticity_setup_parameters;
elasticity_initialize_2D;
%=========================================================================
for lev=1:length(all.levels)
level = all.levels(lev);
for pp=1:length(all.p)
p = all.p(pp);
params.p = p;
params.solver.Gauss_order = 2*params.p; % Gauss quadrature order
elasticity_setup_domain; % n2c + e2n
mesh = elasticity_setup_mesh(n2c,e2n,params);
if params.graphs.Hess
figure; spyown(mesh.Hstr,[mesh.nn mesh.nn]);
end
ndofs = length(mesh.dofsMinim); % number of degrees of freedom
all.ndofs(lev,pp) = ndofs;
all.ne(lev,pp) = mesh.ne;
u_init = zeros(2*mesh.nb,1);
elasticity_setup_u_init;
if params.graphs.mesh
figure(1);
hold on
nDfull = mesh.nodesDirichlet{1};
plot(u_init(2*nDfull-1),u_init(2*nDfull),'ro','MarkerSize',8,'Linewidth',1.5);
hold off
end
fprintf(separate1); fprintf('\n');
for m=1:5
if all.methods_on(m) && all.RUN(lev,m)
fprintf(separate1); fprintf('\n\n');
switch m
case 1
fprintf('ndofs = %d, OPTION 1 : exact gradient, known hessian pattern \n',ndofs);
case 2
fprintf('ndofs = %d, OPTION 2 : approx. gradient, known hessian pattern \n',ndofs);
case 3
fprintf('ndofs = %d, OPTION 3 : approx. gradient, no hessian pattern \n',ndofs);
case 4
fprintf('ndofs = %d, OPTION 4 : no gradient, no hessian pattern \n',ndofs);
case 5
fprintf('ndofs = %d, OPTION 5 : approx. gradient, no hessian pattern \n',ndofs);
end
[u, Ju, time, iters, densities] = approx_solution(mesh,u_init,circshift([0 0 0 0 1],m),params);
all.u{lev,pp,m} = u;
all.Ju(lev,pp,m) = Ju/1e7;
all.times(lev,pp,m) = time;
all.iters(lev,pp,m) = iters;
if ~isempty(params.graphs.solution) && ismember(1,(params.graphs.solution(:,1)==level) ...
& (params.graphs.solution(:,2)==p) & (params.graphs.solution(:,3)==m))
visualize_elasticity_2D;
caxis([0 7.5e6]);
end
end
end
end
end
fprintf('\n'); generate_table_paper_latex_article(all); % main article table
