gusucode.com > fem0 matlab源码程序 > fem0/fem0.m
% program fem0.m % % This solves the PDE problem (D) with homogeneous Dirichlet % conditions using piecewise linear finite elements. clear format long [nodes,elements] = initial_mesh; disp('Initial mesh is displayed in Figure Window') visualise(nodes,elements) disp('Type Return to Continue') pause nref = input('Type how many times you want to refine the initial mesh: '); disp('Refinement Level...') for i = 1:nref [nodes,elements] = refine(nodes,elements); disp(i) end [nnodes,m] = size(nodes); disp('Refined mesh is displayed in Figure Window (if not too large)') if (nnodes < 1000) visualise(nodes,elements) end disp('Type Return to Continue') pause % First assemble the element stiffness matrices elt_matrices = elt_stiffness(elements,nodes); disp('Element stiffness matrices computed') % Next assemble the global stiffness matrix % (including all boundary nodes) Ahat = global_stiffness(elt_matrices,elements,nodes); disp('Global stiffness matrix computed') N0 = find(nodes(:,3)~=1); % Finds the indices of the % interior (non-Dirichlet) nodes A = Ahat(N0,N0); % Select the appropriate blocks of Ahat % to get the matrix A from lectures % Assemble the right hand side vector (the function in the RHS of % (D) has to be specified in source_term) f = rhs('source_term',N0,elements,nodes); U = A\f; % Find coefficient vector U of finite element solution % at the degrees of freedom % Finally for the purpose of drawing graphs of the solution % and for computing errors, it is useful to assemble an extended % vector Uhat which contains the elements of U at the degrees of % freedom together with the Dirichlet boundary condition at the % boundary nodes. The ordering of the elements of Uhat coincides % with the ordering of the nodes. Uhat = zeros(nnodes,1); % Define an extended vector Uhat(N0) = U; % put in the elements of U at the degrees of freedom % Plot the solution using the MATLAB function surf % It is first necessary to interpolate Uhat on a uniform mesh. % We use the MATLAB functions meshgrid and griddata to do this. h_unif = 1/(2^(nref+1)); [XI,YI] = meshgrid([min(nodes(:,1)):h_unif:max(nodes(:,1))],... [min(nodes(:,2)):h_unif:max(nodes(:,2))]); ZI = griddata(nodes(:,1),nodes(:,2),Uhat,XI,YI,'linear'); hold on surf(XI,YI,ZI); hold off disp('Numerical solution is displayed in Figure Window') disp('Type Return to Continue') pause % Evaluate the exact solution at the nodes x1 = nodes(:,1); x2 = nodes(:,2); Uex = x1.*(1-x1).*exp(x1).*x2.*(1-x2); % Calculate error at nodes and plot it err = abs(Uhat - Uex); ZI = griddata(nodes(:,1),nodes(:,2),err,XI,YI,'linear'); surf(XI,YI,ZI); disp('Error is displayed in Figure Window') disp('Type Return to Continue') pause % Calculate a measure for the error (i.e. the energy norm of the difference % between the FE solution u_h and the interpolant of the exact solution u) anorm_err = sqrt((Uhat-Uex)'*Ahat*(Uhat-Uex))