From 92574b5d5dd5cfb5a5d608b61c3074922e98399f Mon Sep 17 00:00:00 2001 From: Nicole Stock Date: Tue, 8 Sep 2020 17:46:47 -0400 Subject: [PATCH] First commit --- basis_functions_weighted_p2.m | 63 +++++++++++++ create_b_p2.m | 67 +++++++++++++ errors_exact_weighted_p2.m | 145 ++++++++++++++++++++++++++++ errors_no_exact_weighted_p2.m | 166 +++++++++++++++++++++++++++++++++ find_midpoints.m | 66 +++++++++++++ main.m | 64 +++++++++++++ main.m~ | 64 +++++++++++++ main_no_exact.m | 73 +++++++++++++++ main_no_exact.m~ | 78 ++++++++++++++++ prolongation_matrix.m | 70 ++++++++++++++ stiffness_matrix_weighted_p2.m | 131 ++++++++++++++++++++++++++ triquad.m | 49 ++++++++++ 12 files changed, 1036 insertions(+) create mode 100644 basis_functions_weighted_p2.m create mode 100644 create_b_p2.m create mode 100644 errors_exact_weighted_p2.m create mode 100644 errors_no_exact_weighted_p2.m create mode 100644 find_midpoints.m create mode 100644 main.m create mode 100644 main.m~ create mode 100644 main_no_exact.m create mode 100644 main_no_exact.m~ create mode 100644 prolongation_matrix.m create mode 100644 stiffness_matrix_weighted_p2.m create mode 100644 triquad.m diff --git a/basis_functions_weighted_p2.m b/basis_functions_weighted_p2.m new file mode 100644 index 0000000..2135309 --- /dev/null +++ b/basis_functions_weighted_p2.m @@ -0,0 +1,63 @@ +function basis = basis_functions_weighted_p2(p,t,p2,t2) +% BASIS_FUNCTIONS_WEIGHTED_P2 - Create a piecewise basis function for each +% node of a triangulation with weight k +% +% Syntax: +% basis = basis_functions_weighted_p2(p,t,p2,t2) +% +% Inputs: +% p - a 2-by-NumNodes matrix representing nodal coordinates. +% t - a matrix representing the element connectivity in terms of node +% IDs. The end row of T represents the geometry face ID to which the +% element belongs +% k - given weight +% +% Outputs: +% basis - a matrix representing piece-wise basis functions for each node +% in each triangle. basis(i,:,T) represents the pieceiwise basis +% function for the ith node in triangle T. +% +% Author: Nicole Stock +% Date: Spring 2020 + +[~,triangles] = size(t); +basis = zeros(6,6,triangles); + +for T = 1:triangles + poly = zeros(6,6); + + % for each row (node) + for i = 1:3 + % 'regular' node + node = t(i, T); + % get r,z coordinates + r = p(1,node); + z = p(2,node); + + % fill in polynomial + poly(i,1) = r^2; + poly(i,2) = r*z; + poly(i,3) = z^2; + poly(i,4) = r; + poly(i,5) = z; + poly(i,6) = 1; + + % midpoint node + node = t2(i,T); + % get r,z coordinates + r = p2(1,node); + z = p2(2,node); + % fill in polynomial + poly(i+3,1) = r^2; + poly(i+3,2) = r*z; + poly(i+3,3) = z^2; + poly(i+3,4) = r; + poly(i+3,5) = z; + poly(i+3,6) = 1; + end + + I = eye(6); + basis(:,:,T) = poly\I; +end + +% end diff --git a/create_b_p2.m b/create_b_p2.m new file mode 100644 index 0000000..283c40b --- /dev/null +++ b/create_b_p2.m @@ -0,0 +1,67 @@ +function b = create_b_p2(p,t,p2,t2,basis,f) +% CREATE_B_P2 - Create vector b +% +% Syntax: +% b = create_b_p2(p,e,t,basis,f) +% +% Inputs: +% p - a 2xNumNodes matrix representing nodal coordinates. +% t - a 4xNumTriangles matrix representing the element connectivity in +% terms of node IDs. The end row of T represents the geometry face ID +% to which the element belongs. +% p2 - a 2xNumNodes matrix representing midpoint nodal coordinates. +% t2 - a 4xNumTriangles matrix representing the element connectivity in +% terms of node IDs. The three node IDs in a column are the three +% midpoints of the node IDS in corresponding column in t. +% basis - a 3x3xNumTriangles matrix representing piece-wise basis +% functions for each node in each triangle. basis(i,:,k) represents +% the pieceiwise basis function for the ith node in triangle k. +% +% Outputs: +% b - vector such that stiffness_matrix * b = solution. +% +% Author: Nicole Stock +% Date: Spring 2020 + +[~,triangles] = size(t); +[~,nodes] = size(p); +i_vec = zeros(1,triangles*6); +j_vec = ones(1,triangles*6); +s_vec = zeros(1,triangles*6); +index = 1; + +for T = 1:triangles + + % get coordinates of triangle T (Tth col of t) + coordinates = zeros(3,2); + for n = 1:3 + node = t(n,T); + % get x,y coordinates + coordinates(n,:) = p(:,node); + end + + [X,Y,Wx,Wy] = triquad(7, coordinates); + + for i = 1:6 + + I = basis(:,i,T); + + integrand =@(x,y) feval(f,x,y).*(I(1).*x.^2 + I(2).*x.*y ... + + I(3).*y.^2 + I(4).*x + I(5).*y + I(6)); + b_i = Wx'*feval(integrand,X,Y)*Wy; + + if i >= 4 + global_i = t2(i-3,T) + nodes; + else + global_i = t(i,T); + end + + i_vec(index) = global_i; + s_vec(index) = b_i; + index = index + 1; + end +end + +[~,n] = size(p); +[~,n2] = size(p2); +b = sparse(i_vec,j_vec,s_vec,n+n2,1); \ No newline at end of file diff --git a/errors_exact_weighted_p2.m b/errors_exact_weighted_p2.m new file mode 100644 index 0000000..1f80004 --- /dev/null +++ b/errors_exact_weighted_p2.m @@ -0,0 +1,145 @@ +function [err,grad_err,max_err] = errors_exact_weighted_p2(p,t,p2,t2,basis,u_h,k) + +% ERRORS_EXACT_WEIGHTED_p2 - Calculate the errors of our solution u_h +% compared to the exact solution u. +% +% Syntax: +% [err,grad_err,max_err] = errors_exact_weighted_p2(p,e,t,u_h_km1,u_h_k) +% +% Inputs: +% p - a 2xNumNodes matrix representing nodal coordinates. +% t - a 4xNumTriangles matrix representing the element connectivity in +% terms of node IDs. The end row of T represents the geometry face ID +% to which the element belongs +% basis - a 3x3xNumTriangles matrix representing piece-wise basis +% functions for each node in each triangle. basis(i,:,k) represents +% the pieceiwise basis function for the ith node in triangle k. +% u_h_km1 - approximate solution for mesh level k-1 +% u_h_k - approximate solution for mesh level k +% k - given weight +% +% Outputs: +% err - L2 error +% grad_err - L2 gradient error +% max_err - max error +% +% Author: Nicole Stock +% Date: Spring 2020 + +[~,triangles] = size(t); +[~,nodes] = size(p); +[~,mid_nodes] = size(p2); + +integral = 0; +grad_integral = 0; + +for T = 1:triangles + + % get coordinates of triangle T + coordinates = zeros(3,2); + for n = 1:3 + node = t(n,T); + % get x,y coordinates + coordinates(n,:) = p(:,node); + end + + [u,grad_u_r,grad_u_z] = get_exact_y_solutions(coordinates(1,1),coordinates(1,2),k); + + [R,Z,Wr,Wz] = triquad(7, coordinates); + + b1 = basis(:,1,T); + b2 = basis(:,2,T); + b3 = basis(:,3,T); + b4 = basis(:,4,T); + b5 = basis(:,5,T); + b6 = basis(:,6,T); + + n1 = t(1,T); + n2 = t(2,T); + n3 = t(3,T); + n4 = t2(1,T); + n5 = t2(2,T); + n6 = t2(3,T); + + if k == 0 + % find L2 Error for u_h_km1 & u_h_k + approx =@(r,z) u_h(n1).*(b1(1).*r.^2 + b1(2).*r.*z ... + + b1(3).*z.^2 + b1(4).*r + b1(5).*z + b1(6)) ... + + u_h(n2).*(b2(1).*r.^2 + b2(2).*r.*z ... + + b2(3).*z.^2 + b2(4).*r + b2(5).*z + b2(6)) ... + + u_h(n3).*(b3(1).*r.^2 + b3(2).*r.*z ... + + b3(3).*z.^2 + b3(4).*r + b3(5).*z + b3(6)) ... + + u_h(n4 + nodes).*(b4(1).*r.^2 + b4(2).*r.*z ... + + b4(3).*z.^2 + b4(4).*r + b4(5).*z + b4(6)) ... + + u_h(n5 + nodes).*(b5(1).*r.^2 + b5(2).*r.*z ... + + b5(3).*z.^2 + b5(4).*r + b5(5).*z + b5(6)) ... + + u_h(n6 + nodes).*(b6(1).*r.^2 + b6(2).*r.*z ... + + b6(3).*z.^2 + b6(4).*r + b6(5).*z + b6(6)); + else + % find L2 Error for u_h_km1 & u_h_k + approx =@(r,z) u_h(n1).*(r./k).*(b1(1).*r.^2 + b1(2).*r.*z ... + + b1(3).*z.^2 + b1(4).*r + b1(5).*z + b1(6)) ... + + u_h(n2).*(r./k).*(b2(1).*r.^2 + b2(2).*r.*z ... + + b2(3).*z.^2 + b2(4).*r + b2(5).*z + b2(6)) ... + + u_h(n3).*(r./k).*(b3(1).*r.^2 + b3(2).*r.*z ... + + b3(3).*z.^2 + b3(4).*r + b3(5).*z + b3(6)) ... + + u_h(n4 + nodes).*(r./k).*(b4(1).*r.^2 + b4(2).*r.*z ... + + b4(3).*z.^2 + b4(4).*r + b4(5).*z + b4(6)) ... + + u_h(n5 + nodes).*(r./k).*(b5(1).*r.^2 + b5(2).*r.*z ... + + b5(3).*z.^2 + b5(4).*r + b5(5).*z + b5(6)) ... + + u_h(n6 + nodes).*(r./k).*(b6(1).*r.^2 + b6(2).*r.*z ... + + b6(3).*z.^2 + b6(4).*r + b6(5).*z + b6(6)); + end + + integrand =@(r,z) ((u(r,z) - approx(r,z)).^2).*r; + + integral = integral + Wr'*feval(integrand,R,Z)*Wz; + + % find Grad L2 Error for u_h_km1 & u_h_k + grad_approx_r =@(r,z) u_h(n1).*(2.*b1(1).*r + b1(2).*z ... + + b1(4)) + u_h(n2).*(2.*b2(1).*r + b2(2).*z + b2(4)) ... + + u_h(n3).*(2.*b3(1).*r + b3(2).*z + b3(4)) ... + + u_h(n4 + nodes).*(2.*b4(1).*r + b4(2).*z + b4(4)) ... + + u_h(n5 + nodes).*(2.*b5(1).*r + b5(2).*z + b5(4)) ... + + u_h(n6 + nodes).*(2.*b6(1).*r + b6(2).*z + b6(4)); + grad_approx_z =@(r,z) u_h(n1).*(b1(2).*r + 2.*b1(3).*z ... + + b1(5)) + u_h(n2).*(b2(2).*r + 2.*b2(3).*z + b2(5)) ... + + u_h(n3).*(b3(2).*r + 2.*b3(3).*z + b3(5)) ... + + u_h(n4 + nodes).*(b4(2).*r + 2.*b4(3).*z + b4(5)) ... + + u_h(n5 + nodes).*(b5(2).*r + 2.*b5(3).*z + b5(5)) ... + + u_h(n6 + nodes).*(b6(2).*r + 2.*b6(3).*z + b6(5)); + + grad_integrand =@(r,z) ((grad_u_r(r,z) - grad_approx_r(r,z)).^2 ... + + (grad_u_z(r,z) - grad_approx_z(r,z)).^2).*r; + + grad_integral = grad_integral + Wr'*feval(grad_integrand,R,Z)*Wz; +end + +% find Max Error +max_err = -Inf; +for n = 1:nodes + r = p(1,n); + z = p(2,n); + + diff = abs(u(r,z) - u_h(n)); + + if diff > max_err + max_err = diff; + end +end + +for n = 1:mid_nodes + r = p2(1,n); + z = p2(2,n); + + diff = abs(u(r,z) - u_h(n + nodes)); + + if diff > max_err + max_err = diff; + end +end + +err = sqrt(integral); +grad_err = sqrt(grad_integral); + +% end \ No newline at end of file diff --git a/errors_no_exact_weighted_p2.m b/errors_no_exact_weighted_p2.m new file mode 100644 index 0000000..3129ddc --- /dev/null +++ b/errors_no_exact_weighted_p2.m @@ -0,0 +1,166 @@ +function [err,grad_err,max_err] = errors_no_exact_weighted_p2(p,t,p2,t2,basis,u_h_km1,u_h_k,k) +% ERRORS_NO_EXACT_WEIGHTED_p2 - Calculate the errors of our solution u_h_km1 +% compared to the approximate solution for the next mesh level (u_h_k). +% +% Syntax: +% [err,grad_err,max_err] = errors_no_exact_weighted_p2(p,e,t,u_h_km1,u_h_k) +% +% Inputs: +% p - a 2xNumNodes matrix representing nodal coordinates. +% t - a 4xNumTriangles matrix representing the element connectivity in +% terms of node IDs. The end row of T represents the geometry face ID +% to which the element belongs +% basis - a 3x3xNumTriangles matrix representing piece-wise basis +% functions for each node in each triangle. basis(i,:,k) represents +% the pieceiwise basis function for the ith node in triangle k. +% u_h_km1 - approximate solution for mesh level k-1 +% u_h_k - approximate solution for mesh level k +% k - given weight +% +% Outputs: +% err - L2 error +% grad_err - L2 gradient error +% max_err - max error +% +% Author: Nicole Stock +% Date: Spring 2020 + +[~,triangles] = size(t); +[~,nodes] = size(p); +[~,mid_nodes] = size(p2); + +integral = 0; +grad_integral = 0; + +for T = 1:triangles + + % get coordinates of triangle T + coordinates = zeros(3,2); + for n = 1:3 + node = t(n,T); + % get x,y coordinates + coordinates(n,:) = p(:,node); + end + + [R,Z,Wr,Wz] = triquad(7, coordinates); + + b1 = basis(:,1,T); + b2 = basis(:,2,T); + b3 = basis(:,3,T); + b4 = basis(:,4,T); + b5 = basis(:,5,T); + b6 = basis(:,6,T); + + n1 = t(1,T); + n2 = t(2,T); + n3 = t(3,T); + n4 = t2(1,T); + n5 = t2(2,T); + n6 = t2(3,T); + + if k == 0 + % find L2 Error for u_h_km1 & u_h_k + approx_km1 =@(r,z) u_h_km1(n1).*(b1(1).*r.^2 + b1(2).*r.*z ... + + b1(3).*z.^2 + b1(4).*r + b1(5).*z + b1(6)) ... + + u_h_km1(n2).*(b2(1).*r.^2 + b2(2).*r.*z ... + + b2(3).*z.^2 + b2(4).*r + b2(5).*z + b2(6)) ... + + u_h_km1(n3).*(b3(1).*r.^2 + b3(2).*r.*z ... + + b3(3).*z.^2 + b3(4).*r + b3(5).*z + b3(6)) ... + + u_h_km1(n4 + nodes).*(b4(1).*r.^2 + b4(2).*r.*z ... + + b4(3).*z.^2 + b4(4).*r + b4(5).*z + b4(6)) ... + + u_h_km1(n5 + nodes).*(b5(1).*r.^2 + b5(2).*r.*z ... + + b5(3).*z.^2 + b5(4).*r + b5(5).*z + b5(6)) ... + + u_h_km1(n6 + nodes).*(b6(1).*r.^2 + b6(2).*r.*z ... + + b6(3).*z.^2 + b6(4).*r + b6(5).*z + b6(6)); + approx_k =@(r,z) u_h_k(n1).*(b1(1).*r.^2 + b1(2).*r.*z ... + + b1(3).*z.^2 + b1(4).*r + b1(5).*z + b1(6)) ... + + u_h_k(n2).*(b2(1).*r.^2 + b2(2).*r.*z ... + + b2(3).*z.^2 + b2(4).*r + b2(5).*z + b2(6)) ... + + u_h_k(n3).*(b3(1).*r.^2 + b3(2).*r.*z ... + + b3(3).*z.^2 + b3(4).*r + b3(5).*z + b3(6)) ... + + u_h_k(n4 + nodes).*(b4(1).*r.^2 + b4(2).*r.*z ... + + b4(3).*z.^2 + b4(4).*r + b4(5).*z + b4(6)) ... + + u_h_k(n5 + nodes).*(b5(1).*r.^2 + b5(2).*r.*z ... + + b5(3).*z.^2 + b5(4).*r + b5(5).*z + b5(6)) ... + + u_h_k(n6 + nodes).*(b6(1).*r.^2 + b6(2).*r.*z ... + + b6(3).*z.^2 + b6(4).*r + b6(5).*z + b6(6)); + else + % find L2 Error for u_h_km1 & u_h_k + approx_km1 =@(r,z) u_h_km1(n1).*(r./k).*(b1(1).*r.^2 + b1(2).*r.*z ... + + b1(3).*z.^2 + b1(4).*r + b1(5).*z + b1(6)) ... + + u_h_km1(n2).*(r./k).*(b2(1).*r.^2 + b2(2).*r.*z ... + + b2(3).*z.^2 + b2(4).*r + b2(5).*z + b2(6)) ... + + u_h_km1(n3).*(r./k).*(b3(1).*r.^2 + b3(2).*r.*z ... + + b3(3).*z.^2 + b3(4).*r + b3(5).*z + b3(6)) ... + + u_h_km1(n4 + nodes).*(r./k).*(b4(1).*r.^2 + b4(2).*r.*z ... + + b4(3).*z.^2 + b4(4).*r + b4(5).*z + b4(6)) ... + + u_h_km1(n5 + nodes).*(r./k).*(b5(1).*r.^2 + b5(2).*r.*z ... + + b5(3).*z.^2 + b5(4).*r + b5(5).*z + b5(6)) ... + + u_h_km1(n6 + nodes).*(r./k).*(b6(1).*r.^2 + b6(2).*r.*z ... + + b6(3).*z.^2 + b6(4).*r + b6(5).*z + b6(6)); + approx_k =@(r,z) u_h_k(n1).*(r./k).*(b1(1).*r.^2 + b1(2).*r.*z ... + + b1(3).*z.^2 + b1(4).*r + b1(5).*z + b1(6)) ... + + u_h_k(n2).*(r./k).*(b2(1).*r.^2 + b2(2).*r.*z ... + + b2(3).*z.^2 + b2(4).*r + b2(5).*z + b2(6)) ... + + u_h_k(n3).*(r./k).*(b3(1).*r.^2 + b3(2).*r.*z ... + + b3(3).*z.^2 + b3(4).*r + b3(5).*z + b3(6)) ... + + u_h_k(n4 + nodes).*(r./k).*(b4(1).*r.^2 + b4(2).*r.*z ... + + b4(3).*z.^2 + b4(4).*r + b4(5).*z + b4(6)) ... + + u_h_k(n5 + nodes).*(r./k).*(b5(1).*r.^2 + b5(2).*r.*z ... + + b5(3).*z.^2 + b5(4).*r + b5(5).*z + b5(6)) ... + + u_h_k(n6 + nodes).*(r./k).*(b6(1).*r.^2 + b6(2).*r.*z ... + + b6(3).*z.^2 + b6(4).*r + b6(5).*z + b6(6)); + end + + integrand =@(r,z) ((approx_km1(r,z) - approx_k(r,z)).^2).*r; + + integral = integral + Wr'*feval(integrand,R,Z)*Wz; + + % find Grad L2 Error for u_h_km1 & u_h_k + grad_approx_r_km1 =@(r,z) u_h_km1(n1).*(2.*b1(1).*r + b1(2).*z ... + + b1(4)) + u_h_km1(n2).*(2.*b2(1).*r + b2(2).*z + b2(4)) ... + + u_h_km1(n3).*(2.*b3(1).*r + b3(2).*z + b3(4)) ... + + u_h_km1(n4 + nodes).*(2.*b4(1).*r + b4(2).*z + b4(4)) ... + + u_h_km1(n5 + nodes).*(2.*b5(1).*r + b5(2).*z + b5(4)) ... + + u_h_km1(n6 + nodes).*(2.*b6(1).*r + b6(2).*z + b6(4)); + grad_approx_z_km1 =@(r,z) u_h_km1(n1).*(b1(2).*r + 2.*b1(3).*z ... + + b1(5)) + u_h_km1(n2).*(b2(2).*r + 2.*b2(3).*z + b2(5)) ... + + u_h_km1(n3).*(b3(2).*r + 2.*b3(3).*z + b3(5)) ... + + u_h_km1(n4 + nodes).*(b4(2).*r + 2.*b4(3).*z + b4(5)) ... + + u_h_km1(n5 + nodes).*(b5(2).*r + 2.*b5(3).*z + b5(5)) ... + + u_h_km1(n6 + nodes).*(b6(2).*r + 2.*b6(3).*z + b6(5)); + + grad_approx_r_k =@(r,z) u_h_k(n1).*(2.*b1(1).*r + b1(2).*z ... + + b1(4)) + u_h_k(n2).*(2.*b2(1).*r + b2(2).*z + b2(4)) ... + + u_h_k(n3).*(2.*b3(1).*r + b3(2).*z + b3(4)) ... + + u_h_k(n4 + nodes).*(2.*b4(1).*r + b4(2).*z + b4(4)) ... + + u_h_k(n5 + nodes).*(2.*b5(1).*r + b5(2).*z + b5(4)) ... + + u_h_k(n6 + nodes).*(2.*b6(1).*r + b6(2).*z + b6(4)); + grad_approx_z_k =@(r,z) u_h_k(n1).*(b1(2).*r + 2.*b1(3).*z ... + + b1(5)) + u_h_k(n2).*(b2(2).*r + 2.*b2(3).*z + b2(5)) ... + + u_h_k(n3).*(b3(2).*r + 2.*b3(3).*z + b3(5)) ... + + u_h_k(n4 + nodes).*(b4(2).*r + 2.*b4(3).*z + b4(5)) ... + + u_h_k(n5 + nodes).*(b5(2).*r + 2.*b5(3).*z + b5(5)) ... + + u_h_k(n6 + nodes).*(b6(2).*r + 2.*b6(3).*z + b6(5)); + + grad_integrand =@(r,z) ((grad_approx_r_km1(r,z) ... + - grad_approx_r_k(r,z)).^2 + (grad_approx_z_km1(r,z) ... + - grad_approx_z_k(r,z)).^2).*r; + + grad_integral = grad_integral + Wr'*feval(grad_integrand,R,Z)*Wz; +end + +% find Max Error +max_err = -Inf; +for n = 1:nodes+mid_nodes + diff = abs(u_h_km1(n) - u_h_k(n)); + + if diff > max_err + max_err = diff; + end +end + +err = sqrt(integral); +grad_err = sqrt(grad_integral); + +% end \ No newline at end of file diff --git a/find_midpoints.m b/find_midpoints.m new file mode 100644 index 0000000..c22e8ff --- /dev/null +++ b/find_midpoints.m @@ -0,0 +1,66 @@ +function [p2,t2] = find_midpoints(p,t) +% FIND_MIDPOINTS - Find the midpoints between adjacent nodes in the mesh +% +% Syntax: +% [p2,t2] = find_midpoints(p,t) +% Inputs: +% p - a 2xNumNodes matrix representing nodal coordinates. +% t - a 4xNumTriangles matrix representing the element connectivity in +% terms of node IDs. The end row of T represents the geometry face +% ID to which the element belongs. +% +% Outputs: +% p2 - a 2xNumNodes matrix representing midpoint nodal coordinates. +% t2 - a 4xNumTriangles matrix representing the element connectivity in +% terms of node IDs. The three node IDs in a column are the three +% midpoints of the node IDS in corresponding column in t. +% +% Author: Nicole Stock +% Date: Spring 2020 + +[~,triangles] = size(t); +[~,nodes] = size(p); +p2 = NaN(2,nodes); +t2 = NaN(3,triangles); +counter_p2 = 1; + +% for each triangle +for T = 1:triangles + % get coordinate points of the 3 nodes in triangle T + xs = NaN(3,1); + ys = NaN(3,1); + for i = 1:3 + node = t(i, T); + % get x,y coordinates + xs(i) = p(1,node); + ys(i) = p(2,node); + end + + % find midpoints for each side + mid_xs = NaN(3,1); + mid_ys = NaN(3,1); + % x values + mid_xs(1) = (xs(1) + xs(2))/2; + mid_xs(2) = (xs(2) + xs(3))/2; + mid_xs(3) = (xs(3) + xs(1))/2; + % y values + mid_ys(1) = (ys(1) + ys(2))/2; + mid_ys(2) = (ys(2) + ys(3))/2; + mid_ys(3) = (ys(3) + ys(1))/2; + + % check for duplicates + [isfound,index] = ismember([mid_xs,mid_ys], p2.', 'rows'); + % add midpoints to p2 and t2 + for i = 1:3 + if ~isfound(i) + p2(1,counter_p2) = mid_xs(i); + p2(2,counter_p2) = mid_ys(i); + t2(i,T) = counter_p2; + counter_p2 = counter_p2 + 1; + else + t2(i,T) = index(i); + end + end +end + +end \ No newline at end of file diff --git a/main.m b/main.m new file mode 100644 index 0000000..e1da870 --- /dev/null +++ b/main.m @@ -0,0 +1,64 @@ +function [err,grad_err,max_err] = main(u,gd,sf,ns,mesh_level,opt) +%MAIN Summary of this function goes here +% Detailed explanation goes here +% +% Author: Nicole Stock +% Date: Fall 2020 + +model=createpde(1); +g=decsg(gd,sf,ns); +geometryFromEdges(model,g); +[p,e,t]=initmesh(g,'hmax',inf); +%pdemesh(p,e,t, 'NodeLabels','on', 'ElementLabels','on'); + +% Init error vectors +err = zeros(1,mesh_level); +grad_err = zeros(1,mesh_level); +max_err = zeros(1,mesh_level); + +if mesh_level > 1 + % To ensure we refine every triangle the same + [~,num_node]=size(p); + it=zeros(1,num_node); + for i=1:num_node + it(i)=i; + end + + [p2,t2] = find_midpoints(p,t); + + [basis,Qh] = solve(p,p2,e,t,t2,f); + [err(1),grad_err(1),max_err(1)] = errors_exact_p2(p,t,p2,t2,basis,Qh,u,grad_u_x,grad_u_y); + % errors_exact_weighted_p2(...) + + for i = 2:mesh_level + [p,e,t]=refinemesh(g,p,e,t,it,'regular'); + %pdemesh(p,e,t, 'NodeLabels','on', 'ElementLabels','on'); + + [p2,t2] = find_midpoints(p,t); + + [basis,Qh] = solve(p,p2,e,t,t2,f); + [err(i),grad_err(i),max_err(i)] = errors_exact_p2(p,t,p2,t2,basis,Qh,u,grad_u_x,grad_u_y); + + end + display_errors(err,grad_err,max_err) + +else + [p2,t2] = find_midpoints(p,t); + + [basis,Qh] = solve(p,p2,e,t,t2,f); + [err(1),grad_err(1),max_err(1)] = errors_exact_p2(p,t,p2,t2,basis,Qh,u,grad_u_x,grad_u_y); +end + +% end main +end + +% subfunction +function [basis,Qh] = solve(p,p2,e,t,t2,u) + basis = basis_functions_weighted_p2(p,t,p2,t2); + S = stiffness_matrix_weighted_p2(p,t,p2,t2,basis); + b = create_b_p2(p,t,p2,t2,basis,u); + Qh = S\b; + + figure(); + pdeplot([p,p2],e,t, 'XYData',Qh, 'ZData', Qh, 'Mesh', 'on'); +end diff --git a/main.m~ b/main.m~ new file mode 100644 index 0000000..13a1696 --- /dev/null +++ b/main.m~ @@ -0,0 +1,64 @@ +function [err,grad_err,max_err] = main(u,gd,sf,ns,mesh_level,opt) +%MAIN Summary of this function goes here +% Detailed explanation goes here +% +% Author: Nicole Stock +% Date: Fall 2020 + +model=createpde(1); +g=decsg(gd,sf,ns); +geometryFromEdges(model,g); +[p,e,t]=initmesh(g,'hmax',inf); +%pdemesh(p,e,t, 'NodeLabels','on', 'ElementLabels','on'); + +% Init error vectors +err = zeros(1,mesh_level); +grad_err = zeros(1,mesh_level); +max_err = zeros(1,mesh_level); + +if mesh_level > 1 + % To ensure we refine every triangle the same + [~,num_node]=size(p); + it=zeros(1,num_node); + for i=1:num_node + it(i)=i; + end + + [p2,t2] = find_midpoints(p,t); + + [basis,Qh] = solve(p,p2,e,t,t2,f); + [err(1),grad_err(1),max_err(1)] = errors_exact_p2(p,t,p2,t2,basis,Qh,u,grad_u_x,grad_u_y); + % errors_no_exact_weighted_p2(p,t,p2,t2,basis,u_h_km1,u_h_k,k) + + for i = 2:mesh_level + [p,e,t]=refinemesh(g,p,e,t,it,'regular'); + %pdemesh(p,e,t, 'NodeLabels','on', 'ElementLabels','on'); + + [p2,t2] = find_midpoints(p,t); + + [basis,Qh] = solve(p,p2,e,t,t2,f); + [err(i),grad_err(i),max_err(i)] = errors_exact_p2(p,t,p2,t2,basis,Qh,u,grad_u_x,grad_u_y); + + end + display_errors(err,grad_err,max_err) + +else + [p2,t2] = find_midpoints(p,t); + + [basis,Qh] = solve(p,p2,e,t,t2,f); + [err(1),grad_err(1),max_err(1)] = errors_exact_p2(p,t,p2,t2,basis,Qh,u,grad_u_x,grad_u_y); +end + +% end main +end + +% subfunction +function [basis,Qh] = solve(p,p2,e,t,t2,u) + basis = basis_functions_weighted_p2(p,t,p2,t2); + S = stiffness_matrix_weighted_p2(p,t,p2,t2,basis); + b = create_b_p2(p,t,p2,t2,basis,u); + Qh = S\b; + + figure(); + pdeplot([p,p2],e,t, 'XYData',Qh, 'ZData', Qh, 'Mesh', 'on'); +end diff --git a/main_no_exact.m b/main_no_exact.m new file mode 100644 index 0000000..5f4f90f --- /dev/null +++ b/main_no_exact.m @@ -0,0 +1,73 @@ +function [err,grad_err,max_err] = main(f,gd,sf,ns,mesh_level,opt) +%MAIN Summary of this function goes here +% Detailed explanation goes here +% +% Author: Nicole Stock +% Date: Fall 2020 + +model=createpde(1); +g=decsg(gd,sf,ns); +geometryFromEdges(model,g); +[p,e,t]=initmesh(g,'hmax',inf); +%pdemesh(p,e,t, 'NodeLabels','on', 'ElementLabels','on'); + +% Init error vectors +err = zeros(1,mesh_level); +grad_err = zeros(1,mesh_level); +max_err = zeros(1,mesh_level); + +if mesh_level > 1 + % To ensure we refine every triangle the same + [~,num_node]=size(p); + it=zeros(1,num_node); + for i=1:num_node + it(i)=i; + end + + [p2,t2] = find_midpoints(p,t); + + [basis,Qh] = solve(p,p2,e,t,t2,f); + % errors_no_exact_weighted_p2(p,t,p2,t2,basis,u_h_km1,u_h_k,k) + + for i = 2:mesh_level + basis_1 = basis; + Qh_1 = Qh; + t_1 = t; + t2_1 = t2; + [~,nodes_1] = size(p); + [~,mid_nodes_1] = size(p2); + + % Refine mesh to next level + [p,e,t]=refinemesh(g,p,e,t,it,'regular'); + + % Find the midpoints for P2 nodal points + [p2,t2] = find_midpoints(p,t); + + % Create prolongation matrix to extend mesh level - 1 vectors + P = prolongation_matrix(nodes_1,mid_nodes_1,p,p2,t_1,t2_1,basis_1,k); + + % Extend Qh + Qh_extended = P*Qh_1; + + % Solve + [basis,Qh] = solve(p,p2,e,t,t2,f); + [err(i),grad_err(i),max_err(i)] = errors_no_exact_weighted_p2(p,t,p2,t2,basis,Qh_extended,Qh,k); + end + + display_errors(err,grad_err,max_err) +end +% mesh level must be greater than 1 + +% end main +end + +% subfunction +function [basis,Qh] = solve(p,p2,e,t,t2,u) + basis = basis_functions_weighted_p2(p,t,p2,t2); + S = stiffness_matrix_weighted_p2(p,t,p2,t2,basis); + b = create_b_p2(p,t,p2,t2,basis,u); + Qh = S\b; + + figure(); + pdeplot([p,p2],e,t, 'XYData',Qh, 'ZData', Qh, 'Mesh', 'on'); +end diff --git a/main_no_exact.m~ b/main_no_exact.m~ new file mode 100644 index 0000000..61563e9 --- /dev/null +++ b/main_no_exact.m~ @@ -0,0 +1,78 @@ +function [err,grad_err,max_err] = main(f,gd,sf,ns,mesh_level,opt) +%MAIN Summary of this function goes here +% Detailed explanation goes here +% +% Author: Nicole Stock +% Date: Fall 2020 + +model=createpde(1); +g=decsg(gd,sf,ns); +geometryFromEdges(model,g); +[p,e,t]=initmesh(g,'hmax',inf); +%pdemesh(p,e,t, 'NodeLabels','on', 'ElementLabels','on'); + +% Init error vectors +err = zeros(1,mesh_level); +grad_err = zeros(1,mesh_level); +max_err = zeros(1,mesh_level); + +if mesh_level > 1 + % To ensure we refine every triangle the same + [~,num_node]=size(p); + it=zeros(1,num_node); + for i=1:num_node + it(i)=i; + end + + [p2,t2] = find_midpoints(p,t); + + [basis,Qh] = solve(p,p2,e,t,t2,f); + % errors_no_exact_weighted_p2(p,t,p2,t2,basis,u_h_km1,u_h_k,k) + + for i = 2:mesh_level + basis_1 = basis; + Qh_1 = Qh; + t_1 = t; + t2_1 = t2; + [~,nodes_1] = size(p); + [~,mid_nodes_1] = size(p2); + + % Refine mesh to next level + [p,e,t]=refinemesh(g,p,e,t,it,'regular'); + + % Find the midpoints for P2 nodal points + [p2,t2] = find_midpoints(p,t); + + % Create prolongation matrix to extend mesh level - 1 vectors + P = prolongation_matrix(nodes_1,mid_nodes_1,p,p2,t_1,t2_1,basis_1,k); + + % Extend Qh + Qh_extended = P*Qh_1; + + % Solve + [basis,Qh] = solve(p,p2,e,t,t2,f); + [err(i),grad_err(i),max_err(i)] = + errors_no_exact_weighted_p2(p,t,p2,t2,basis,Qh_extended,Qh,k); + end + display_errors(err,grad_err,max_err) + +else + [p2,t2] = find_midpoints(p,t); + + [basis,Qh] = solve(p,p2,e,t,t2,f); + [err(1),grad_err(1),max_err(1)] = errors_exact_p2(p,t,p2,t2,basis,Qh,u,grad_u_x,grad_u_y); +end + +% end main +end + +% subfunction +function [basis,Qh] = solve(p,p2,e,t,t2,u) + basis = basis_functions_weighted_p2(p,t,p2,t2); + S = stiffness_matrix_weighted_p2(p,t,p2,t2,basis); + b = create_b_p2(p,t,p2,t2,basis,u); + Qh = S\b; + + figure(); + pdeplot([p,p2],e,t, 'XYData',Qh, 'ZData', Qh, 'Mesh', 'on'); +end diff --git a/prolongation_matrix.m b/prolongation_matrix.m new file mode 100644 index 0000000..e49a9ed --- /dev/null +++ b/prolongation_matrix.m @@ -0,0 +1,70 @@ +function prolongation_matrix = prolongation_matrix(nodes_km1,mid_nodes_km1,p,p2,t_km1,t2_km1,basis_km1,k) +% PROLONGATION_MATRIX - Create prolongation matrix +% +% Syntax: +% prolongation_matrix = prolongation_matrix() +% +% Inputs: +% basis - a 3x3xNumTriangles matrix representing piece-wise basis +% functions for each node in each triangle. basis(i,:,k) represents +% the pieceiwise basis function for the ith node in triangle k. +% Outputs: +% prolongation_matrix - prolongation matrix +% +% Author: Nicole Stock +% Date: Spring 2020 + +[~,nodes] = size(p); +[~,mid_nodes] = size(p2); +total_nodes = nodes + mid_nodes; +total_nodes_km1 = nodes_km1 + mid_nodes_km1; + +% P is a total_nodes x total_nodes_km1 matrix +% ij entry of P : basis_km1_j evaluated at node i in p / p2 + +[~,triangles_km1] = size(t_km1); +i_vec = zeros(1,triangles_km1*6*total_nodes); +j_vec = zeros(1,triangles_km1*6*total_nodes); +s_vec = zeros(1,triangles_km1*6*total_nodes); +index = 1; + +for T = 1:triangles_km1 + for j = 1:6 + J = basis_km1(:,j,T); + if k == 0 + basis_fn_j =@(r,z) J(1).*r.^2 + J(2).*r.*z + J(3).*z.^2 ... + + J(4).*r + J(5).*z + J(6); + else + basis_fn_j =@(r,z) (r./k).*(J(1).*r.^2 + J(2).*r.*z ... + + J(3).*z.^2 + J(4).*r + J(5).*z + J(6)); + end + + if j >= 4 + global_j = t2_km1(j-3,T) + nodes_km1; + else + global_j = t_km1(j,T); + end + + for i = 1:nodes + ij_val = feval(basis_fn_j,p(1,i),p(2,i)); + + i_vec(index) = i; + j_vec(index) = global_j; + s_vec(index) = ij_val; + index = index + 1; + end + + for i = 1:mid_nodes + ij_val = feval(basis_fn_j,p2(1,i),p2(2,i)); + + i_vec(index) = i + nodes; + j_vec(index) = global_j; + s_vec(index) = ij_val; + index = index + 1; + end + end +end + +prolongation_matrix = sparse(i_vec,j_vec,s_vec,total_nodes,total_nodes_km1); + +% end \ No newline at end of file diff --git a/stiffness_matrix_weighted_p2.m b/stiffness_matrix_weighted_p2.m new file mode 100644 index 0000000..b068a29 --- /dev/null +++ b/stiffness_matrix_weighted_p2.m @@ -0,0 +1,131 @@ +function stiffness_matrix = stiffness_matrix_weighted_p2(p,t,p2,t2,basis,k) +% STIFFNESS_MATRIX_WEIGHTED_P2 - Create stiffness matrix with weight k +% +% Syntax: +% A = stiffness_matrix_weighted_p2(p,t,p2,t2,basis,k) +% +% Inputs: +% p - a 2xNumNodes matrix representing nodal coordinates. +% t - a 4xNumTriangles matrix representing the element connectivity in +% terms of node IDs. The end row of T represents the geometry face ID +% to which the element belongs +% basis - a 3x3xNumTriangles matrix representing piece-wise basis +% functions for each node in each triangle. basis(i,:,k) represents +% the pieceiwise basis function for the ith node in triangle k. +% k - given weight +% +% Outputs: +% stiffness_matrix - stiffness matrix +% +% Author: Nicole Stock +% Date: Spring 2020 + +[~,triangles] = size(t); +[~,nodes] = size(p); +i_vec = zeros(1,triangles*36); +j_vec = zeros(1,triangles*36); +s_vec = zeros(1,triangles*36); +index = 1; + +for T = 1:triangles + + % get coordinates of triangle T + coordinates = zeros(3,2); + for n = 1:3 + node = t(n,T); + % get x,y coordinates + coordinates(n,:) = p(:,node); + end + + [R,Z,Wr,Wz] = triquad(7, coordinates); + + % integrate for each pair of nodes in the triangle + for i = 1:6 + for j = i:6 + I = basis(:,i,T); + J = basis(:,j,T); + + % integrate grad(I) * grad(J) + I * J for (x,y) + if k == 0 + grad_i_r =@(r,z) 2.*I(1).*r + I(2).*z + I(4); + grad_i_z =@(r,z) I(2).*r + 2.*I(3).*z + I(5); + grad_j_r =@(r,z) 2.*J(1).*r + J(2).*z + J(4); + grad_j_z =@(r,z) J(2).*r + 2.*J(3).*z + J(5); + + grad_integrand =@(r,z) (grad_i_r(r,z).*grad_j_r(r,z) ... + + grad_i_z(r,z).*grad_j_z(r,z)).*r; + % integral grad(I) * grad(J) + Q1 = Wr'*feval(grad_integrand,R,Z)*Wz; + + integrand =@(r,z) (I(1).*r.^2 + I(2).*r.*z + I(3).*z.^2 ... + + I(4).*r + I(5).*z + I(6)).*(J(1).*r.^2 + J(2).*r.*z ... + + J(3).*z.^2 + J(4).*r + J(5).*z + J(6)).*r; + % integral I * J + Q2 = Wr'*feval(integrand,R,Z)*Wz; + else + % grad^k_rz(v) = [ partial_deriv_r(v) + % (-k/r)*v + % partial_deriv_z(v) ] + % phi_k_i = (r/k)(ar^2 + brz + cz^2 + dr + ez + f) + grad_i_r =@(r,z) (1./k).*(3.*I(1).*r.^2 + 2.*I(2).*r.*z ... + + I(3).*z.^2 + 2.*I(4).*r + I(5).*z + I(6)); + second_row_i =@(r,z) (I(1).*r.^2 + I(2).*r.*z ... + + I(3).*z.^2 + I(4).*r + I(5).*z + I(6)); + grad_i_z = @(r,z) (1./k).*(I(2).*r.^2 + 2.*I(3).*r.*z ... + + I(5).*r); + + grad_j_r =@(r,z) (1./k).*(3.*J(1).*r.^2 + 2.*J(2).*r.*z ... + + J(3).*z.^2 + 2.*J(4).*r + J(5).*z + J(6)); + second_row_j =@(r,z) (J(1).*r.^2 + J(2).*r.*z ... + + J(3).*z.^2 + J(4).*r + J(5).*z + J(6)); + grad_j_z =@(r,z) (1./k).*(J(2).*r.^2 + 2.*J(3).*r.*z ... + + I(5).*r); + + grad_integrand =@(r,z) (grad_i_r(r,z).*grad_j_r(r,z) ... + + second_row_i(r,z).*second_row_j(r,z) ... + + grad_i_z(r,z).*grad_j_z(r,z)).*r; + Q1 = Wr'*feval(grad_integrand,R,Z)*Wz; + + integrand =@(r,z) ((r./k).^2).*(I(1).*r.^2 + I(2).*r.*z ... + + I(3).*z.^2 + I(4).*r + I(5).*z + I(6)).*(J(1).*r.^2 ... + + J(2).*r.*z + J(3).*z.^2 + J(4).*r + J(5).*z ... + + J(6)).*r; + % integral I * J + Q2 = Wr'*feval(integrand,R,Z)*Wz; + end + + Q = Q1+Q2; + + % if i or j == 4,5,6, they are midpoint nodes, so their global + % id is found in t2 added to the number of original nodes + if i >= 4 + global_i = t2(i-3,T) + nodes; + else + global_i = t(i,T); + end + + if j >=4 + global_j = t2(j-3,T) + nodes; + else + global_j = t(j,T); + end + + i_vec(index) = global_i; + j_vec(index) = global_j; + s_vec(index) = Q; + index = index + 1; + if global_j ~= global_i + i_vec(index) = global_j; + j_vec(index) = global_i; + s_vec(index) = Q; + index = index + 1; + end + end + end +end + +[~,n2] = size(p2); +n = nodes + n2; +stiffness_matrix = sparse(i_vec,j_vec,s_vec,n,n); + +% end \ No newline at end of file diff --git a/triquad.m b/triquad.m new file mode 100644 index 0000000..1334b9e --- /dev/null +++ b/triquad.m @@ -0,0 +1,49 @@ +function [X,Y,Wx,Wy]=triquad(N,v) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% triquad.m - Gaussian Quadrature for a triangular domain +% +% This scripts computes the N^2 nodes and weights for a triangle with +% vertices given by the 3x2 vector v. The nodes are produced by collapsing +% the square to a triangle. +% +% Sample Usage: +% +% >>[X,Y,Wx,Wy]=triquad(8,[0 0; 0 2; 2 1]) +% >>f=@(x,y) exp(x+y); +% >>Q=Wx'*feval(f,X,Y)*Wy; +% +% Reference: J.N. Lyness, Ronald Cools, A Survey of Numerical Cubature +% over Triangles (1994) +% http://citeseer.ist.psu.edu/lyness94survey.html +% +% Written by: Greg von Winckel +% Contact: gregvw(at)math(dot)unm(dot)edu +% http://math.unm.edu/~gregvw +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +n=1:N; nnk=2*n+1; A=[1/3 repmat(1,1,N)./(nnk.*(nnk+2))]; +n=2:N; nnk=nnk(n); B1=2/9; nk=n+1; nnk2=nnk.*nnk; +B=4*(n.*nk).^2./(nnk2.*nnk2-nnk2); ab=[A' [2; B1; B']]; s=sqrt(ab(2:N,2)); +[V,X]=eig(diag(ab(1:N,1),0)+diag(s,-1)+diag(s,1)); +[X,I]=sort(diag(X)); x=(X+1)/2; wx=ab(1,2)*V(1,I)'.^2/4; + +N=N-1; N1=N+1; N2=N+2; y=cos((2*(N:-1:0)'+1)*pi/(2*N+2)); +L=zeros(N1,N2); y0=2; iter=0; +while max(abs(y-y0))>eps + L(:,1)=1; L(:,2)=y; + for k=2:N1 + L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k; + end + Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2); + y0=y; y=y0-L(:,N2)./Lp; iter=iter+1; +end + +cd=[ 1, 0, 0; -1, 0, 1; 0, 1,-1]*v; +t=(1+y)/2; Wx=abs(det(cd(2:3,:)))*wx; Wy=1./((1-y.^2).*Lp.^2)*(N2/N1)^2; +[tt,xx]=meshgrid(t,x); yy=tt.*xx; +X=cd(1,1)+cd(2,1)*xx+cd(3,1)*yy; Y=cd(1,2)+cd(2,2)*xx+cd(3,2)*yy; + +