diff --git a/main/elasticity.m b/main/elasticity.m index 8f3de38..f1afaf6 100644 --- a/main/elasticity.m +++ b/main/elasticity.m @@ -131,12 +131,4 @@ end axis([0 1 0 1]) -axis("square") - -[VI,eigen] = eigs(AII,MII,20,'smallestabs'); -%h+)!$Q99$mzXCdT - -V(internal_dofs) = VI(:,1); -V(boundary_dofs) = 0; - -eigen = diag(eigen) \ No newline at end of file +axis("square") \ No newline at end of file diff --git a/main/loop_vem_eig_elasticity.m b/main/loop_vem_eig_elasticity.m new file mode 100644 index 0000000..b6102bd --- /dev/null +++ b/main/loop_vem_eig_elasticity.m @@ -0,0 +1,186 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% FUNCTION: vem_eig +% +% Created by : U. Zerinati +% +%--------------------------------------------------------------------------------------------------- +% Purpose +% ======= +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +clear, close, clc + +fix_path(); + + +%% PARAMETERS OF THE PDE +matProps.E = 70; +matProps.nu = 0.2; +matProps.mu = matProps.E / (2 * (1 + matProps.nu)); +matProps.lambda = matProps.E * matProps.nu / ( (1 + matProps.nu) * (1 - 2*matProps.nu)); + +fprintf('\n\nParameters: ') %Print the parameters of the PDE +fprintf('\nYung modulus = %.2f\n',matProps.E) +fprintf('\nPoisson ratio = %.2f\n',matProps.nu) +fprintf('\nLame mu = %.2f\n',matProps.mu) +fprintf('\nLame lambda = %.2f\n',matProps.lambda) + +%% DEFINITION OF THE FUNCTIONS +[f1, f2, u1, u1_x, u1_y, u2, u2_x, u2_y] = problem_test_elasticity(4,matProps); %Obtain problem functions + +%% INFORMATION ON THE POLYNOMIALS +k = 1; %Degree of the polynomials used to solve the equation + +fprintf('\nPolynomials degree for solving the equation: %d',k) + +if (k ~= 1) + error("Actually, the method only works with k = 1") +end + +mesh = ["polygon_16.txt" "polygon_64.txt" "polygon_256.txt" "polygon_1024.txt" "polygon_4096.txt"]; +%mesh = ["star_2_0.1.txt" "star_4_0.1.txt" "star_8_0.1.txt" "star_16_0.1.txt"]; + +%% INITIALIZE MEMORY +eigLightning = zeros(numel(mesh),10); +errLightning = zeros(numel(mesh),10); +rateLightning = zeros(numel(mesh),10); +diamVec = zeros(numel(mesh),1); + +%% INFORMATION ON THE EIGENVALUE PROBLEM +eigExact = [1007.875293633368,... + 1007.875295841894,... + 1492.376174897711,... + 2165.996543952849,... + 2755.460012318369,... + 2882.723160229796,... + 2882.723160442838,... + 3529.696562815136,... + 4082.412001289926,... + 4082.412020260300]; + +numEigs = 10; meshIndex = 1; + +for mesh_filename = mesh + +tic + +%% PRINT INIT MESSAGE TO SCREEN +fprintf('\n\n[%.2f] Starting the method... \n',toc); + +%% READ THE MESH +fprintf('[%.2f] Reading mesh [%12s]...\n',toc, mesh_filename); +domainMesh = read_mesh(mesh_filename); %Read mesh + +%% OBTAIN INFORMATION ON THE MESH +fprintf('[%.2f] Obtaining information on the mesh... \n',toc); + +domainMesh = add_edges(domainMesh, k); %Add to the domainMesh structure the edges +boundary_vertex = domainMesh.boundary_nodes.all; %Extract boundary nodes + +domainMesh.boundary_edges = get_boundary_edges(domainMesh); %Get the boundary edges +domainMesh.internal_edges = setdiff(1:domainMesh.nedges, domainMesh.boundary_edges); %Get the internal edges + +[boundary_dofs, boundary_intdofs] = get_boundary_dofs(domainMesh, boundary_vertex, k); + +%% ASSEMBLYING LIGHTNING VEM MATRICES +fprintf('[%.2f] Assemblying element matrices...\n',toc); +[K_global, M_global, f_global, domainMesh] = elasticity_assembly(domainMesh, matProps, f1, f2, k); + +boundary_dofs = [boundary_vertex; boundary_vertex + domainMesh.nvertex]; +internal_dofs = setdiff(1:size(K_global,1), boundary_dofs); %Get the indexes of the internal edges + +diamVec(meshIndex) = domainMesh.diameter; + +AII = K_global(internal_dofs,internal_dofs); %Matrix internal - internal +MII = M_global(internal_dofs,internal_dofs); + +fprintf('[%.2f] Solving the eigenvalue problem...\n',toc); + +eigLightning(meshIndex,:) = eigs(AII,MII,numEigs,'smallestabs'); +errLightning(meshIndex,:) = abs(eigLightning(meshIndex,:) - eigExact); + +if meshIndex > 1 + for eigIndex = 1:numEigs + rate = polyfit(log(diamVec([meshIndex-1,meshIndex])), log(errLightning([meshIndex-1,meshIndex],eigIndex)),1); + rateLightning(meshIndex,eigIndex) = rate(1); + end +end + +meshIndex = meshIndex + 1; + + +end + +nus = [0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.4999]; +%% INITIALIZE MEMORY +eigPoisson = zeros(numel(mesh),numel(nus)); +errPoisson = zeros(numel(mesh),numel(nus)); +diamVec = zeros(numel(mesh),1); + +%% INFORMATION ON THE EIGENVALUE PROBLEM +eigExact = [1007.875293633368, 1043.4500230981541, 1109.5419789873215, 1230.7475057513004,... + 1295.3480604990973, 1256.4697090150764, 1221.442537732282] + +poissonIndex = 1; meshIndex = 1; + +for nu = nus + + matProps.E = 70; + matProps.nu = nu; + matProps.mu = matProps.E / (2 * (1 + matProps.nu)); + matProps.lambda = matProps.E * matProps.nu / ( (1 + matProps.nu) * (1 - 2*matProps.nu)); + + for mesh_filename = mesh + + [f1, f2, u1, u1_x, u1_y, u2, u2_x, u2_y] = problem_test_elasticity(4,matProps); + + %% INFORMATION ON THE POLYNOMIALS + k = 1; + %% PRINT INIT MESSAGE TO SCREEN + fprintf('\n\n[%.2f] Starting the method... \n',toc); + + %% READ THE MESH + fprintf('[%.2f] Reading mesh [%12s]...\n',toc, mesh_filename); + domainMesh = read_mesh(mesh_filename); %Read mesh + + %% OBTAIN INFORMATION ON THE MESH + fprintf('[%.2f] Obtaining information on the mesh... \n',toc); + + domainMesh = add_edges(domainMesh, k); %Add to the domainMesh structure the edges + boundary_vertex = domainMesh.boundary_nodes.all; %Extract boundary nodes + + domainMesh.boundary_edges = get_boundary_edges(domainMesh); %Get the boundary edges + domainMesh.internal_edges = setdiff(1:domainMesh.nedges, domainMesh.boundary_edges); %Get the internal edges + + [boundary_dofs, boundary_intdofs] = get_boundary_dofs(domainMesh, boundary_vertex, k); + + %% ASSEMBLYING LIGHTNING VEM MATRICES + fprintf('[%.2f] Assemblying element matrices...\n',toc); + [K_global, M_global, f_global, domainMesh] = elasticity_assembly(domainMesh, matProps, f1, f2, k); + + boundary_dofs = [boundary_vertex; boundary_vertex + domainMesh.nvertex]; + internal_dofs = setdiff(1:size(K_global,1), boundary_dofs); %Get the indexes of the internal edges + + diamVec(meshIndex) = domainMesh.diameter; + + AII = K_global(internal_dofs,internal_dofs); %Matrix internal - internal + MII = M_global(internal_dofs,internal_dofs); + + fprintf('[%.2f] Solving the eigenvalue problem...\n',toc); + + eigPoisson(meshIndex,poissonIndex) = eigs(AII,MII,1,'smallestabs'); + errPoisson(meshIndex,poissonIndex) = abs(eigPoisson(meshIndex,poissonIndex) - eigExact(poissonIndex)); + meshIndex = meshIndex + 1; + end + poissonIndex = poissonIndex + 1; +end + +loglog(diamVec,errPoisson(:,1), 'k*-','LineWidth',2, 'MarkerSize', 10) +hold on +loglog(diamVec,errPoisson(:,3), 'r*-','LineWidth',2, 'MarkerSize', 10) +loglog(diamVec,errPoisson(:,5),'b*-','LineWidth',2, 'MarkerSize', 10) +loglog(diamVec,errPoisson(:,7),'m*-','LineWidth',2, 'MarkerSize', 10) +grid on +axis square +legend('$\nu=0.2$', '$\nu=0.3$', '$\nu=0.4$', '$\nu=0.4999$', 'interpreter', 'latex') +xlabel('$h$','interpreter','latex') +ylabel('$|\lambda_1-\lambda_{1,h}|$','interpreter','latex') \ No newline at end of file diff --git a/main/loop_vem_eig_lighting.m b/main/loop_vem_eig_lighting.m new file mode 100644 index 0000000..c01707c --- /dev/null +++ b/main/loop_vem_eig_lighting.m @@ -0,0 +1,142 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% FUNCTION: vem_eig +% +% Created by : U. Zerinati +% +%--------------------------------------------------------------------------------------------------- +% Purpose +% ======= +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +clear, close, clc + +fix_path(); + + +%% PARAMETERS OF THE PDE +matProps.sigma = 0; %Reaction coefficient +matProps.epsilon = 1; %Diffusion coefficient +matProps.beta{1} = @(x,y) 0.*x + 0.*y; +matProps.beta{2} = @(x,y) 0.*x + 0.*y; + +beta1 = func2str(matProps.beta{1}); %Convert functions to string +beta2 = func2str(matProps.beta{2}); +beta1(1:6) = []; %Remove @(x,y) +beta2(1:6) = []; + +fprintf('\n\nParameters: ') %Print the parameters of the PDE +fprintf('\nEpsilon = %.2f\n',matProps.epsilon) +disp(['Beta = [ ' beta1 ' ; ' beta2 ' ]']) +fprintf('Sigma = %.2f\n',matProps.sigma) + +%% DEFINITION OF THE FUNCTIONS +[f, g, u, grad_u_x, grad_u_y] = problem_test_lighting(1,matProps); %Obtain problem functions + +%% INFORMATION ON THE POLYNOMIALS +k = 1; %Degree of the polynomials used to solve the equation + +%polynomial = get_polynomial_info(k); +fprintf('\nPolynomials degree for solving the equation: %d',k) + +if (k ~= 1) + error("Actually, the method only works with k = 1") +end + +%mesh = ["polygon_16.txt" "polygon_64.txt" "polygon_256.txt" "polygon_1024.txt"]; +%mesh = ["star_2_0.1.txt" "star_4_0.1.txt" "star_8_0.1.txt" "star_16_0.1.txt"]; +mesh = ["polygon_1024.txt"]; + +%% INFORMATION ON THE EIGENVALUE PROBLEM +eigExact = [2 5 5 8 10 10 13 13 17 17 18 20 20 25 25 26 26 29 29 32]; +numEigs = 20; meshIndex = 1; + +%% INITIALIZE MEMORY +eigLightning = zeros(numel(mesh),numEigs); +errLightning = zeros(numel(mesh),numEigs); +rateLightning = zeros(numel(mesh),numEigs); +diamVec = zeros(numel(mesh),1); + + +for mesh_filename = mesh + +tic + +%% PRINT INIT MESSAGE TO SCREEN +fprintf('\n\n[%.2f] Starting the method... \n',toc); + +%% READ THE MESH +fprintf('[%.2f] Reading mesh [%12s]...\n',toc, mesh_filename); +domainMesh = read_mesh(mesh_filename); %Read mesh + +%% OBTAIN INFORMATION ON THE MESH +fprintf('[%.2f] Obtaining information on the mesh... \n',toc); + +domainMesh = add_edges(domainMesh, k); %Add to the domainMesh structure the edges +boundary_vertex = domainMesh.boundary_nodes.all; %Extract boundary nodes + +domainMesh.boundary_edges = get_boundary_edges(domainMesh); %Get the boundary edges +domainMesh.internal_edges = setdiff(1:domainMesh.nedges, domainMesh.boundary_edges); %Get the internal edges + +[boundary_dofs, boundary_intdofs] = get_boundary_dofs(domainMesh, boundary_vertex, k); + +%% ASSEMBLYING LIGHTNING VEM MATRICES +fprintf('[%.2f] Assemblying element matrices...\n',toc); +[K_global, M_global, f_global, domainMesh] = vem_lighting_assembly(domainMesh, matProps, f, k); + +internal_dofs = setdiff(1:size(K_global,1), boundary_dofs); %Get the indexes of the internal edges + +diamVec(meshIndex) = domainMesh.diameter; + +AII = K_global(internal_dofs,internal_dofs); %Matrix internal - internal +MII = M_global(internal_dofs,internal_dofs); + +fprintf('[%.2f] Solving the eigenvalue problem...\n',toc); + +eigLightning(meshIndex,:) = eigs(AII,MII,numEigs,'smallestabs') ./ pi ./ pi; +errLightning(meshIndex,:) = abs(eigLightning(meshIndex,:) - eigExact); + +if meshIndex > 1 + for eigIndex = 1:numEigs + rate = polyfit(log(diamVec([meshIndex-1,meshIndex])), log(errLightning([meshIndex-1,meshIndex],eigIndex)),1); + rateLightning(meshIndex,eigIndex) = rate(1); + end +end + +meshIndex = meshIndex + 1; + + +end + +polynomial = get_polynomial_info(k); +%% ASSEMBLYING VEM MATRICES +fprintf('[%.2f] Assemblying element matrices...\n',toc); +[A1, A2, B1, B2, f_global, domainMesh] = vem_adr_assembly(domainMesh, matProps, polynomial, f); + +A1 = A1(internal_dofs,internal_dofs); +A2 = A2(internal_dofs,internal_dofs); +B1 = B1(internal_dofs,internal_dofs); +B2 = B2(internal_dofs,internal_dofs); + +numTest = 51; +a = linspace(0.00,10,numTest); +beta = 5; + +eigVEM = zeros(numEigs,numTest); +for i = 1:numTest + + + A = A1 + a(i)*A2; + B = B1 + beta*B2; + + eigVEM(:,i) = eigs(A,B,numEigs,'smallestabs') ./ pi ./ pi; +end + +for i = 1:numEigs + grid on + plot(a,eigVEM(i,:),'o-', 'Linewidth', 1.5, 'MarkerSize', 5); + hold on + plot(a(end)+10/numTest,eigLightning(1,i),'k*') + plot(a(end)+10/numTest,eigExact(i),'ro') + ylim([0,40]) + xlabel('$\alpha$','interpreter','latex') + ylabel('$\lambda_i$','interpreter','latex') +end \ No newline at end of file diff --git a/main/vem_adr.m b/main/vem_adr.m index b4cf067..cf809db 100644 --- a/main/vem_adr.m +++ b/main/vem_adr.m @@ -102,15 +102,10 @@ A = A1 + A2; AII = A(internal_dofs,internal_dofs); %Matrix internal - internal - A1 = A1(internal_dofs,internal_dofs); A2 = A2(internal_dofs,internal_dofs); B1 = B1(internal_dofs,internal_dofs); B2 = B2(internal_dofs,internal_dofs); -save("A1.mat",'A1') -save("A2.mat",'A2') -save("B1.mat",'B1') -save("B2.mat",'B2') AIB = A(internal_dofs, boundary_dofs); %Matrix internal - boundary fI = f_global(internal_dofs); %Interior load term diff --git a/main/vem_lighting.m b/main/vem_lighting.m index 0014676..72fd39a 100644 --- a/main/vem_lighting.m +++ b/main/vem_lighting.m @@ -1,17 +1,17 @@ -% vem_lighting: This function computes the numerical solution of the PDE -% -eps * div(grad(u)) + beta * grad(u) + sigma * u = f -% using virtual element method (VEM). +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% FUNCTION: vem_lighting % -% The basis functions of the virtual element space are computed using the lighting technique. +% Created by : M. Trezzi % -% The user can set the parameters of the PDE in lines 29-32. +%--------------------------------------------------------------------------------------------------- +% Purpose +% ======= +% vem_lighting: this function computes the numerical solution of the pde +% -eps * div(grad(u)) + beta * grad(u) + sigma * u = f +% using virtual element method (vem). % -% The test problem is set in line 43. -% -% The mesh is selected in line 61. The meshes are constructed using VEMLAB 2.4 -% (https://camlab.cl/software/vemlab/) and the functions plot_mesh2d and -% read_mesh are from that software. - +% the basis functions of the virtual element space are computed using the lighting technique. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% INITIALIZATION clear; close; clc; @@ -111,7 +111,4 @@ fprintf("\n[%.2f] Errors computed: ",toc) fprintf("\nL2 norm : %f", errL2) -fprintf("\nH1 seminorm: %f\n", errH1) - -save("A.mat",'AII') -save("B.mat",'MII') \ No newline at end of file +fprintf("\nH1 seminorm: %f\n", errH1) \ No newline at end of file