# Two step Taylor Galerkin algorithm

### Xose Manuel Carreira (view profile)

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

```