www.gusucode.com > pde 案例源码 matlab代码程序 > pde/cCoefficientLagrangePlaneStress.m

    function c = cCoefficientLagrangePlaneStress(E, nu, loc, state)
%cCoefficientLagrangePlaneStress Calculate c-coefficient for nonlinear plane stress
% Calculate the c-coefficient for a geometrically nonlinear Lagrangian formulation 
% of plane stress elasticity. The strain measure is the Green-Lagrange strain
% tensor. The stress is the second Piola-Kirchoff stress tensor. The material
% is assumed to be isotropic with linear behavior (Hooke's law applies).
% 
% E  - Young's modulus of the linear isotropic material
% nu - Poisson's ratio for the material
% p  - matrix of point (node) locations
% t  - element connectivity matrix
% u  - current displacement vector

%    This function was generated by the Symbolic Math Toolbox version 6.0.
%    31-Jan-2014 09:50:09

% Copyright 2014-2015 The MathWorks, Inc.


ux = reshape(state.ux,2,[]);
uy = reshape(state.uy,2,[]);



dudx=ux(1,:); 
dvdx=ux(2,:);
dudy=uy(1,:);  
dvdy=uy(2,:);
 
% if(~isempty(u))
%   [ux,uy] = pdegrad(p,t,u);
%   dudx=ux(1,:); dudy=uy(1,:); dvdx=ux(2,:); dvdy=uy(2,:);
% else
%   dudx = zeros(1, size(t,2)); dudy=dudx; dvdx=dudx; dvdy=dudx;
% end

t4 = 1/(nu^2-1);
t6 = 1/(1+nu);
t7 = E*dudy.*t4*.25;
t8 = dudx+1.0;
t9 = E*dudy.*t4.*t8*.25;
t10 = dvdy+1.0;
t11 = t7+t9-E*dvdx.*t6.*t10*.25;
t12 = dvdy.*2.0;
t13 = dudx.^2;
t14 = dudy.^2;
t15 = dvdy.^2;
t16 = dvdx.^2;
t17 = E*dvdx.*t4.*(1.0./2.0);
t18 = E*dudx.*dvdx.*t4*.25;
t19 = t17+t18-E*dudy.*t6*.25-E*dudy.*dvdy.*t6.*(1.0./8.0);
t20 = E*dudy.*dvdx.*nu.*t4*.25;
t21 = t20-E*t6.*(1.0./2.0)-E*dudx.*t6*.25-E*dvdy.*t6*.25-E*dudx.*dvdy.*t6.*(1.0./8.0);
t22 = dudx.*2.0;
t23 = dvdy+2.0;
t24 = nu-1.0;
t25 = E*nu.*t4;
t26 = E*dudx.*nu.*t4.*(1.0./2.0);
t27 = E*dvdy.*nu.*t4.*(1.0./2.0);
t28 = E*dudx.*dvdy.*nu.*t4*.25;
t29 = t25+t26+t27+t28-E*dudy.*dvdx.*t6.*(1.0./8.0);
t30 = E*dudy.*t4.*t23.*(1.0./8.0);
t31 = E*dudy.*dvdy.*t4.*(1.0./8.0);
t32 = t7+t30+t31-E*dvdx.*t6.*(1.0./8.0)-E*dvdx.*t4.*t8.*t24.*(1.0./8.0);
t33 = dudy.*2.0;
t34 = dvdx.*2.0;
t35 = dudx.*dudy.*2.0;
t36 = dvdx.*dvdy;
t37 = t33+t34+t35+t36;
t38 = 1.0./t24;
t39 = E*dvdx.*t23.*t38.*(1.0./8.0);
t40 = t39-E*t6.*t37.*(1.0./8.0);
out1 = [E*t4.*(dudx.*6.0+t13.*2.0+t14+t16+4.0)*.25+E*nu.*t4.*(t12+t15)*.25;
  t11;
  t19;
  t29;
  t11;
  E*t4.*(t12+t13+t14.*2.0+t15+t22+2.0)*.25+E*nu.*t4.*(t16-2.0)*.25;
  t21;
  t32;
  t19;
  t21;
  E*t4.*(t12+t13+t15+t16.*2.0+t22+2.0)*.25+E*nu.*t4.*(t14-2.0)*.25;
  t40;
  t29;
  t32;
  t40;
  E*t4.*(dvdy.*6.0+t14+t15.*2.0+t16+4.0)*.25+E*nu.*t4.*(t13+t22)*.25];

c = -out1([1 5 6 9 10 13 14 11 15 16], :);

end