fem3d_tetrah_cyl.m

FEM tools for caculation of nonlinear problems

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];