No BSD License  

Highlights from
Two step Taylor Galerkin algorithm

image thumbnail

Two step Taylor Galerkin algorithm

by

 

18 Jun 2007 (Updated )

The two step Taylor Galerkin algorithm applied to nonlinear structural dynamics

...
function [DJ0,DJ1,DJ2,DJ3,DJ4,B0,B1,B2,B3,B4,Normal,Length,N1,N2] = ...
          Geometry(NE,X,NOC,NBC,NBE,NB1,NB2)
%
% SPECIFICATION: This function serves to create the geometric 
%                properties of the element
%
% MAIN VARIABLES
%
% OUTPUT
% DJ:        Jacobian of the local-global transformation DJ=2A 
%            if counter clock wise notation of nodes DJ>0
% X:         Coordinates
% NOC:       Number of components
% N:         Number of the element
% Normal:    Normal to the side on the boundary element
% Length:    Length of the side on the boundary of the element
% N1,N2:     Local coodinates of the nodes on teh boundary NB1,NB2 
%
% INPUT
% NE:        Number of elements
% X:         Global coordinates of nodes
% NOC:       Local numeration of nodes
% NBC:       Number of Boundary Conditions
% NBE:       Number of Element on the boundary
% NB1,NB2:   Nodes on the Boundary side

DJ0 = zeros(1,NE);
DJ1 = zeros(1,NE);
DJ2 = zeros(1,NE);
DJ3 = zeros(1,NE);
DJ4 = zeros(1,NE);

B0  = zeros(4,2,NE);
B1  = zeros(4,2,NE);
B2  = zeros(4,2,NE);
B3  = zeros(4,2,NE);
B4  = zeros(4,2,NE);

Normal = zeros(2,1,NBC);
Length = zeros(1,NBC);
N1 = zeros(1,NBC);
N2 = zeros(1,NBC);

XI0=0.0;            ETA0=0.0; 
XI1=-1.0/sqrt(3.0); ETA1=-1.0/sqrt(3.0); 
XI2= 1.0/sqrt(3.0); ETA2=-1.0/sqrt(3.0); 
XI3= 1.0/sqrt(3.0); ETA3= 1.0/sqrt(3.0); 
XI4=-1.0/sqrt(3.0); ETA4= 1.0/sqrt(3.0); 

for N = 1:NE
      
   [DJ0(1,N),B0(1:4,1:2,N)] = Voigt(N,X,NOC,XI0,ETA0);
   [DJ1(1,N),B1(1:4,1:2,N)] = Voigt(N,X,NOC,XI1,ETA1);
   [DJ2(1,N),B2(1:4,1:2,N)] = Voigt(N,X,NOC,XI2,ETA2);
   [DJ3(1,N),B3(1:4,1:2,N)] = Voigt(N,X,NOC,XI3,ETA3);
   [DJ4(1,N),B4(1:4,1:2,N)] = Voigt(N,X,NOC,XI4,ETA4);   
   
end;

for side = 1:NBC
   
   [Normal(1:2,1,side),Length(1,side),N1(1,side),N2(1,side)]=...
      Get_Normal(X,NOC,NB1(side),NB2(side),NBE(side));

end;

% End of function Geometry



%------------------------  function Voigt  ---------------------------
      
function [DJ,B] = Voigt(N,X,NOC,XI,ETA);

% SPECIFICATION: This procedure serves to create the Jacobian and Voigt B matrix

% MAIN VARIABLES

% DJ:        Jacobian of the local-global transformation DJ=2A 
%            if counter clock wise notation of nodes DJ>0
% B(i,j):    dNi/dXj 
% X:         Coordinates
% NOC:       Number of components
% N:         Number of the element
% I1,I2,I3:  Global notation of node

%--- Get coordinates

      I1 = NOC(N, 1);
      I2 = NOC(N, 2);
      I3 = NOC(N, 3);
      I4 = NOC(N, 4);      
      
      X1 = X(I1, 1);
      Y1 = X(I1, 2);
      X2 = X(I2, 1);
      Y2 = X(I2, 2);
      X3 = X(I3, 1);
      Y3 = X(I3, 2);
      X4 = X(I4, 1);
      Y4 = X(I4, 2);

%  --- Determinant of the JACOBIAN

      J11 = ((1 - ETA) * (X2 - X1) + (1 + ETA) * (X3 - X4)) / 4;
      J12 = ((1 - ETA) * (Y2 - Y1) + (1 + ETA) * (Y3 - Y4)) / 4;
      J21 = ((1 - XI) * (X4 - X1) + (1 + XI) * (X3 - X2)) / 4;
      J22 = ((1 - XI) * (Y4 - Y1) + (1 + XI) * (Y3 - Y2)) / 4;      

      DJ = J11 * J22 - J12 * J21;

%--- Definition of B() Matrix  (Voigt matrix)
%     B(i,j)= partial of shape function Ni respect to Xj

%  --- A(3,4) Matrix relates Strains to
%  --- Local Derivatives of u
      A = zeros(3, 4);

      A(1, 1) = J22 / DJ;
      A(2, 1) = 0;
      A(3, 1) = -J21 / DJ;
      A(1, 2) = -J12 / DJ;
      A(2, 2) = 0;
      A(3, 2) = J11 / DJ;
      A(1, 3) = 0;
      A(2, 3) = -J21 / DJ;
      A(3, 3) = J22 / DJ;
      A(1, 4) = 0;
      A(2, 4) = J11 / DJ;
      A(3, 4) = -J12 / DJ;
      
%  --- G(4,8) Matrix relates Local Derivatives of u
%  --- to Local Nodal Displacements q(8)
      G = zeros(4, 8);
      
      G(1, 1) = -(1 - ETA) / 4;
      G(2, 1) = -(1 - XI) / 4;
      G(3, 2) = -(1 - ETA) / 4;
      G(4, 2) = -(1 - XI) / 4;
      G(1, 3) = (1 - ETA) / 4;
      G(2, 3) = -(1 + XI) / 4;
      G(3, 4) = (1 - ETA) / 4;
      G(4, 4) = -(1 + XI) / 4;
      G(1, 5) = (1 + ETA) / 4;
      G(2, 5) = (1 + XI) / 4;
      G(3, 6) = (1 + ETA) / 4;
      G(4, 6) = (1 + XI) / 4;
      G(1, 7) = -(1 + ETA) / 4;
      G(2, 7) = (1 - XI) / 4;
      G(3, 8) = -(1 + ETA) / 4;
      G(4, 8) = (1 - XI) / 4;
      
      AG=zeros(3,8);
      AG=A*G;

      B(1,1)=AG(1,1);
      B(2,1)=AG(1,3);
      B(3,1)=AG(1,5);
      B(4,1)=AG(1,7);
      
      B(1,2)=AG(2,2);
      B(2,2)=AG(2,4);
      B(3,2)=AG(2,6);
      B(4,2)=AG(2,8);
             
      
% End of function Voigt   
   
   
%------------------------ function Get_Normal  ---------------------------
%
function [Normal,Length,N1,N2]=Get_Normal(X,NOC,NB1,NB2,N);
%
% SPECIFICATION: This function serves to calculate:
%                - The normal between NB1 NB2 for the element N,
%                - The length of the side 
%                - The local coordinates N1,N2 of NB1,NB2

         if     (NB1==NOC(N,1)); N1=1; 
         elseif (NB1==NOC(N,2)); N1=2; 
         elseif (NB1==NOC(N,3)); N1=3; 
         elseif (NB1==NOC(N,4)); N1=4;             
         end;
         
         if     (NB2==NOC(N,1)); N2=1;
         elseif (NB2==NOC(N,2)); N2=2; 
         elseif (NB2==NOC(N,3)); N2=3; 
         elseif (NB2==NOC(N,4)); N2=4; 
         end;
         
         if (N1==N2)            
            Normal=zeros(2,1);                   
            Length=0;         
         else   
                 
         X1 = X(NB1,1);
         Y1 = X(NB1,2);
         X2 = X(NB2,1);
         Y2 = X(NB2,2);
         
         X12 = X1 - X2;
         Y12 = Y1 - Y2;
                 
         Length=sqrt(X12*X12+Y12*Y12);
         
         Normal(1,1)= -Y12/Length; 
         Normal(2,1)=  X12/Length;  
                                             
        end;  % End of if 
        
% End of function Get_Normal     

Contact us