%---------------- function blatt5_2D %---------------- % exact_energy = 1/45; mmax = 5; energy = zeros(mmax,1); h = zeros(mmax,1); for mm = 1:mmax mm m = 2^mm; h(mm) = 1.0/m; energy(mm) = fem(m); end; %compute extrapolated energies for mm=3:mmax %BESTIMME extrapolierte Energien mittels Aitken Deltaverfahren %aus energy(mm), energy(mm-1), energy(mm-2) end figure(1) loglog(h,exact_energy - energy, '*-',... h,0.1*h.^2,... h(3:mmax),abs(exact_energy-extrapolated_energy(3:mmax)),... 'LineWidth',2) set(gca,'FontSize',18) xlabel('Gitterweite h') ylabel('Fehler in der Energie') legend('Fehler','O(h^2)','Fehler der Extrapolation') title('FEM auf Einheitsquadrat; gleichmaessiges Dreiecksgitter, glatte Loesung') %axis([1.e-3 1 1*10^(-8) 5*10^(-2)]) return end %--------------------------------------------------------- function [energy] = fem(m) % %define the geometry, i.e., the % list of nodes % list of elements % list of edges (with boundary flags) % % define geometry % [nodes,elements,edges ] =geometry_uniform_mesh(m); %[nodes,elements,edges ] =geometry_square; %---- dirichlet_nodes = determine_dirichlet_nodes(edges); N = length(nodes(:,1)); %we implicitly assume that the node list starts with 1 [B,l] = assemble_neumann_problem(nodes,elements) ; [BD,lD] = apply_Dirichlet_bc(B,l,dirichlet_nodes) ; % %COMPUTE FEM SOLUTION AND FEM ENERGY % return; end % %--------------------------------------------------------- %--------------------------------------------------------- %subroutines for FEM %--------------------------------------------------------- %--------------------------------------------------------- function [BD,lD] = apply_Dirichlet_bc(B,l,dirichlet_nodes); %simply delete the rows and columns of B that correspond to Dirichlet nodes %simply delete the row of l that correspond to Dirichlet nodes N = size(B,1); active_nodes = determine_active_nodes(N,dirichlet_nodes); BD = B(active_nodes,active_nodes); lD = l(active_nodes); return end %--------------------------------------------------------- function [B,l] = assemble_neumann_problem(nodes,elements) %assembles the *unconstrained* stiffness matrix, i.e., the fact that Dirichlet nodes %may exist is ignored N = length(nodes(:,1)); B = sparse(N,N); % d.h. B = zeros(N,N); l = zeros(N,1); % %ASSEMBLE STIFFNESS MATRIX AND LOAD VECTOR % return end % % % % %--------------------------------------------------------- function BK = element_stiffness_matrix(vertex_coords); %input: vertex_coords(2:3) contains the coordinates of the three vertices of the triangle BK = zeros(3,3); v1 = vertex_coords(:,1) ; v2 = vertex_coords(:,2) ; v3 = vertex_coords(:,3) ; % %COMPUTE ELEMENT STIFFNESS MATRIX % return end %--------------------------------------------------------- function lK = element_load_vector(vertex_coords); %input: vertex_coords(2:3) contains the coordinates of the three vertices of the triangle BK = zeros(3,1); v1 = vertex_coords(:,1); v2 = vertex_coords(:,2); v3 = vertex_coords(:,3); % %COMPUTE ELEMENT LOAD VECTOR % return end %--------------------------------------------------------- function dirichlet_nodes = determine_dirichlet_nodes(edges); %determine Dirichlet nodes, i.e., those that are associated with DoF on the Dirichlet part of the bdy % %find all edges that are Dirichlet edges Dedges = find(edges(:,4)); %find spits out indices of non-zero entries if length(Dedges) == 0 dirichlet_nodes = []; return; end; % Dnodes = [edges(Dedges,2); edges(Dedges,3)]; %grab all nodes that are endpoints of Dirichlet edges % remove possible duplicates n = max(Dnodes); %get highest node number tmp = zeros(n,1); tmp(Dnodes) = 1; %tmp has non-zeros entries only at indices corresponding %to Dirichlet nodes dirichlet_nodes = find(tmp); %find gives the indices of the nonzero entries return end %--------------------------------------------------------- function active_nodes = determine_active_nodes(N,dirichlet_nodes) %returns a list of nodes that are the active nodes, i.e., those that are not on the %Dirichlet part of the boundary %N = total number of nodes %dirichlet_nodes = list of nodes that are on the Dirichlet part % if length(dirichlet_nodes) == 0 active_nodes = [1:N]; return end; tmp = ones(N,1); tmp(dirichlet_nodes) = 0; %tmp(i) = 0 iff i is Dirichlet node active_nodes = find(tmp); % find returns indices with nonzero entries return end %--------------------------------------------------------- function y = N1(x) xi = x(1); eta = x(2); y = 1 - xi - eta; return end %--------------------------------------------------------- function g = grad_N1(x) xi = x(1); eta = x(2); g = [-1; -1]; return end %--------------------------------------------------------- %--------------------------------------------------------- function y = N2(x) xi = x(1); eta = x(2); y = xi; return end %--------------------------------------------------------- function g = grad_N2(x) xi = x(1); eta = x(2); g = [1; 0]; return end %--------------------------------------------------------- %--------------------------------------------------------- function y = N3(x) xi = x(1); eta = x(2); y = eta; return end %--------------------------------------------------------- function g = grad_N3(x) xi = x(1); eta = x(2); g = [0; 1]; return end %--------------------------------------------------------- function y = f(x) %the right-hand side function x1 = x(1); x2 = x(2); y = 2*(x1*(1-x1)+x2*(1-x2)); %this corresponds to the exact solution u(x,y) = x*(1-x)*y*(1-y) with %energy on unit square: 1/45 return end % %---------------------------- function [nodes,elements,edges] = geometry_triangle %---------------------------- %simplest geometry: single triangle %all edges are Neumann edges % nodes = [ 1 0.0 0.0 2 1.0 0.0 3 1.0 1.0]; elements = [ 1 1 2 3 ]; edges = [ 1 1 2 0 2 2 3 0 3 1 3 0 ]; return end %---------------------------- function [nodes,elements,edges ] = geometry_square %---------------------------- %simple geometry: square split into 2 triangles %all edges are Neumann edges % nodes = [ 1 0.0 0.0 2 1.0 0.0 3 1.0 1.0 4 0.0 1.0]; elements = [ 1 1 2 3 2 1 3 4 ]; edges = [ 1 1 2 0 2 2 3 0 3 1 3 0 4 1 4 0 5 4 3 0 ]; return end %--------------------------------------------------------- function [nodes,elements,edges] = geometry_uniform_mesh(n) % %generates a uniform mesh on [0,1]^2 with (n+1)^2 nodes %we prescribe Dirichlet bc on the full boundary % h = 1.0/n; nnodes = (n+1)*(n+1); nele = 2*n*n; nedges = 3*n*n+n+n; nodes = zeros(nnodes,3); elements = zeros(nele,4); edges = zeros(nedges,4); %define nodes counter = 1; for i=0:n for j=0:n nodes(counter,:) = [counter j*h i*h]; counter = counter+1; end; end; %define elements counter = 1; for i=0:n-1 for j=0:n-1 v1 = i*(n+1)+j+1 ; v2 = i*(n+1)+j+1+1 ; v3 = (i+1)*(n+1)+j+1 ; v4 = (i+1)*(n+1)+j+1+1 ; elements(counter,:) = [counter v1 v2 v4]; counter = counter+1; elements(counter,:) = [counter v1 v4 v3]; counter = counter+1; end; end; % %define edges counter = 1; flag2 = 0; for i=0:n-1 if (i==0) flag1 = 1; else flag1 = 0; end; % for j=0:n-1 if (j==0) flag3 = 1; else flag3 = 0; end; % v1 = i*(n+1)+j+1 ; v2 = i*(n+1)+j+1+1 ; v3 = (i+1)*(n+1)+j+1 ; v4 = (i+1)*(n+1)+j+1+1 ; edges(counter,:) = [counter v1 v2 flag1]; counter = counter+1; edges(counter,:) = [counter v1 v4 flag2]; counter = counter+1; edges(counter,:) = [counter v1 v3 flag3]; counter = counter+1; end; edges(counter,:) = [counter v2 v4 1]; counter = counter+1; end; i = n-1; for j=0:n-1 v1 = i*(n+1)+j+1 ; v2 = i*(n+1)+j+1+1 ; v3 = (i+1)*(n+1)+j+1 ; v4 = (i+1)*(n+1)+j+1+1 ; edges(counter,:) = [counter v3 v4 1]; counter = counter+1; end; return end % %---------------------------- function [nodes,elements,edges ] = geometry_A %---------------------------- % the letter 'A' as produced by the code "triangle" % one edge is Dirichlet edge %---------------------------- % number x y nodes = [ 1 0.20000000000000001 -0.77639999999999998 2 0.22 -0.7732 3 0.24560000000000001 -0.75639999999999996 4 0.27760000000000001 -0.70199999999999996 5 0.48880000000000001 -0.20760000000000001 6 0.50480000000000003 -0.20760000000000001 7 0.74080000000000001 -0.73960000000000004 8 0.75600000000000001 -0.76119999999999999 9 0.77439999999999998 -0.77239999999999998 10 0.80000000000000004 -0.77639999999999998 11 0.80000000000000004 -0.79239999999999999 12 0.57920000000000005 -0.79239999999999999 13 0.57920000000000005 -0.77639999999999998 14 0.62160000000000004 -0.77159999999999995 15 0.63360000000000005 -0.76280000000000003 16 0.63919999999999999 -0.74439999999999995 17 0.62080000000000002 -0.68440000000000001 18 0.58720000000000006 -0.60440000000000005 19 0.36080000000000001 -0.60440000000000005 20 0.31919999999999998 -0.70679999999999998 21 0.312 -0.73960000000000004 22 0.31840000000000002 -0.76119999999999999 23 0.33439999999999998 -0.77159999999999995 24 0.37119999999999997 -0.77639999999999998 25 0.37119999999999997 -0.79239999999999999 26 0.37440000000000001 -0.56999999999999995 27 0.57440000000000002 -0.56999999999999995 28 0.47360000000000002 -0.33079999999999998 29 0.20000000000000001 -0.79239999999999999 ]; % element number vertex 1 vertex 2 vertex 3 elements = [ 1 29 2 1 2 2 29 23 3 25 24 23 4 23 22 2 5 25 23 29 6 2 22 3 7 3 21 4 8 21 3 22 9 4 21 20 10 5 4 26 11 19 26 4 12 26 19 18 13 19 4 20 14 5 26 28 15 12 14 13 16 14 12 11 17 11 10 9 18 8 14 9 19 8 15 14 20 9 14 11 21 6 27 7 22 26 18 27 23 5 28 6 24 27 18 7 25 28 27 6 26 15 7 16 27 7 15 8 28 17 7 18 29 7 17 16 ]; %edge number, first node, second node, Dirichlet bc marker edges = [ 1 29 2 1 2 2 1 0 3 1 29 0 4 29 23 0 5 23 2 0 6 25 24 0 7 24 23 0 8 23 25 0 9 23 22 0 10 22 2 0 11 29 25 0 12 22 3 0 13 3 2 0 14 3 21 0 15 21 4 0 16 4 3 0 17 22 21 0 18 21 20 0 19 20 4 0 20 5 4 0 21 4 26 0 22 26 5 0 23 19 26 0 24 4 19 0 25 19 18 0 26 18 26 0 27 20 19 0 28 26 28 0 29 28 5 0 30 12 14 0 31 14 13 0 32 13 12 0 33 12 11 0 34 11 14 0 35 11 10 0 36 10 9 0 37 9 11 0 38 8 14 0 39 14 9 0 40 9 8 0 41 8 15 0 42 15 14 0 43 6 27 0 44 27 7 0 45 7 6 0 46 18 27 0 47 27 26 0 48 28 6 0 49 6 5 0 50 18 7 0 51 28 27 0 52 15 7 0 53 7 16 0 54 16 15 0 55 8 7 0 56 17 7 0 57 18 17 0 58 17 16 0 ]; % return end