Skip to content

Commit

Permalink
First commit
Browse files Browse the repository at this point in the history
  • Loading branch information
@NEStock
NEStock committed Sep 8, 2020
0 parents commit 92574b5
Show file tree
Hide file tree
Showing 12 changed files with 1,036 additions and 0 deletions.
63 changes: 63 additions & 0 deletions basis_functions_weighted_p2.m
Original file line number Diff line number Diff line change
@@ -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
67 changes: 67 additions & 0 deletions create_b_p2.m
Original file line number Diff line number Diff line change
@@ -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);
145 changes: 145 additions & 0 deletions errors_exact_weighted_p2.m
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit 92574b5

Please sign in to comment.