www.pudn.com > 87361026FEA.rar > lagrange_basis.m, change:2002-10-01,size:10967b

```function [Nv,dNdxi]=lagrange_basis(type,coord,dim)

% returns the lagrange interpolant basis and its gradients w.r.t the
% parent coordinate system.
%
%         [N(xi),dNdxi(xi)]=lagrange_basis(type-order,coord,dim)
%
%   type is the toplogical class of finite element it is in the general
%   form 'topology-#of nodes' ie a three node triangel is T3 a four
%   node quadralateral is Q4 a 4 node tetrahedra is H4 a 27 node brick
%   is B27 etc
%
%   coord is the parent coordinates at which the basis and its
%   gradients are to be evaluated at.
%
%   presently defined are L2, L3, T3, T4(cubic bubble), T6, Q4, Q9,
%   H4, H10, B8 and B27
%
%   If dim is set to 2 then the vector representation of the N
%   matrix is returned.
%
% written by Jack Chessa
%            j-chessa@northwestern.edu
% Department of Mechanical Engineering
% Northwestern University

if ( nargin == 2 )
dim=1;
end

switch type
case 'L2'
%%%%%%%%%%%%%%%%%%%%% L2 TWO NODE LINE ELEMENT %%%%%%%%%%%%%%%%%%%%%
%
%    1---------2
%
if size(coord,2) < 1
disp('Error coordinate needed for the L2 element')
else
xi=coord(1);
N=([1-xi,1+xi]/2)';
dNdxi=[-1;1]/2;
end

case 'L3'
%%%%%%%%%%%%%%%%%%% L3 THREE NODE LINE ELEMENT %%%%%%%%%%%%%%%%%%%%%
%
%    1---------2----------3
%
if size(coord,2) < 1
disp('Error two coordinates needed for the L3 element')
else
xi=coord(1);
N=[(1-xi)*xi/(-2);(1+xi)*xi/2;1-xi^2];
dNdxi=[xi-.5;xi+.5;-2*xi];
end

case 'T3'
%%%%%%%%%%%%%%%% T3 THREE NODE TRIANGULAR ELEMENT %%%%%%%%%%%%%%%%%%
%
%               3
%             /  \
%            /    \
%           /      \
%          /        \
%         /          \
%        /            \
%       /              \
%      /                \
%     /                  \
%    1--------------------2
%
if size(coord,2) < 2
disp('Error two coordinates needed for the T3 element')
else
xi=coord(1); eta=coord(2);
N=[1-xi-eta;xi;eta];
dNdxi=[-1,-1;1,0;0,1];
end

case 'T3fs'
if size(coord,2) < 2
disp('Error two coordinates needed for the T3fs element')
else
xi=coord(1); eta=coord(2);
N=[1-xi-eta;xi;eta];
dNdxi=[-1,-1;1,0;0,1];
end

case 'T4'
%%%%%%%%%% T4 FOUR NODE TRIANGULAR CUBIC BUBBLE ELEMENT %%%%%%%%%%%%
%
%               3
%             /  \
%            /    \
%           /      \
%          /        \
%         /          \
%        /      4     \
%       /              \
%      /                \
%     /                  \
%    1--------------------2
%
if size(coord,2) < 2
disp('Error two coordinates needed for the T4 element')
else
xi=coord(1); eta=coord(2);
N=[1-xi-eta-3*xi*eta;xi*(1-3*eta);eta*(1-3*xi);9*xi*eta];
dNdxi=[-1-3*eta,-1-3*xi;
1-3*eta, -3*xi;
-3*eta,   1-3*xi;
9*eta,   9*xi ];
end

case 'T6'
%%%%%%%%%%%%%%%%%% T6 SIX NODE TRIANGULAR ELEMENT %%%%%%%%%%%%%%%%%%
%
%               3
%             /  \
%            /    \
%           /      \
%          /        \
%         6          5
%        /            \
%       /              \
%      /                \
%     /                  \
%    1---------4----------2
%
if size(coord,2) < 2
disp('Error two coordinates needed for the T6 element')
else
xi=coord(1); eta=coord(2);
N=[1-3*(xi+eta)+4*xi*eta+2*(xi^2+eta^2);
xi*(2*xi-1);
eta*(2*eta-1);
4*xi*(1-xi-eta);
4*xi*eta;
4*eta*(1-xi-eta)];

dNdxi=[4*(xi+eta)-3   4*(xi+eta)-3;
4*xi-1              0;
0        4*eta-1;
4*(1-eta-2*xi)          -4*xi;
4*eta           4*xi;
-4*eta  4*(1-xi-2*eta)];
end

case 'Q4'
%%%%%%%%%%%%%%% Q4 FOUR NODE QUADRILATERIAL ELEMENT %%%%%%%%%%%%%%%%
%
%    4--------------------3
%    |                    |
%    |                    |
%    |                    |
%    |                    |
%    |                    |
%    |                    |
%    |                    |
%    |                    |
%    |                    |
%    1--------------------2
%
if size(coord,2) < 2
disp('Error two coordinates needed for the Q4 element')
else
xi=coord(1); eta=coord(2);
N=1/4*[ (1-xi)*(1-eta);
(1+xi)*(1-eta);
(1+xi)*(1+eta);
(1-xi)*(1+eta)];
dNdxi=1/4*[-(1-eta), -(1-xi);
1-eta,    -(1+xi);
1+eta,      1+xi;
-(1+eta),   1-xi];
end

case 'Q9'
%%%%%%%%%%%%%%% Q9 NINE NODE QUADRILATERIAL ELEMENT %%%%%%%%%%%%%%%%
%
%    4---------7----------3
%    |                    |
%    |                    |
%    |                    |
%    |                    |
%    8          9         6
%    |                    |
%    |                    |
%    |                    |
%    |                    |
%    1----------5---------2
%
if size(coord,2) < 2
disp('Error two coordinates needed for the Q9 element')
else
xi=coord(1); eta=coord(2);
N=1/4*[xi*eta*(xi-1)*(eta-1);
xi*eta*(xi+1)*(eta-1);
xi*eta*(xi+1)*(eta+1);
xi*eta*(xi-1)*(eta+1);
-2*eta*(xi+1)*(xi-1)*(eta-1);
-2*xi*(xi+1)*(eta+1)*(eta-1);
-2*eta*(xi+1)*(xi-1)*(eta+1);
-2*xi*(xi-1)*(eta+1)*(eta-1);
4*(xi+1)*(xi-1)*(eta+1)*(eta-1)];
dNdxi=1/4*[eta*(2*xi-1)*(eta-1),xi*(xi-1)*(2*eta-1);
eta*(2*xi+1)*(eta-1),xi*(xi+1)*(2*eta-1);
eta*(2*xi+1)*(eta+1),xi*(xi+1)*(2*eta+1);
eta*(2*xi-1)*(eta+1),xi*(xi-1)*(2*eta+1);
-4*xi*eta*(eta-1),   -2*(xi+1)*(xi-1)*(2*eta-1);
-2*(2*xi+1)*(eta+1)*(eta-1),-4*xi*eta*(xi+1);
-4*xi*eta*(eta+1),   -2*(xi+1)*(xi-1)*(2*eta+1);
-2*(2*xi-1)*(eta+1)*(eta-1),-4*xi*eta*(xi-1);
8*xi*(eta^2-1),      8*eta*(xi^2-1)];
end

case 'H4'
%%%%%%%%%%%%%%%% H4 FOUR NODE TETRAHEDRAL ELEMENT %%%%%%%%%%%%%%%%%%
%
%             4
%           / | \
%          /  |  \
%         /   |   \
%        /    |    \
%       /     |     \
%      1 -----|------3
%         -   2  -
if size(coord,2) < 3
disp('Error three coordinates needed for the H4 element')
else
xi=coord(1); eta=coord(2); zeta=coord(3);
N=[1-xi-eta-zeta;
xi;
eta;
zeta];
dNdxi=[-1  -1  -1;
1   0   0;
0   1   0;
0   0   1];
end

case 'H10'
%%%%%%%%%%%%%%%% H10 TEN NODE TETRAHEDRAL ELEMENT %%%%%%%%%%%%%%%%%%
disp(['Element ',type,' not yet supported'])
if size(coord,2) < 3
disp('Error three coordinates needed for the H10 element')
else
xi=coord(1); eta=coord(2); zeta=coord(3);
N=zeros(10,1);
dNdxi=zeros(10,3);
end

case 'B8'
%%%%%%%%%%%%%%%%%%% B8 EIGHT NODE BRICK ELEMENT %%%%%%%%%%%%%%%%%%%%
%
%                  8
%               /    \
%            /          \
%         /                \
%      5                     \
%      |\                     7
%      |   \                / |
%      |     \     4    /     |
%      |        \    /        |
%      |           6          |
%      1           |          |
%       \          |          3
%          \       |        /
%            \     |     /
%               \  |  /
%                  2
%
if size(coord,2) < 3
disp('Error three coordinates needed for the B8 element')
else
xi=coord(1); eta=coord(2); zeta=coord(3);
I1=1/2-coord/2;
I2=1/2+coord/2;
N=[   I1(1)*I1(2)*I1(3);
I2(1)*I1(2)*I1(3);
I2(1)*I2(2)*I1(3);
I1(1)*I2(2)*I1(3);
I1(1)*I1(2)*I2(3);
I2(1)*I1(2)*I2(3);
I2(1)*I2(2)*I2(3);
I1(1)*I2(2)*I2(3)   ];
dNdxi=[   -1+eta+zeta-eta*zeta   -1+xi+zeta-xi*zeta  -1+xi+eta-xi*eta;
1-eta-zeta+eta*zeta   -1-xi+zeta+xi*zeta  -1-xi+eta+xi*eta;
1+eta-zeta-eta*zeta    1+xi-zeta-xi*zeta  -1-xi-eta-xi*eta;
-1-eta+zeta+eta*zeta    1-xi-zeta+xi*zeta  -1+xi-eta+xi*eta;
-1+eta-zeta+eta*zeta   -1+xi-zeta+xi*zeta   1-xi-eta+xi*eta;
1-eta+zeta-eta*zeta   -1-xi-zeta-xi*zeta   1+xi-eta-xi*eta;
1+eta+zeta+eta*zeta    1+xi+zeta+xi*zeta   1+xi+eta+xi*eta;
-1-eta-zeta-eta*zeta    1-xi+zeta-xi*zeta   1-xi+eta-xi*eta  ]/8;
end

case 'B27'
%%%%%%%%%%%%%% B27 TWENTY SEVEN NODE BRICK ELEMENT %%%%%%%%%%%%%%%%%
disp(['Element ',type,' not yet supported'])
if size(coord,2) < 3
disp('Error three coordinates needed for the B27 element')
else
xi=coord(1); eta=coord(2); zeta=coord(3);
N=zeros(27,1);
dNdxi=zeros(27,3);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
otherwise
disp(['Element ',type,' not yet supported'])
N=[]; dNdxi=[];
end

I=eye(dim);
Nv=[];
for i=1:size(N,1)
Nv=[Nv;I*N(i)];
end

if ( dim == 1 )
B=dNdxi;
elseif ( dim == 2 )
B=zeros(dim*size(N,1),3);

B(1:dim:dim*size(N,1)-1,1) = dNdxi(:,1);
B(2:dim:dim*size(N,1),2)   = dNdxi(:,2);

B(1:dim:dim*size(N,1)-1,3) = dNdxi(:,2);
B(2:dim:dim*size(N,1),3)   = dNdxi(:,1);
elseif ( dim == 3 )
B=zeros(dim*size(N,1),6);

disp('Error: need to add 3D N and dNdxi')

B(1:dim:dim*size(N,1)-2,1) = dNdxi(:,1);
B(2:dim:dim*size(N,1)-1,2) = dNdxi(:,2);
B(3:dim:dim*size(N,1),3)   = dNdxi(:,3);

B(2:dim:dim*size(N,1)-1,4) = dNdxi(:,3);
B(3:dim:dim*size(N,1),4)   = dNdxi(:,2);

B(3:dim:dim*size(N,1),5)   = dNdxi(:,1);
B(1:dim:dim*size(N,1)-2,5) = dNdxi(:,3);

B(1:dim:dim*size(N,1)-2,6) = dNdxi(:,2);
B(2:dim:dim*size(N,1)-1,6) = dNdxi(:,1);

end

% end of function```