81 lines
2.3 KiB
Mathematica
81 lines
2.3 KiB
Mathematica
|
function [] = FESolve()
|
||
|
|
% MATLAB based FE Solver
|
||
|
|
% J. Grogan (2012)
|
||
|
|
clear all
|
||
|
|
% Define Geometry
|
||
|
|
len=1.;
|
||
|
|
% Define Section Properties
|
||
|
|
rho=1.;
|
||
|
|
spec=1.;
|
||
|
|
cond=1.;
|
||
|
|
% Generate Mesh
|
||
|
|
numElem=10;
|
||
|
|
ndCoords=linspace(0,len,numElem+1);
|
||
|
|
numNodes=size(ndCoords,2);
|
||
|
|
indx=1:numElem;
|
||
|
|
elemNodes(:,1)=indx;
|
||
|
|
elemNodes(:,2)=indx+1;
|
||
|
|
% Initialize Conditions
|
||
|
|
Tnew=zeros(numNodes,1);
|
||
|
|
Tnew(1)=1.;
|
||
|
|
% Define Time Step
|
||
|
|
dtime=.1;
|
||
|
|
tsteps=10;
|
||
|
|
time=0.;
|
||
|
|
% Loop through time steps
|
||
|
|
for ts=1:tsteps;
|
||
|
|
time=time+dtime;
|
||
|
|
K=zeros(numNodes,numNodes);
|
||
|
|
M=zeros(numNodes,numNodes);
|
||
|
|
% Loop Through Elements
|
||
|
|
for e=1:numElem;
|
||
|
|
Ke=zeros(2);
|
||
|
|
Me=zeros(2);
|
||
|
|
gpx(1)=-1./sqrt(3.);
|
||
|
|
gpx(2)=1./sqrt(3.);
|
||
|
|
ajacob=abs(ndCoords(elemNodes(e,2))-ndCoords(elemNodes(e,1)))/2.;
|
||
|
|
% Loop Through Int Points
|
||
|
|
for i=1:2;
|
||
|
|
c=gpx(i);
|
||
|
|
phi(1)=(1.-c)/2.;
|
||
|
|
phi(2)=(1.+c)/2.;
|
||
|
|
phic(1)=-0.5;
|
||
|
|
phic(2)=0.5;
|
||
|
|
phix(1)=phic(1)/ajacob;
|
||
|
|
phix(2)=phic(2)/ajacob;
|
||
|
|
we=ajacob;
|
||
|
|
Ke=Ke+we*cond*phix'*phix;
|
||
|
|
Me=Me+(we*rho*spec*phi'*phi)/dtime;
|
||
|
|
end
|
||
|
|
% Assemble Global Matrices
|
||
|
|
K(elemNodes(e,1),elemNodes(e,1))=K(elemNodes(e,1),elemNodes(e,1))+Ke(1,1);
|
||
|
|
K(elemNodes(e,1),elemNodes(e,2))=K(elemNodes(e,1),elemNodes(e,2))+Ke(1,2);
|
||
|
|
K(elemNodes(e,2),elemNodes(e,1))=K(elemNodes(e,2),elemNodes(e,1))+Ke(2,1);
|
||
|
|
K(elemNodes(e,2),elemNodes(e,2))=K(elemNodes(e,2),elemNodes(e,2))+Ke(2,2);
|
||
|
|
M(elemNodes(e,1),elemNodes(e,1))=M(elemNodes(e,1),elemNodes(e,1))+Me(1,1);
|
||
|
|
M(elemNodes(e,1),elemNodes(e,2))=M(elemNodes(e,1),elemNodes(e,2))+Me(1,2);
|
||
|
|
M(elemNodes(e,2),elemNodes(e,1))=M(elemNodes(e,2),elemNodes(e,1))+Me(2,1);
|
||
|
|
M(elemNodes(e,2),elemNodes(e,2))=M(elemNodes(e,2),elemNodes(e,2))+Me(2,2);
|
||
|
|
end
|
||
|
|
%Apply Boundary Conditions (Reduce Matrices)
|
||
|
|
T1=1;
|
||
|
|
RHS=M*Tnew;
|
||
|
|
for i=1:numNodes-1;
|
||
|
|
for j=1:numNodes-1;
|
||
|
|
Kred(i,j)=K(i+1,j+1);
|
||
|
|
Mred(i,j)=M(i+1,j+1);
|
||
|
|
end
|
||
|
|
Subr(i)=(K(i+1,1)+M(i+1,1))*T1;
|
||
|
|
RHSr(i)=RHS(i+1);
|
||
|
|
end
|
||
|
|
%Solve
|
||
|
|
StiffI=(Mred+Kred)^-1;
|
||
|
|
Tnewr=StiffI*(RHSr'-Subr');
|
||
|
|
for i=1:numNodes-1;
|
||
|
|
Tnew(i+1)=Tnewr(i);
|
||
|
|
end
|
||
|
|
Tnew(1)=1.;
|
||
|
|
Tnew
|
||
|
|
end
|
||
|
|
|