circle_area.m

来自「有限元计算中」· M 代码 · 共 20 行

exact_area =pi*25^2;
errors = []
mesh_sizes = [32, 16, 8, 4, 2];
for mesh_size = mesh_sizes
    [fens,gcells,groups,edge_gcells,edge_groups]=targe2_mesher({...
        'curve 1 circle  Center 13 4.2 radius 25',...
        'subregion 1  property 1 boundary 1',...
        ['m-ctl-point constant ' num2str(mesh_size)]
        }, 1.0);
    drawmesh({fens,gcells},'gcells','shrink', 1.0);
    view(2)
    geom = field(struct ('name',['geom'], 'dim', 2, 'fens',fens));
    feb = feblock (struct ('mater',[], 'gcells',gcells,...
        'integration_rule', tri_rule (1)));
    approximate_area =measure(feb,geom,inline('x(2)','x'));
    errors = [errors abs(exact_area-approximate_area)]
    disp([' The area is = ' num2str(approximate_area) ', to be compared with ' num2str(exact_area)])
    pause(1); close all
end
loglog(mesh_sizes, errors, 'bo-'); grid on