Code covered by the BSD License  

Highlights from
FEM toolbox for solid mechanics

image thumbnail
from FEM toolbox for solid mechanics by Anton Zaicenco
The finite element toolbox for solid mechanics with GUI.

fem3D_tetrah_cyl (in_data)
function in_data = fem3D_tetrah_cyl  (in_data)

t=load('tetrah_2.txt');
p=load('points_2.txt');

cnt = 1;
d=ones(length(p)*length(p),1);
for j=1:length(p)
            d(j)=sqrt( (0-p(j,1))^2 + (0.5-p(j,2))^2 + (1-p(j,3))^2   );
end


pl = length(p);
tl = size(t);

in_data.ND(1:pl,1) = [1:pl];
in_data.ND(1:pl,2) = p(:,1);
in_data.ND(1:pl,3) = p(:,2);
in_data.ND(1:pl,4) = p(:,3);


in_data.mater.E    = 7e3;
in_data.mater.miu  = 0.3;
in_data.mater.rho  = 2000; % density of material - X

E = in_data.mater.E; rho = in_data.mater.rho; miu = in_data.mater.miu;

in_data.EL(1:tl(1),1) = [1:tl(1)];
in_data.EL(1:tl(1),2) = 10;
in_data.EL(1:tl(1),3:6) = t;
in_data.EL(1:tl(1),7) = E;
in_data.EL(1:tl(1),8) = miu;
in_data.EL(1:tl(1),9) = rho;

in_data.mater.E=0; in_data.mater.rho=0; in_data.mater.miu=0;

for i=1:length(in_data.EL)
    node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
    node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
    node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
    node4 = find(in_data.ND(:,1)==in_data.EL(i,6));
    x1=in_data.ND(node1,2); y1=in_data.ND(node1,3);  z1=in_data.ND(node1,4);
    x2=in_data.ND(node2,2); y2=in_data.ND(node2,3);  z2=in_data.ND(node2,4);
    x3=in_data.ND(node3,2); y3=in_data.ND(node3,3);  z3=in_data.ND(node3,4);
    x4=in_data.ND(node4,2); y4=in_data.ND(node4,3);  z4=in_data.ND(node4,4);
    if det([1 1 1 1; x1 x2 x3 x4; y1 y2 y3 y4; z1 z2 z3 z4])<0

        ss=in_data.EL(i,5); in_data.EL(i,5)=in_data.EL(i,4); in_data.EL(i,4)=ss; 
    end;
end;
disp(['-------- negative DET computation ... done']);
for i=1:length(in_data.EL)
    node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
    node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
    node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
    node4 = find(in_data.ND(:,1)==in_data.EL(i,6));
    x1=in_data.ND(node1,2); y1=in_data.ND(node1,3);  z1=in_data.ND(node1,4);
    x2=in_data.ND(node2,2); y2=in_data.ND(node2,3);  z2=in_data.ND(node2,4);
    x3=in_data.ND(node3,2); y3=in_data.ND(node3,3);  z3=in_data.ND(node3,4);
    x4=in_data.ND(node4,2); y4=in_data.ND(node4,3);  z4=in_data.ND(node4,4);
    if det([1 1 1 1; x1 x2 x3 x4; y1 y2 y3 y4; z1 z2 z3 z4])<0
        disp(['negative determinant:  ' num2str(i)]);
    end;
end



indx = 0;
for i=1:pl
    if in_data.ND(i,4)==-1;
        indx=indx+1;
        in_data.CON(indx,1)=in_data.ND(i,1);
        in_data.CON(indx,2)=0;
        in_data.CON(indx,3)=0;
        in_data.CON(indx,4)=0;
    end
end

indx = 0;
for i=1:pl
    if in_data.ND(i,4)>(-.2) & in_data.ND(i,4)<(.2)
        if in_data.ND(i,2)<(-.8)
        indx=indx+1;
        in_data.LOAD_(indx,1)=in_data.ND(i,1);
        in_data.LOAD_(indx,2)=100;
        in_data.LOAD_(indx,3)=0;
        in_data.LOAD_(indx,4)=0;
        end
    end
end



in_data.dynam.TIMEH    = [ 'pulse_a2.txt' ];
in_data.dynam.delta_tm = [0.005];       
in_data.dynam.TIMEHDIR = [0 0 -1];      
in_data.dynam.TIMEHM   = [941];%1:length(in_data.EL)*3]; 
in_data.dynam.TIMEHPL  = [150*3-0];         
in_data.dynam.DAMP_C   = [0.02 0.02];       
in_data.dynam.DAMP_F   = [3];               
in_data.dynam.ab       = [0.4 0.0001];      
in_data.dynam.MODA     = [1];               

Contact us at files@mathworks.com