# 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 [Error_RF,DE] = ...
Mesh_Error(NOC,NE,UG,B0,Young,Poisson)
%
% SPECIFICATION: This function serves to calculate some parameter to
%                control the errors
%
% MAIN VARIABLES
%
% Error_RF:  Average rotational of F (deformation gradient) per element
%            It must be zero becase the rotational of a gradient is 0.
% NOC:       Nodal Coordinates
% NE:        Number of elements
% UG:        Unknowns in global coordinates
% B:         Bij is the partial of Ni respect to Xj (Voigt matrix)
% Young,Poisson: elastic parameters

% Initialize
Error_RF=0;
DE=0;

for N=1:NE

Node1 = NOC(N,1);
Node2 = NOC(N,2);
Node3 = NOC(N,3);
Node4 = NOC(N,4);

% Voigt matrix
BE = B0(1:4,1:2,N);

% Local unknowns
[UE] = Local(UG,Node1,Node2,Node3,Node4);

F1(1,1)=UE(3);   F1(1,2)=UE(4);  % F Node1
F1(2,1)=UE(5);   F1(2,2)=UE(6);
F2(1,1)=UE(9);   F2(1,2)=UE(10); % F Node2
F2(2,1)=UE(11);  F2(2,2)=UE(12);
F3(1,1)=UE(15);  F3(1,2)=UE(16); % F Node3
F3(2,1)=UE(17);  F3(2,2)=UE(18);
F4(1,1)=UE(21);  F4(1,2)=UE(22); % F Node4
F4(2,1)=UE(23);  F4(2,2)=UE(24);

F_Ave=(F1+F2+F3+F4)/4;

% 1st Piola-Kirchhoff stress tensor
Piola1 = Piola(Young,Poisson,F1);
Piola2 = Piola(Young,Poisson,F2);
Piola3 = Piola(Young,Poisson,F3);
Piola4 = Piola(Young,Poisson,F4);
Piola_Ave=(Piola1+Piola2+Piola3+Piola4)/4;

% Divergence of 1st P-K
[DP] = Div_P(Piola1,Piola2,Piola3,Piola4,BE);

% Linear momentum
p1 = UE(1:2,1);
p2 = UE(7:8,1);
p3 = UE(13:14,1);
p4 = UE(19:20,1);
p_Ave=(p1+p2+p3+p4)/4;

% Rotational error averaged
[ErrorE]   = Error_Rot_E(F1,F2,F3,F4,BE);
Error_RF  = Error_RF+(ErrorE/NE);

% Total derivative of entropy
DE        = DE + EntropyE;

end; % End of for N=1:NE

% End of function Error

% ----------- function Error_RotF_E --------------

function [ErrorE]=Error_Rot_E(F1,F2,F3,F4,B);
%
% SPECIFICATION: This Function serves to calculate the error
%                of the vector of Unknowns
% DESCRIPTION OF MAIN VARIABLES:
%
% INPUT
%
% Fi:        Deformation gradient of node i
%             | Fxx  Fxy |
%             | Fyx  Fyy |
% dt:        Increment of ime step
% BE:         Bij is the partial of Ni respect to Xj (Voigt matrix)
%
% OUTPUT
%
% RotF:     Rotational of the deformation gradient F
%              | Fxy,x - Fxx,y |
%              | Fyy,x - Fyx,y |

RotF=zeros(2,1);

DxFxy = F1(1,2)*B(1,1) + F2(1,2)*B(2,1) +...
F3(1,2)*B(3,1) + F4(1,2)*B(4,1);
DyFxx = F1(1,1)*B(1,2) + F2(1,1)*B(2,2) +...
F3(1,1)*B(3,2) + F4(1,1)*B(4,2);
DxFyy = F1(2,2)*B(1,1) + F2(2,2)*B(2,1) +...
F3(2,2)*B(3,1) + F4(2,2)*B(4,1);
DyFyx = F1(2,1)*B(1,2) + F2(2,1)*B(2,2) +...
F3(2,1)*B(3,2) + F4(2,1)*B(4,2);

RotF(1,1)=DxFxy-DyFxx;
RotF(2,1)=DxFyy-DyFyx;

ErrorE=sqrt(RotF(1)*RotF(1)+RotF(2)*RotF(2));

% End of function Error_Rot_F_E

%------------------------  function Local  ---------------------------
%
function [UE]=Local(UG,Node1,Node2,Node3,Node4);
%
% SPECIFICATION: This function serves to obtain
%                the local unknowns from global coordinates
%
% UG: Global vector
% UE: Local vector for the element

UE=zeros(24,1);

% Node 1
Start1 = 6*(Node1-1)+1;
End1   = Start1+5;
UE(1:6,1)   = UG(Start1:End1,1);

% Node 2
Start2 = 6*(Node2-1)+1;
End2   = Start2+5;
UE(7:12,1)  = UG(Start2:End2,1);

% Node 3
Start3 = 6*(Node3-1)+1;
End3   = Start3+5;
UE(13:18,1) = UG(Start3:End3,1);

% Node 4
Start4 = 6*(Node4-1)+1;
End4   = Start4+5;
UE(19:24,1) = UG(Start4:End4,1);

% End of function Local

%
% SPECIFICATION: This Fluxction serves to calculate the gradient
%                of the vector of Unknowns
%
% DESCRIPTION OF MAIN VARIABLES:
%
% INPUT
%
% pi=   linear momentum of node i
%            p=[ px
%                py ]
% dt:   Step of time
% B:    Bij is the partial of Ni respect to Xj (Voigt matrix)
%
% OUTPUT
%
%            | px,x    px,y |
%            | py,x    py,y |

% First column:  derivatives with respect to X

% Second column: derivatives with respect to Y

%------- function Div_P -----------------
%
function [DP] = Div_P(P1,P2,P3,P4,B)
%
% SPECIFICATION: This function serves to calculate the divergence
%                of the 1st Piola Kirchhoff stress tensor

DP=zeros(2,1);

DP=P1(1:2,1)*B(1,1)+P2(1:2,1)*B(2,1)+P3(1:2,1)*B(3,1)+P4(1:2,1)*B(4,1)+...
P1(1:2,2)*B(1,2)+P2(1:2,2)*B(2,2)+P3(1:2,2)*B(3,2)+P4(1:2,2)*B(4,2);

% End of Div_P function

%------- function Entropy_E -----------------
%