基于Matlab的非线性能量的快速 FEM 评估

✅作者简介：热爱科研的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

⛄ 运行结果

⛄ 参考文献

​

❤️部分理论引用网络文献，若有侵权联系博主删除

❤️ 关注我领取海量matlab电子书和数学建模资料