dyn_fem.m

FEM tools for caculation of nonlinear problems

function [xc, t] = dyn_FEM ( in_data, obj, L, LK) 

% =========================================================================
% FEM analysis toolbox for solid mechanics. Project started: Aug 2004
% Anton Zaicenco <[e-mail]>
% =========================================================================

global dof h CH CH2 C K zdd delta_t

dimn = length(L);

[dof_per_node] = type_of_elem (in_data);   

dof = size(in_data.ND,1)*dof_per_node;     
hh  = zeros(dof(1),1);                     

for i=1:length(in_data.dynam.TIMEHM)    
    t4 = in_data.dynam.TIMEHM(i)*dof_per_node-(dof_per_node-1);
    for k=1:dof_per_node
        hh(t4+(k-1))=1;  
    end;
end;

for i=1:length(obj.M)        
   if obj.M(i,i)>0
      if hh(i)==1  hh(i)= 1;  else hh(i)= 0; end;
   end;
end;
 
for i=1:dof_per_node:dof(1)
    for k=1:dof_per_node
        if hh(i+(k-1))~=0  & in_data.dynam.TIMEHDIR(k)==0  
            hh(i+(k-1))=0;   end; 
        if hh(i+(k-1))~=0  & in_data.dynam.TIMEHDIR(k)~=0  
            hh(i+(k-1))=hh(i+(k-1))*in_data.dynam.TIMEHDIR(k); end;
    end;
end;

if  isfield(in_data,'slaves')
    if length(in_data.slaves)>0
        [u,T] = MFC_master_slave (obj.Ksys.Kgl,L,123,in_data.slaves,in_data.master,0);
        obj.M_old = obj.M;
        clear u
        obj.M = T'*obj.M*T;
        K = T'*obj.Ksys.Kgl*T;
        hh = T'*hh;
        h = hh;
    end
else
    h=hh(L);
end

clear hh
zdd     = load(in_data.dynam.TIMEH);
points  = length(zdd);
delta_t = in_data.dynam.delta_tm;

disp(['     Dynamic analysis - start Cholesky factorization. Time (sec): ' num2str(toc) ]);

if  isfield(in_data,'slaves')
    if length(in_data.slaves)>0
        CH = chol(obj.M);
    end
else
    CH = chol(obj.M(L,L));
end

disp(['     Dynamic analysis - Cholesky done. Time (sec): ' num2str(toc) ]);
CH2 = CH';

switch 2
    case 1
        
        damping = in_data.dynam.DAMP_C;
        N=in_data.dynam.DAMP_F;
        if  isfield(in_data,'slaves')
            if length(in_data.slaves)>0
                [alp, bet] = rayleigh_coef (obj.Ksys.Kgl,obj.M,damping,N);
            end
        else
            [alp, bet] = rayleigh_coef (obj.Ksys.Kgl(L,L),obj.M(L,L),damping,N);
        end
    case 2

        alp = in_data.dynam.ab(1); bet = in_data.dynam.ab(2);
end

disp(['  Rayleigh coef.: alpha =  ' num2str(alp) ' .  beta = ' num2str(bet)]);

if  isfield(in_data,'slaves') 
    if length(in_data.slaves)>0
        C = alp*obj.M + bet*K; 
    end
else
    K = obj.Ksys.Kgl(L,L);
    C = alp*obj.M(L,L) + bet*obj.Ksys.Kgl(L,L);
end

% =========================================================================
clear i k t4 damping
dof = length(C);
x0 = zeros(2*dof,1);
t      = [1:points-1]*delta_t;
disp(['     Dynamic analysis - START. Time (sec): ' num2str(toc) ]);

[t,xc] = odeRK6('dxdtU', t, x0, zdd', 1.0e-9);

disp(['     Dynamic analysis - END.   Time (sec): ' num2str(toc) ])
clear CH  CH2 h zdd x0 K
% =========================================================================
disp(['     Dynamic analysis : recovering displ., veloc., and accel' ])

sd   = length(t);
h=[1:LK]; h(L)=0;

if  isfield(in_data,'slaves')
    if length(in_data.slaves)>0
        xc_d = T*xc(:,1:size(T,2))';
        xc_v = T*xc(:,size(T,2)+1:end)';
        xc   = [xc_d' xc_v'];
        dof = length(L);  
        C = alp*obj.M_old(L,L) + bet*obj.Ksys.Kgl(L,L);
    end
else
    xc   = [ xc   zeros(sd,2*(LK-dof)) ];
    xc(:,LK+L) = xc(:,dof+1:2*dof);            
    xc(:,L)    = xc(:,1:dof);                  
    xc(:,find(h))=0; xc(:,LK+find(h))=0;
end

disp(['     Dynamic analysis : velocities and displacements recovered' ])

xc   = [ xc   zeros(sd,LK) ];

xc(:,LK*2+1:LK*2+dof) = (C*(xc(:,LK+L))')';

clear C
disp(['     Dynamic analysis : C erased' ])
if  isfield(in_data,'slaves')
    if length(in_data.slaves)>0
        xc(:,LK*2+1:LK*2+dof) = xc(:,LK*2+1:LK*2+dof) + (obj.Ksys.Kgl(L,L)*(xc(:,L))')';
    end
else
    xc(:,LK*2+1:LK*2+dof) = xc(:,LK*2+1:LK*2+dof) + (obj.Ksys.Kgl(L,L)*(xc(:,L))')';
end

clear K
disp(['     Dynamic analysis : K erased' ])
if  isfield(in_data,'slaves')
    if length(in_data.slaves)>0
        xc(:,LK*2+1:LK*2+dof) = (-obj.M_old(L,L)\xc(:,LK*2+1:LK*2+dof)')';
    end
else
    xc(:,LK*2+1:LK*2+dof) = (-obj.M(L,L)\xc(:,LK*2+1:LK*2+dof)')';
end

clear MM M
disp(['     Dynamic analysis : MM erased' ])
xc(:,2*LK+L) = xc(:,LK*2+1:LK*2+dof);

xc(:,2*LK+find(h))=0;

clear h

time_now = toc;
time_min = floor(time_now/60); time_hrs = floor(time_min/60);
time_min = floor(time_min-time_hrs*60);
time_sec = floor(time_now-time_hrs*60*60-time_min*60);

disp([     '  ...   Dynamic analysis: OK.        Time:  ' ...
    num2str(time_hrs) ' hrs : ' num2str(time_min) ' min : ' ...
    num2str(time_sec) ' sec ' ])
% =========================================================================