FEM-Course-Matlab/9.四面体单元悬臂梁Matlab有限元编程/main_for_C3D4.m

clear all;clc;close all;
E=210000; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><EFBFBD>
u=0.3;    %<EFBFBD><EFBFBD><EFBFBD>ɱ<EFBFBD>
D=E/((1+u)*(1-2*u))*[1-u u u 0 0 0;......
    u 1-u u 0 0 0;......
    u u 1-u 0 0 0;......
    0 0 0 (1-2*u)/2 0 0;......
    0 0 0 0 (1-2*u)/2 0;......
    0 0 0 0 0 (1-2*u)/2];%elastic modules matrix
load nodes_and_elements.mat
Forces=[7 2 -100;10 2 -100;12 2 -100;13 2 -100;14 2 -100;]; %[nodeID direction value]
ConID=[2,   3,   5,   8,  36,  37,  38,  61,  62,  63,  64,  87,  88,  89,  90,  91......
  92,  93, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260];
Constraints=zeros(size(ConID,2)*3,3);
for i=1:size(ConID,2)
Constraints(3*i-2:3*i,:)=[ConID(i) 1 0;ConID(i) 2 0;ConID(i) 3 0;];%[nodeID direction value];
end
Dof=3; %ÿ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶ<EFBFBD>
NodeNum = size(Nodes,1);  % <EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
ElementNum= size(Elements,1);  %<EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Dofs = Dof*NodeNum; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶ<EFBFBD><EFBFBD><EFBFBD>
U=sparse(Dofs,1);	% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD>ṹλ<EFBFBD>ƣ<EFBFBD><EFBFBD>վ<EFBFBD><EFBFBD>󣩣<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>߼<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ٶ<EFBFBD>
K = sparse(Dofs,Dofs);  %<EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ն<EFBFBD><EFBFBD>󣨿վ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
F = sparse(Dofs,1);  %<EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>վ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>

for I=1:ElementNum
    % <EFBFBD><EFBFBD>Ԫ<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    ElementNodeCoordinate=Nodes(Elements(I,:),:);
    % <EFBFBD><EFBFBD><EFBFBD>㵥Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
    ElementStiffnessMatrix=Ke_C3D4(D,ElementNodeCoordinate);
    % <EFBFBD><EFBFBD><EFBFBD>㵥Ԫ<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶȱ<EFBFBD><EFBFBD><EFBFBD>
    ElementNodeDOF=zeros(1,12);
    for J=1:4
        II=(J-1)*Dof+1;
        ElementNodeDOF(II:II+2)=(Elements(I,J)-1)*Dof+1:(Elements(I,J)-1)*Dof+3;
    end
    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>װ
    K(ElementNodeDOF,ElementNodeDOF)=K(ElementNodeDOF,ElementNodeDOF)+ElementStiffnessMatrix;%<EFBFBD><EFBFBD><EFBFBD>ݵ<EFBFBD>Ԫ<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶ<EFBFBD>  <EFBFBD><EFBFBD>װ<EFBFBD>ܸ<EFBFBD>
end
% ʩ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
if size(Forces,1)>0
    ForceDOF = Dof*(Forces(:,1)-1)+Forces(:,2);   %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶȱ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD>Ϊ<EFBFBD>ڵ<EFBFBD><EFBFBD>ţ<EFBFBD><EFBFBD>ڶ<EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    F(ForceDOF) = F(ForceDOF) + Forces(:,3);
end
% <EFBFBD>˴<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʩ<EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD>
BigNumber=1e8;
ConsNum=size(Constraints,1);
if ConsNum~=0
    FixedDof=Dof*(Constraints(:,1)-1)+Constraints(:,2);  %<EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɶȱ<EFBFBD><EFBFBD><EFBFBD>(<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>)
    for i=1:ConsNum
        K(FixedDof(i),FixedDof(i))=K(FixedDof(i),FixedDof(i))*BigNumber;
        F(FixedDof(i))=Constraints(i,3)*K(FixedDof(i),FixedDof(i));
    end
end
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>
U = K\F;
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>˹<EFBFBD>㡢<EFBFBD>ڵ<EFBFBD>Ӧ<EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD>ֵ
[NodeStrain,NodeStress]=NodeStressStrain(U,D,Nodes,Elements);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>MISESӦ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
ElementNodeCount=4;% ÿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>
MISES=zeros(1,ElementNum*ElementNodeCount);
for I=1:ElementNum*ElementNodeCount
    MISES(I)=sqrt(0.5)*sqrt((NodeStress(1,I)-NodeStress(2,I))^2+(NodeStress(1,I)-NodeStress(3,I))^2+....
        (NodeStress(2,I)-NodeStress(3,I))^2+6*(NodeStress(4,I)^2+NodeStress(5,I)^2+NodeStress(6,I)^2));
end

newNodes=Nodes';
newNodes=newNodes(:);
SF=5.0e2; %<EFBFBD><EFBFBD><EFBFBD>ηŴ<EFBFBD>ϵ<EFBFBD><EFBFBD>
newNodes=newNodes+SF*U;
newNodes=reshape(newNodes,[3,size(Nodes,1)]);
newNodes=newNodes';
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ԭʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͼ
  for i=1:1:size(Elements,1)
    points=Nodes(Elements(i,:),:);
    deformed_points=newNodes(Elements(i,:),:);;
    mesh=1:1:4;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϣ
    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>嵥Ԫ<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    vertices_matrix = [points(mesh(1,:),1),points(mesh(1,:),2),points(mesh(1,:),3)];
    deformed_vertices_matrix=[deformed_points(mesh(1,:),1),deformed_points(mesh(1,:),2),deformed_points(mesh(1,:),3)];
    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>嵥Ԫ<EFBFBD>ڵ<EFBFBD>˳<EFBFBD><EFBFBD>
    faces_matrix= [1 2 4;2 3 4;3 1 4;1 3 2];%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ÿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ľڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    patch('vertices', vertices_matrix,'faces',faces_matrix,'facecolor','g','FaceAlpha',.5);
    view(3);hold on%<EFBFBD><EFBFBD>ͼ
    patch('vertices', deformed_vertices_matrix,'faces',faces_matrix,'facecolor','r','FaceAlpha',.5);
    
  end 
axis equal
view(3);
alpha(1);


%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Umag<EFBFBD><EFBFBD>ͼ
Umag=zeros(NodeNum,1);
for i=1:NodeNum
    Umag(i)=sqrt(U(3*i-2)^2+U(3*i-1)^2+U(3*i)^2);
end
PostContour(Nodes,Elements,U,Umag)
title('Umag')
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>MISES<EFBFBD><EFBFBD>ͼ
PostContour(Nodes,Elements,U,MISES)
title('Mises')