慕工大数值分析第五课第二题

clc;
clear;

%%  Mehrdimensionale Gauß-Quadratur

n = input('n = ');
gaussx = gx2dref(n);
gaussw = gw2dref(n);
s = 0;
nodes = [2, 1; 4, 1; 4, 3; 2, 2];

if n == 1
    for i = 1:n*n
        [J, det_J, inv_J] = getJacobian(nodes, gaussx(1, 1), gaussx(1, 2));
        N1 = 1/4*(1-gaussx(1,1))*(1-gaussx(1,2));
        N2 = 1/4*(1+gaussx(1,1))*(1-gaussx(1,2));
        s = N1*N2*det_J*gaussw;
    end

elseif (n == 2)
    for i = 1:n*n
        [J, det_J, inv_J] = getJacobian(nodes, gaussx(i, 1), gaussx(i, 2));
        N1=1/4*(1-gaussx(i,1))*(1-gaussx(i,2));
        N2=1/4*(1+gaussx(i,1))*(1-gaussx(i,2));
        s= s + N1*N2*det_J*gaussw(i);
    end

elseif n==3
    for i=1:n*n
        [J,det_J,invJ] = getJacobian([2, 1; 4, 1; 4, 3; 2, 2], gaussx(i,1), gaussx(i,2));
        N1=1/4*(1-gaussx(i,1))*(1-gaussx(i,2));
        N2=1/4*(1+gaussx(i,1))*(1-gaussx(i,2));
        s=s + N1*N2*det_J*gaussw(i);
    end
end
display(s)

function gaussw = gw(n)

if n == 1

gaussw = 2;

elseif n == 2

gaussw = [1, 1];

elseif n == 3

gaussw = [5/9, 8/9, 5/9];

else

disp('false n')

end

function gaussx = gx(n)

if n == 1

gaussx = 0;

elseif n == 2

gaussx = [-1/sqrt(3), 1/sqrt(3)];

elseif n == 3

gaussx = [-sqrt(3/5), 0, sqrt(3/5)];

else

disp('false n')

end

function gaussx = gx2dref(n)
if n == 1
    gaussx = [0, 0];
elseif n == 2
    gaussx = [-1/sqrt(3), -1/sqrt(3); -1/sqrt(3), 1/sqrt(3); 1/sqrt(3), -1/sqrt(3); 1/sqrt(3), 1/sqrt(3)];
elseif n == 3
    k = sqrt(3/5);
    gaussx = [-k, -k; -k, 0; -k, k; 0, -k; 0, 0; 0, k; k, -k; k, 0; k, k];
else
    disp('false n')
end

function gaussw = gw2dref(n)
w_1 = 5/9;
w_2 = 8/9;
if n == 1
    gaussw = [4.0];
elseif n == 2
    gaussw = [1.0; 1.0; 1.0; 1.0];
elseif n == 3
    gaussw = [w_1*w_1; w_1*w_2; w_1*w_1; w_1*w_2; w_2*w_2; w_1*w_2; w_1*w_1; w_1*w_2; w_1*w_1];
else
    disp('false n')
end

function x = getxPos(nodes, x_i, eta)

x = linquadref(x_i, eta)' * nodes;

end

function [J, det_J, inv_J] = getJacobian(nodes, x_i, eta)
d_J = linquadderiref(x_i, eta)' * nodes;
J = d_J';
det_J = det(J);
inv_J = inv(J);
end

function deriv = linquadderiref(x_i, eta)

deriv(1, 1) = -1/4 * (1 - eta);

deriv(1, 2) = -1/4 * (1 - x_i);

deriv(2, 1) = 1/4 * (1 - eta);

deriv(2, 2) = -1/4 * (1 + x_i);

deriv(3, 1) = 1/4 * (1 + eta);

deriv(3, 2) = 1/4 * (1 + x_i);

deriv(4, 1) = -1/4 * (1 + eta);

deriv(4, 2) = 1/4 * (1 - x_i);

end

function val = linquadref(x_i, eta)

val(1, 1) = 1/4 * (1 - x_i) * (1 - eta);

val(2, 1) = 1/4 * (1 + x_i) * (1 - eta);

val(3, 1) = 1/4 * (1 + x_i) * (1 + eta);

val(4, 1) = 1/4 * (1 - x_i) * (1 + eta);

end