circle.rar_Circle_circle matlab 网格_finite element_有限元_自适应有限元

function circle % CIRCLE simulates adaptive grids for a problem with moving circlular % singularities. % % L. Chen & C. Zhang 11-15-2006 %-------------------------------------------------------------------------- % Initialize Figure Window %-------------------------------------------------------------------------- figure(1); set(gcf,'Units','normal'); set(gcf,'Position',[0.02,0.1,0.8,0.6]); %-------------------------------------------------------------------------- % Parameters %-------------------------------------------------------------------------- theta = 0.3; theta_c = 0.1; t0 = 1.1; %-------------------------------------------------------------------------- % PreStep 0: generate initial mesh %-------------------------------------------------------------------------- node = [-1,-1; 0,-1; 1,-1; -1,0; 0,0; 1,0; -1,1; 0,1; 1,1]; elem = [2,5,1; 3,6,2; 4,1,5; 5,2,6; 5,8,4; 7,4,8; 6,9,5; 8,5,9]; Dirichlet = []; Neumann = []; mesh = getmesh(node,elem,Dirichlet,Neumann); %-------------------------------------------------------------------------- % PreStep 1: uniform mesh refinement %-------------------------------------------------------------------------- for i=1:4 mesh = bisection(mesh,ones(size(mesh.elem,1),1),1); end mesh.solu = u(mesh.node,t0); % now we pretent we know exact solution %-------------------------------------------------------------------------- % PreStep 2: mesh refinement of initial data %-------------------------------------------------------------------------- for i=1:6 eta = estimate(mesh); [mesh,eta] = bisection(mesh,eta.^2,theta); end %-------------------------------------------------------------------------- % AFEM %-------------------------------------------------------------------------- for t = t0:-0.1:0.25 % Coarsen for i = 1:5, mesh.solu = u(mesh.node,t); % now we pretent we know exact solution eta = estimate(mesh).^2; [mesh,eta] = coarsening(mesh,eta,theta_c); end % Refine for i = 1:5, mesh.solu = u(mesh.node,t); % now we pretent we know exact solution eta = estimate(mesh).^2; [mesh,eta] = bisection(mesh,eta,theta); end mesh.solu = u(mesh.node,t); % Graphic representation subplot(1,2,2); hold off; trisurf(mesh.elem,mesh.node(:,1),mesh.node(:,2),zeros(size(mesh.node,1),1)); view(2), axis equal, axis off; title(sprintf('Mesh at time %5.4f',t), 'FontSize', 14) subplot(1,2,1); hold off; trisurf(mesh.elem, mesh.node(:,1), mesh.node(:,2), mesh.solu', ... 'FaceColor', 'interp', 'EdgeColor', 'interp'); view(60,50), axis equal, axis off; title('Solution to the circle problem', 'FontSize', 14) pause(0.2) end mesh %-------------------------------------------------------------------------- % End of function CIRCLE %-------------------------------------------------------------------------- %-------------------------------------------------------------------------- % Sub functions called by CIRCLE %-------------------------------------------------------------------------- function z = u(p,t) % exact solution % % NOTE % Singularities is on a moving circle % % epsilon is the width of circle epsilon = 0.05; r = max(sum(p.^2,2).^(1/2),eps); z = exp(-((r-0.5*t)/epsilon).^2); %--------------------------------------------------------------------------