Utente:Kiban/Solution1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % solution.m % linear elastic analyses for plane problems %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all clear all

global displ analysis coor connec stress

disp 'reading input'; [val_dbc,val_tbc,tbc_number,pname,mnane]=read_input_analysis;  % read input data

[element_number,node_number,node_bc,node_force]=read_input_mesh(val_dbc,val_tbc,tbc_number,pname,mnane); %read input mesh %% Lettura delle proprietà del materiale E1=analysis.young1; E2=analysis.young2; nu1=analysis.poisson1; nu2=analysis.poisson2; G2=analysis.tangenziale2; beta=analysis.beta; type=analysis.type;

%% Calcolo la matrice di rigidezza disp 'Calcolo della matrice di rigidezza'; x=zeros(3,1); y=zeros(3,1); Ke=zeros(6,6); KK=zeros(2*node_number,2*node_number); FF=zeros(2*node_number,1);

for k=1:element_number,  % loop over the elements

   n=connec(k,:);                  % element nodes
   for i=1:3,
      x(i)=coor(n(i),1);           % coordinates of the nodes of the element
      y(i)=coor(n(i),2);
   end
   [Ke]=stiff_t3(x,y, E1, E2, nu1, nu2, G2, beta);             % computes element stiffness matrix
   
   for i=1:3,                      % matrix assemblage
      for j=1:3,
          KK(2*n(i)-1,2*n(j)-1)=KK(2*n(i)-1,2*n(j)-1)+Ke(2*i-1,2*j-1);
          KK(2*n(i),2*n(j))=KK(2*n(i),2*n(j))+Ke(2*i,2*j);
          KK(2*n(i)-1,2*n(j))=KK(2*n(i)-1,2*n(j))+Ke(2*i-1,2*j);
          KK(2*n(i),2*n(j)-1)=KK(2*n(i),2*n(j)-1)+Ke(2*i,2*j-1);
      end
   end

end

%%%% BC imposition disp 'computing BC condition and traction forces'; for j=1:node_number,

  if (node_bc(j,1)==1),
      KK(2*j-1,2*j-1)=1e10;
      FF(2*j-1)=1e10*node_bc(j,3);
  end
  if (node_bc(j,2)==1),
      KK(2*j,2*j)=1e10;
      FF(2*j)=1e10*node_bc(j,4);
  end

end

%%% computing nodal forces for j=1:node_number,

   if (node_force(j,1)==1),
      FF(2*j-1)= FF(2*j-1)+node_force(j,3);
  end
  if (node_force(j,2)==1),
      FF(2*j)=FF(2*j)+node_force(j,4);
  end
  

end

%% solving linear system disp 'solving linear system'; u=KK\FF;

%% compute stresses disp 'computing stress'; stress=stress_function(u,node_number,element_number, E1, E2, nu1, nu2, G2, beta, type);

for i=1:node_number,

  displ(i,1)=u(2*i-1); 
  displ(i,2)=u(2*i);

end

%%% results interface interface('Initialize',[' mesh-undef| mesh-def| mesh-both| u1|'...

         ' u2| s11| s22| s12'],'output_matlab(infopost)');