fem1cub.m

来自「"introduction to the finite element meth」· M 代码 · 共 84 行

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%          FEM-1D Cubic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear % Clear all variables
%
% Define the Input Variables
% ==========================
L=8*10^-2; % Length of the domain
rho0=1*10^-8; % Charge density
epsr=1.0; % Dielectric constant of the domain
Va=1; % Boundary condition at the leftmost node
Vb=0; % Boundary condition at the rightmost node
Ne=3; % Number of quadratic elements
% ==========================
eps=epsr*8.85*10^-12;
Nn=3*Ne+1;
for e=1:Ne
    elmconn(e,1)=3*e-2;
    elmconn(e,2)=3*e;
    elmconn(e,3)=3*e+1;
    elmconn(e,4)=3*e-1;
end
le=L/Ne;

Ke=eps/(40*le)*[148 -13 -189 54; 
                -13 148 54 -189; 
                -189 54 432 -297;
                54 -189 -297 432];
fe=-le*rho0/8*[1;1;3;3];

K=zeros(Nn);
f=zeros(Nn,1);
for e=1:Ne
    for i=1:4
        for j=1:4
            K(elmconn(e,i),elmconn(e,j))=K(elmconn(e,i),elmconn(e,j))+Ke(i,j);
        end
        f(elmconn(e,i))=f(elmconn(e,i))+fe(i);
    end
end
for i=2:Nn
    f(i)=f(i)-K(i,1)*Va;
end
% K(:,1)=0;
% K(1,:)=0;
% K(1,1)=1;
% f(1)=Va;

for i=2:Nn-1
    f(i)=f(i)-K(i,Nn)*Vb;
end
% K(:,Nn)=0;
% K(Nn,:)=0;
% K(Nn,Nn)=1;
% f(Nn)=Vb;


V(2:Nn-1)=K(2:Nn-1,2:Nn-1)\f(2:Nn-1);

V(1)=Va;V(Nn)=Vb;
plot(V)