www.pudn.com > 87361026FEA.rar > static.m, change:2006-02-28,size:7168b


% bimaterial_enriched.m
%
% by Jack Chessa
% Northwestern University
%
% A small displacement static code for analysisng a bimaterial plate
% We use an discontinously enriched assumed strain finite element.
clear
colordef black
state = 0;

%******************************************************************************
%***                             I N P U T                                  ***
%******************************************************************************
tic;
disp('************************************************')
disp('***          S T A R T I N G    R U N        ***')
disp('************************************************')

disp([num2str(toc),'   START'])
% HERE WE DEFINE SOME SYSTEM PARPMETERS
homeDir='/home/jack';   % home directory
addpath(homeDir,[homeDir,'/MatlabFEMCode']);
addpath(homeDir,[homeDir,'/MatlabFEMCode/MeshGeneration']);

% GMSH FILE PARPMETERS
gmshFileName = 'plate.msh';  % name of gmsh input file
domainID     = 9;            % zone id for domain mesh
topEdgeID    = 8;            % zone id for displacement boundary mesh
botEdgeID    = 7;            % zone id for traction boundary mesh
sideEdgeID   = 10;           % zone id for traction boundary mesh

% MATERIAL PROPERTIES
lambda = 10;    % Lame constants
mu     = 3;   
lambda = 1;  
Ht      = 5.333; % y-coordinate of material interface
Height  = 10;    % plate hieght
Width   = 5;

% BOUNDARY CONDITIONS
trac_y=10;
%******************************************************************************
%***                    P R E - P R O C E S S I N G                         ***
%******************************************************************************

disp([num2str(toc),'   READING INPUT FILE'])
% READ GMSH FILE
[node,elements,elemType]=msh2mlab(gmshFileName);

% CHECK FOR ONLY TRIANGULAR ELEMENTS

% GET NODE AND CONNECTIVITY MATRICIES
node=node(:,1:2);                % a nodal coordinate matix (reference config)
element  = elements{domainID};   % connectivity matrix for interior
topEdge = elements{topEdgeID};   % connectivity matrix for top boundary
botEdge = elements{botEdgeID};   % connectivity matrix for bottom boundary
sides   = elements{sideEdgeID};  % connectivity matrix for side boundary

% CHECK FOR BAD JACOBIANS
disp([num2str(toc),'   CHECKING MESH FILE'])
tricheck(node,element);

% % STRUCTURED MESH
% numx=6;
% numy=2*numx;
% [node,elements]=make_cross_mesh([0 0],[Width Height],numx,numy);
% topEdge = elements{3};   % connectivity matrix for top boundary
% botEdge = elements{1};   % connectivity matrix for bottom boundary
% sides   = [elements{2};
%            elements{4}]; % connectivity matrix for side boundary
% element = elements{5};   % connectivity matrix for interior

% GET NODES ON DISPLACEMENT BOUNDARY
botEdgeNodes=unique(botEdge);
sideEdgeNodes=unique(sides);

numnode=size(node,1);    % number of nodes
numelem=size(element,1); % number of elements

% DEFINE SYSTEM DATA STRUCTURES
disp([num2str(toc),'   INITIALIZING DATA STRUCTURES'])
U=zeros(4*numnode,1);          % nodal displacement vector
X=node;                        % nodal coordinates in current configureation
f=zeros(2*numnode,1);          % external load vector
K=sparse(2*numnode,2*numnode); % stiffness matrix

xs=1:numnode;                  % x portion of u and v vectors
ys=(numnode+1):2*numnode;      % y portion of u and v vectors

%******************************************************************************
%***                          P R O C E S S I N G                           ***
%******************************************************************************

ind2voit=[1;4;2];   % the are used to switch between indicie and voigt form
voit2ind=[1 3;3 2];

disp([num2str(toc),'   COMPUTING LOAD'])
% COMPUTE EXTERNAL FORCE 
[W,Q]=quadrature( 1, 'GAUSS', 1 );  % define quadrature rule
for e=1:size(topEdge,1)
  
  sctr=topEdge(e,:);
  sctrx=sctr;
  sctry=sctrx+numnode;
  
  for q=1:size(W,1)                      % quadrature loop
    pt=Q(q,:);                           % quadrature point
    wt=W(q);                             % quadrature weight
    N=lagrange_basis('L2',pt);           % element shape functions
    J0=abs( X(sctr(2))-X(sctr(1)) )/2;
    
    f(sctry)=f(sctry)+N*trac_y*det(J0)*wt;
    
  end % of quadrature loop
end  % of element loop


%%%%%%%%%%%%%%%%%%%%% COMPUTE STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp([num2str(toc),'   COMPUTING STIFFNESS MATRIX'])
for e=1:numelem                          % start of element loop
  
  sctr=element(e,:);
  sctrB0=[sctr(1)         sctr(1)+numnode sctr(2) ...
          sctr(2)+numnode sctr(3)         sctr(3)+numnode ];
  
  [W,Q]=quadrature(1,'TRIANGULAR',2);      % define discontionus quadrature rule
  for q=1:size(W,1)                        % quadrature loop
    pt=Q(q,:);                             % quadrature point
    wt=W(q);                               % quadrature weight
    [N,dNdxi]=lagrange_basis('T3',pt);     % element shape functions
    
    J0=X(sctr,:)'*dNdxi;  % element Jacobian matrix w.r.t ref configuration       
    invJ0=inv(J0);
    dNdX=dNdxi*invJ0;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % COMPUTE B BAR MATRIX
    B0=[dNdX(1,1)         0 dNdX(2,1)         0 dNdX(3,1)         0;
                0 dNdX(1,2)         0 dNdX(2,2)         0 dNdX(3,2);
        dNdX(1,2) dNdX(1,1) dNdX(2,2) dNdX(2,1) dNdX(3,2) dNdX(3,1)];
      
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % COMPUTE COMPLIANCE MATRIX
    C=[lambda+2*mu  lambda      0;
       lambda       lambda+2*mu 0;
       0            0           mu];
     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % COMPUTE ELEMENT STIFFNESS AT QUADRATURE POINT
    K(sctrB0,sctrB0)=K(sctrB0,sctrB0)+B0'*C*B0*W(q)*det(J0);
    
  end                                    % of quadrature loop
end   % of element loop
%%%%%%%%%%%%%%%%%%% END OF STIFFNESS MATRIX COMPUTATION %%%%%%%%%%%%%%%%%%%%%%

% APPLY BOUNDARY CONDITIONS ( only for zero essential bcs )
disp([num2str(toc),'   APPLYING BOUNDARY CONDITIONS'])
bcwt=mean(diag(K));
ydofs=botEdgeNodes+numnode;
xdofs=sideEdgeNodes;     
K(xdofs,:)=0;
K(ydofs,:)=0;
K(:,xdofs)=0;
K(:,ydofs)=0;
K(xdofs,xdofs)=bcwt*speye(size(xdofs,1));
K(ydofs,ydofs)=bcwt*speye(size(ydofs,1));
f(xdofs)=0;
f(ydofs)=0;

% SOLVE SYSTEM
disp([num2str(toc),'   SOLVING SYSTEM'])
U=K\f;

%******************************************************************************
%***                     P O S T  -  P R O C E S S I N G                    ***
%******************************************************************************
disp([num2str(toc),'   POST-PROCESSING'])

trisurf(element,X(:,1),X(:,2),U(ys),U(ys));

disp([num2str(toc),'   RUN FINISHED'])
% ***************************************************************************
% ***                     E N D   O F    P R O G R A M                    ***
% ***************************************************************************
disp('************************************************')
disp('***            E N D    O F    R U N         ***')
disp('************************************************')