Getting a warning while running the program.Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.394922e-018..

2 次查看（过去 30 天）
Deviprakash Upadhyaya
I am running the below given program which gives the warning message.How to remove the warning. ---------------------------------------------------------------------
A barrel vault has radius r=25 ft, subtended angle of 80 degrees,
length 50 f, and thickness 3in.Elastic modulus is 3 msi,Poisson's ratio is 0,and weight of the shell is 90 lb per square ft. Two curved edges are
supported by rigid diaphragm and the other two edges are free.
Solve using 4 by 4 elements of a quarter of the shell.
Variable descriptions
% k = element matrix in the local axes
% ke = element matrix in the global axes
% kb = element matrix for bending stiffness
% ks = element matrix for shear stiffness
% km = element matrix for membrane stiffness
% f = element vector
% kk = system matrix
% ff = system vector
% disp = system nodal displacement vector
% gcoord = coordinate values of each node
% nodes = nodal connectivity of each element
% index = a vector containing system dofs associated with each element
% pointb = matrix containing sampling points for bending term
% weightb = matrix containing weighting coefficients for bending term
% points = matrix containing sampling points for shear term
% weights = matrix containing weighting coefficients for shear term
% bcdof = a vector containing dofs associated with boundary conditions
% bcval = a vector containing boundary condition values associated with
% the dofs in 'bcdof'
% kinmtsb = matrix for kinematic equation for bending
% matmtsb = matrix for material property for bending
% kinmtss = matrix for kinematic equation for shear
% matmtss = matrix for material property for shear
% kinmtsm = matrix for kinematic equation for membrane
% matmtsm = matrix for material property for membrane
% tr3d = transformation matrix from local to global axes
%-------------------------- % input data for control parameters %---------------------------
clear
nel=16; % number of elements
nnel=4; % number of nodes per element
ndof=6; % number of dofs per node
nnode=25; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
edof=nnel*ndof; % degrees of freedom per element
emodule=3e6; % elastic modulus
poisson=0.0; % Poisson's ratio
t=3; % plate thickness
nglxb=2; nglyb=2; % 2x2 Gauss-Legendre quadrature for bending
nglb=nglxb*nglyb; % number of sampling points per element for bending
nglxs=1; nglys=1; % 1x1 Gauss-Legendre quadrature for shear
ngls=nglxs*nglys; % number of sampling points per element for shear
%--------------------------------- % input data for nodal coordinate values %---------------------------------
gcoord=[0.0 0.0 0.0; 0.0 6.25 0.0; 0.0 12.5 0.0; 0.0 18.75 0.0;0.0 25.0 0.0;
4.34 0.0 -0.38;4.34 6.25 -0.38;4.34 12.5 -0.38;4.34 18.75 -0.38;4.34 25.0 -0.38;
8.55 0.0 -1.51;8.55 6.25 -1.51;8.55 12.5 -1.51;8.55 18.75 -1.51;8.55 25.0 -1.51;
12.5 0.0 -3.35;12.5 6.25 -3.35;12.5 12.5 -3.35;12.5 18.75 -3.35;12.5 25.0 -3.35;
16.1 0.0 -5.85;16.1 6.25 -5.85;16.1 12.5 -5.85;16.1 18.75 -5.85;16.1 25.0 -5.85];
gcoord=12.0*gcoord;
%----------------------- % input data for nodal connectivity for each element %-----------------------
nodes=[6 7 2 1; 7 8 3 2; 8 9 4 3; 9 10 5 4;
11 12 7 6; 12 13 8 7; 13 14 9 8; 14 15 10 9;
16 17 12 11; 17 18 13 12; 18 19 14 13; 19 20 15 14;
21 22 17 16; 22 23 18 17; 23 24 19 18; 24 25 20 19];
%------------------------------------- % input data for boundary conditions %-------------------------------------
bcdof=[1 3 5 6 7 11 12 13 17 18 19 23 24 25 26 28 29 30 ...
31 33 56 58 60 61 63 86 88 90 91 93 116 118 120 ...
121 122 123 146 148 150]; % constrained dofs
bcval=zeros(size(bcdof)); % whose described values are 0
%-------------------------------------- % initialization of matrices and vectors %------------------------------------
ff=zeros(sdof,1); % system force vector
kk=zeros(sdof,sdof); % system matrix
disp=zeros(sdof,1); % system displacement vector
index=zeros(edof,1); % index vector
kinmtsb=zeros(3,edof); % kinematic matrix for bending
matmtsb=zeros(3,3); % constitutive matrix for bending
kinmtsm=zeros(3,edof); % kinematic matrix for membrane
matmtsm=zeros(3,3); % constitutive matrix for membrane
kinmtss=zeros(2,edof); % kinematic matrix for shear
matmtss=zeros(2,2); % constitutive matrix for shear
tr3d=zeros(edof,edof); % transformation matrix
%---------------------------- % force vector %----------------------------
ff(3)=-613.6; % transverse force at node 1
ff(9)=-1227.2; % transverse force at node 2
ff(15)=-1227.2; % transverse force at node 3
ff(21)=-1227.2; % transverse force at node 4
ff(27)=-613.6; % transverse force at node 5
ff(33)=-1227.2; % transverse force at node 6
ff(39)=-2454.4; % transverse force at node 7
ff(45)=-2454.4; % transverse force at node 8
ff(51)=-2454.4; % transverse force at node 9
ff(57)=-1227.2; % transverse force at node 10
ff(63)=-1227.2; % transverse force at node 11
ff(69)=-2454.4; % transverse force at node 12
ff(75)=-2454.4; % transverse force at node 13
ff(81)=-2454.4; % transverse force at node 14
ff(87)=-1227.2; % transverse force at node 15
ff(93)=-1227.2; % transverse force at node 16
ff(99)=-2454.4; % transverse force at node 17
ff(105)=-2454.4; % transverse force at node 18
ff(111)=-2454.4; % transverse force at node 19
ff(117)=-1227.2; % transverse force at node 20
(123)=-613.6; % transverse force at node 21
ff(129)=-1227.2; % transverse force at node 22
ff(135)=-1227.2; % transverse force at node 23
ff(141)=-1227.2; % transverse force at node 24
ff(147)=-613.6; % transverse force at node 25
%-------------------- % computation of element matrices and vectors and their assembly %---------------------
% for bending and membrane stiffness
[pointb,weightb]=feglqd2(nglxb,nglyb); % sampling points & weights
matmtsm=fematiso(1,emodule,poisson)*t; % membrane material property
matmtsb=fematiso(1,emodule,poisson)*t^3/12; % bending material property
[points,weights]=feglqd2(nglxs,nglys); % sampling points & weights
shearm=0.5*emodule/(1.0+poisson); % shear modulus
shcof=5/6; % shear correction factor
matmtss=shearm*shcof*t*[1 0; 0 1]; % shear material property
for iel=1:nel % loop for the total number of elements
for i=1:nnel
nd(i)=nodes(iel,i); % extract connected node for (iel)-th element
xcoord(i)=gcoord(nd(i),1); % extract x value of the node
ycoord(i)=gcoord(nd(i),2); % extract y value of the node
zcoord(i)=gcoord(nd(i),3); % extract z value of the node
end
%
% compute the local direction cosines and local axes
%
[tr3d,xprime,yprime]=fetransh(xcoord,ycoord,zcoord,nnel);
%
k=zeros(edof,edof); % element matrix in local axes
ke=zeros(edof,edof); % element matrix in global axes
km=zeros(edof,edof); % element membrane matrix
kb=zeros(edof,edof); % element bending matrix
ks=zeros(edof,edof); % element shear matrix
%------------------------------------------------------
% numerical integration for bending term
%------------------------------------------------------
for intx=1:nglxb
x=pointb(intx,1); % sampling point in x-axis
wtx=weightb(intx,1); % weight in x-axis
for inty=1:nglyb
y=pointb(inty,2); % sampling point in y-axis
wty=weightb(inty,2) ; % weight in y-axis
[shape,dhdr,dhds]=feisoq4(x,y); % compute shape functions and
% derivatives at sampling point
jacob2=fejacob2(nnel,dhdr,dhds,xprime,yprime); % compute Jacobian
detjacob=det(jacob2); % determinant of Jacobian
invjacob=inv(jacob2); % inverse of Jacobian matrix
[dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob); % derivatives w.r.t.
% physical coordinate
kinmtsb=fekinesb(nnel,dhdx,dhdy); % bending kinematic matrix
kinmtsm=fekinesm(nnel,dhdx,dhdy); % membrane kinematic matrix
%--------------------------------------------
% compute bending element matrix
%--------------------------------------------
kb=kb+kinmtsb'*matmtsb*kinmtsb*wtx*wty*detjacob;
km=km+kinmtsm'*matmtsm*kinmtsm*wtx*wty*detjacob;
end
end % end of numerical integration loop for bending term
%------------------------------------------------------
% numerical integration for bending term
%------------------------------------------------------
for intx=1:nglxs
x=points(intx,1); % sampling point in x-axis
wtx=weights(intx,1); % weight in x-axis
for inty=1:nglys
y=points(inty,2); % sampling point in y-axis
wty=weights(inty,2) ; % weight in y-axis
[shape,dhdr,dhds]=feisoq4(x,y); % compute shape functions and
% derivatives at sampling point
jacob2=fejacob2(nnel,dhdr,dhds,xprime,yprime); % compute Jacobian
detjacob=det(jacob2); % determinant of Jacobian
invjacob=inv(jacob2); % inverse of Jacobian matrix
[dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob); % derivatives w.r.t.
% physical coordinate
kinmtss=fekiness(nnel,dhdx,dhdy,shape); % shear kinematic matrix
%----------------------------------------
% compute shear element matrix
%----------------------------------------
ks=ks+kinmtss'*matmtss*kinmtss*wtx*wty*detjacob;
end
end % end of numerical integration loop for shear term
%--------------------------------
% compute element matrix
%--------------------------------
k=km+kb+ks;
%-----------------------------------------------
% transform from local to global systems
%-----------------------------------------------
ke=tr3d'*k*tr3d;
index=feeldof(nd,nnel,ndof);% extract system dofs associated with element
kk=feasmbl1(kk,ke,index); % assemble element matrices
end
%-------------------------------------
% check the singular drilling dof
%-------------------------------------
for i=1:sdof
if(abs(kk(i,i)) < 1e-5)
sum=0.0;
for j=1:sdof
sum=sum+abs(kk(i,j));
end
if (sum < 1e-5)
kk(i,i)=1;
end
end
end
%-----------------------------
% apply boundary conditions
%-----------------------------
[kk,ff]=feaplyc2(kk,ff,bcdof,bcval);
%----------------------------
% solve the matrix equation
%----------------------------
disp=kk\ff;
num=1:1:sdof;
displace=[num' disp] % print nodal displacements
%--------------------------------------------------------------
% Sub_routines
%%%%%%%%%%%%%%%%%
function [point2,weight2]=feglqd2(nglx,ngly)
%-------------------------------------------------------------------
% Purpose:
% determine the integration points and weighting coefficients
% of Gauss-Legendre quadrature for two-dimensional integration
%
% Synopsis:
% [point2,weight2]=feglqd2(nglx,ngly)
%
% Variable Description:
% nglx - number of integration points in the x-axis
% ngly - number of integration points in the y-axis
% point2 - vector containing integration points
% weight2 - vector containing weighting coefficients
%-------------------------------------------------------------------
% determine the largest one between nglx and ngly
if nglx > ngly
ngl=nglx;
else
ngl=ngly;
end
% initialization
point2=zeros(ngl,2);
weight2=zeros(ngl,2);
% find corresponding integration points and weights
[pointx,weightx]=feglqd1(nglx); % quadrature rule for x-axis
[pointy,weighty]=feglqd1(ngly); % quadrature rule for y-axis
% quadrature for two-dimension
for intx=1:nglx % quadrature in x-axis
point2(intx,1)=pointx(intx);
weight2(intx,1)=weightx(intx);
end
for inty=1:ngly % quadrature in y-axis
point2(inty,2)=pointy(inty);
weight2(inty,2)=weighty(inty);
end
_______________________________________________________________________
function [matmtrx]=fematiso(iopt,elastic,poisson)
%------------------------------------------------------------------------
% Purpose:
% determine the constitutive equation for isotropic material
%
% Synopsis:
% [matmtrx]=fematiso(iopt,elastic,poisson)
%
% Variable Description:
% elastic - elastic modulus
% poisson - Poisson's ratio
% iopt=1 - plane stress analysis
% iopt=2 - plane strain analysis
% iopt=3 - axisymmetric analysis
% iopt=4 - three dimensional analysis
%------------------------------------------------------------------------
if iopt==1 % plane stress
matmtrx= elastic/(1-poisson*poisson)* ...
[1 poisson 0; ...
poisson 1 0; ...
0 0 (1-poisson)/2];
elseif iopt==2 % plane strain
matmtrx= elastic/((1+poisson)*(1-2*poisson))* ...
[(1-poisson) poisson 0;
poisson (1-poisson) 0;
0 0 (1-2*poisson)/2];
elseif iopt==3 % axisymmetry
matmtrx= elastic/((1+poisson)*(1-2*poisson))* ...
[(1-poisson) poisson poisson 0;
poisson (1-poisson) poisson 0;
poisson poisson (1-poisson) 0;
0 0 0 (1-2*poisson)/2];
else % three-dimension
matmtrx= elastic/((1+poisson)*(1-2*poisson))* ...
[(1-poisson) poisson poisson 0 0 0;
poisson (1-poisson) poisson 0 0 0;
poisson poisson (1-poisson) 0 0 0;
0 0 0 (1-2*poisson)/2 0 0;
0 0 0 0 (1-2*poisson)/2 0;
0 0 0 0 0 (1-2*poisson)/2];
end
______________________________________________________________________
function [tr3d,xprime,yprime]=fetransh(xcoord,ycoord,zcoord,n)
%--------------------------------------------------------------
% Purpose:
% Compute direction cosines between three-dimensional
% local and global coordinate axes
%
% Synopsis:
% [tr3d,xprime,yprime]=fetransh(xcoord,ycoord,zcoord,n)
%
% Variable Description:
% xcoord - nodal x coordinates (4x1)
% ycoord - nodal y coordinates (4X1)
% zcoord - nodal z coordinates (4X1)
% n - number of nodes per element
% tr3d - 3d transformation matrix from local to global axes
% xprime - coordinate in terms of the local axes (4x1)
% yprime - coordinate in terms of the global axes (4X1)
%
% Note:
% The local x-axis is defined in the direction from the first node
% to the second node. Nodes 1, 2 and 4 define the local xy-plane.
% The local z-axis is defined normal to the local xy-plane.
% The local y-axis is defined normal to the x and z axes.
%--------------------------------------------------------------------------
%
% compute direction cosines
%
v12x=xcoord(2)-xcoord(1);
v12y=ycoord(2)-ycoord(1);
v12z=zcoord(2)-zcoord(1);
l12=sqrt(v12x^2+v12y^2+v12z^2);
v23x=xcoord(3)-xcoord(2);
v23y=ycoord(3)-ycoord(2);
v23z=zcoord(3)-zcoord(2);
l23=sqrt(v23x^2+v23y^2+v23z^2);
v34x=xcoord(4)-xcoord(3);
v34y=ycoord(4)-ycoord(3);
v34z=zcoord(4)-zcoord(3);
l34=sqrt(v34x^2+v34y^2+v34z^2);
v14x=xcoord(4)-xcoord(1);
v14y=ycoord(4)-ycoord(1);
v14z=zcoord(4)-zcoord(1);
l14=sqrt(v14x^2+v14y^2+v14z^2);
v13x=xcoord(3)-xcoord(1);
v13y=ycoord(3)-ycoord(1);
v13z=zcoord(3)-zcoord(1);
l13=sqrt(v13x^2+v13y^2+v13z^2);
v1tx=v12y*v14z-v12z*v14y;
v1ty=v12z*v14x-v12x*v14z;
v1tz=v12x*v14y-v12y*v14x;
v1yx=v1ty*v12z-v1tz*v12y;
v1yy=v1tz*v12x-v1tx*v12z;
v1yz=v1tx*v12y-v1ty*v12x;
vxx=v12x/l12;
vxy=v12y/l12;
vxz=v12z/l12;
vyx=v1yx/sqrt(v1yx^2+v1yy^2+v1yz^2);
vyy=v1yy/sqrt(v1yx^2+v1yy^2+v1yz^2);
vyz=v1yz/sqrt(v1yx^2+v1yy^2+v1yz^2);
vzx=v1tx/sqrt(v1tx^2+v1ty^2+v1tz^2);
vzy=v1ty/sqrt(v1tx^2+v1ty^2+v1tz^2);
vzz=v1tz/sqrt(v1tx^2+v1ty^2+v1tz^2);
%
% transformation matrix
%
for i=1:2*n
i1=(i-1)*3+1;
i2=i1+1;
i3=i2+1;
tr3d(i1,i1)=vxx;
tr3d(i1,i2)=vxy;
tr3d(i1,i3)=vxz;
tr3d(i2,i1)=vyx;
tr3d(i2,i2)=vyy;
tr3d(i2,i3)=vyz;
tr3d(i3,i1)=vzx;
tr3d(i3,i2)=vzy;
tr3d(i3,i3)=vzz;
end
%
% compute nodal values in terms of local axes
%
alpa213=acos((l12^2+l13^2-l23^2)/(2*l12*l13));
alpa314=acos((l13^2+l14^2-l34^2)/(2*l13*l14));
alpa41y=2*atan(1)-alpa213-alpa314;
xprime(1)=0; yprime(1)=0;
xprime(2)=l12; yprime(2)=0;
xprime(3)=l13*cos(alpa213); yprime(3)=l13*sin(alpa213);
xprime(4)=l14*sin(alpa41y); yprime(4)=l14*cos(alpa41y);
_____________________________________________________________________
function [shapeq4,dhdrq4,dhdsq4]=feisoq4(rvalue,svalue)
%------------------------------------------------------------------------
% Purpose:
% compute isoparametric four-node quadilateral shape functions
% and their derivatves at the selected (integration) point
% in terms of the natural coordinate
%
% Synopsis:
% [shapeq4,dhdrq4,dhdsq4]=feisoq4(rvalue,svalue)
%
% Variable Description:
% shapeq4 - shape functions for four-node element
% dhdrq4 - derivatives of the shape functions w.r.t. r
% dhdsq4 - derivatives of the shape functions w.r.t. s
% rvalue - r coordinate value of the selected point
% svalue - s coordinate value of the selected point
%
% Notes:
% 1st node at (-1,-1), 2nd node at (1,-1)
% 3rd node at (1,1), 4th node at (-1,1)
%------------------------------------------------------------------------
% shape functions
shapeq4(1)=0.25*(1-rvalue)*(1-svalue);
shapeq4(2)=0.25*(1+rvalue)*(1-svalue);
shapeq4(3)=0.25*(1+rvalue)*(1+svalue);
shapeq4(4)=0.25*(1-rvalue)*(1+svalue);
% derivatives
dhdrq4(1)=-0.25*(1-svalue);
dhdrq4(2)=0.25*(1-svalue);
dhdrq4(3)=0.25*(1+svalue);
dhdrq4(4)=-0.25*(1+svalue);
dhdsq4(1)=-0.25*(1-rvalue);
dhdsq4(2)=-0.25*(1+rvalue);
dhdsq4(3)=0.25*(1+rvalue);
dhdsq4(4)=0.25*(1-rvalue);
________________________________________________________________________
function [jacob2]=fejacob2(nnel,dhdr,dhds,xcoord,ycoord)
%------------------------------------------------------------------------
% Purpose:
% determine the Jacobian for two-dimensional mapping
%
% Synopsis:
% [jacob2]=fejacob2(nnel,dhdr,dhds,xcoord,ycoord)
%
% Variable Description:
% jacob2 - Jacobian for one-dimension
% nnel - number of nodes per element
% dhdr - derivative of shape functions w.r.t. natural coordinate r
% dhds - derivative of shape functions w.r.t. natural coordinate s
% xcoord - x axis coordinate values of nodes
% ycoord - y axis coordinate values of nodes
%------------------------------------------------------------------------
jacob2=zeros(2,2);
for i=1:nnel
jacob2(1,1)=jacob2(1,1)+dhdr(i)*xcoord(i);
jacob2(1,2)=jacob2(1,2)+dhdr(i)*ycoord(i);
jacob2(2,1)=jacob2(2,1)+dhds(i)*xcoord(i);
jacob2(2,2)=jacob2(2,2)+dhds(i)*ycoord(i);
end
________________________________________________________________________
function [dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob)
%------------------------------------------------------------------------
% Purpose:
% determine derivatives of 2-D isoparametric shape functions with
% respect to physical coordinate system
%
% Synopsis:
% [dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob)
%
% Variable Description:
% dhdx - derivative of shape function w.r.t. physical coordinate x
% dhdy - derivative of shape function w.r.t. physical coordinate y
% nnel - number of nodes per element
% dhdr - derivative of shape functions w.r.t. natural coordinate r
% dhds - derivative of shape functions w.r.t. natural coordinate s
% invjacob - inverse of 2-D Jacobian matrix
%------------------------------------------------------------------------
for i=1:nnel
dhdx(i)=invjacob(1,1)*dhdr(i)+invjacob(1,2)*dhds(i);
dhdy(i)=invjacob(2,1)*dhdr(i)+invjacob(2,2)*dhds(i);
end
_________________________________________________________________________
function [kinmtsb]=fekinesb(nnel,dhdx,dhdy)
%--------------------------------------------------------------------------
% Purpose:
% determine the kinematic matrix expression relating bending curvatures
% to rotations and displacements for shear deformable plate bending
%
% Synopsis:
% [kinmtsb]=fekinesb(nnel,dhdx,dhdy)
%
% Variable Description:
% nnel - number of nodes per element
% dhdx - derivatives of shape functions with respect to x
% dhdy - derivatives of shape functions with respect to y
%--------------------------------------------------------------------------
for i=1:nnel
i1=(i-1)*6+1;
i2=i1+1;
i3=i2+1;
i4=i3+1;
i5=i4+1;
i6=i5+1;
kinmtsb(1,i5)=dhdx(i);
kinmtsb(2,i4)=-dhdy(i);
kinmtsb(3,i5)=dhdy(i);
kinmtsb(3,i4)=-dhdx(i);
kinmtsb(3,i6)=0;
end
_________________________________________________________________________
function [kinmtsm]=fekinesm(nnel,dhdx,dhdy)
%------------------------------------------------------------------------
% Purpose:
% determine the kinematic equation between strains and displacements
% for two-dimensional solids
%
% Synopsis:
% [kinmtsm]=fekinesm(nnel,dhdx,dhdy)
%
% Variable Description:
% nnel - number of nodes per element
% dhdx - derivatives of shape functions with respect to x
% dhdy - derivatives of shape functions with respect to y
%------------------------------------------------------------------------
for i=1:nnel
i1=(i-1)*6+1;
i2=i1+1;
i3=i2+1;
i4=i3+1;
i5=i4+1;
i6=i5+1;
kinmtsm(1,i1)=dhdx(i);
kinmtsm(2,i2)=dhdy(i);
kinmtsm(3,i1)=dhdy(i);
kinmtsm(3,i2)=dhdx(i);
kinmtsm(3,i6)=0.0;
end
__________________________________________________________________________
function [kinmtss]=fekiness(nnel,dhdx,dhdy,shape)
%------------------------------------------------------------------------
% Purpose:
% determine the kinematic matrix expression relating shear strains
% to rotations and displacements for shear deformable plate bending
%
% Synopsis:
% [kinmtss]=fekiness(nnel,dhdx,dhdy,shape)
%
% Variable Description:
% nnel - number of nodes per element
% dhdx - derivatives of shape functions with respect to x
% dhdy - derivatives of shape functions with respect to y
% shape - shape function
%------------------------------------------------------------------------
for i=1:nnel
i1=(i-1)*6+1;
i2=i1+1;
i3=i2+1;
i4=i3+1;
i5=i4+1;
i6=i5+1;
kinmtss(1,i3)=dhdx(i);
kinmtss(1,i5)=shape(i);
kinmtss(2,i3)=dhdy(i);
kinmtss(2,i4)=-shape(i);
kinmtss(2,i6)=0.0;
end
_________________________________________________________________________
function [index]=feeldof(nd,nnel,ndof)
%----------------------------------------------------------
% Purpose:
% Compute system dofs associated with each element
%
% Synopsis:
% [index]=feeldof(nd,nnel,ndof)
%
% Variable Description:
% index - system dof vector associated with element "iel"
% iel - element number whose system dofs are to be determined
% nnel - number of nodes per element
% ndof - number of dofs per node
%-----------------------------------------------------------
edof = nnel*ndof;
k=0;
for i=1:nnel
start = (nd(i)-1)*ndof;
for j=1:ndof
k=k+1;
index(k)=start+j;
end
end
________________________________________________________________________
function [kk]=feasmbl1(kk,k,index)
%----------------------------------------------------------
% Purpose:
% Assembly of element matrices into the system matrix
%
% Synopsis:
% [kk]=feasmbl1(kk,k,index)
%
% Variable Description:
% kk - system matrix
% k - element matri
% index - d.o.f. vector associated with an element
%-----------------------------------------------------------
edof = length(index);
for i=1:edof
ii=index(i);
for j=1:edof
jj=index(j);
kk(ii,jj)=kk(ii,jj)+k(i,j);
end
end
________________________________________________________________________
function [kk,ff]=feaplyc2(kk,ff,bcdof,bcval)
%----------------------------------------------------------
% Purpose:
% Apply constraints to matrix equation [kk]{x}={ff}
%
% Synopsis:
% [kk,ff]=feaplybc(kk,ff,bcdof,bcval)
%
% Variable Description:
% kk - system matrix before applying constraints
% ff - system vector before applying constraints
% bcdof - a vector containging constrained d.o.f
% bcval - a vector containing contained value
%
% For example, there are constraints at d.o.f=2 and 10
% and their constrained values are 0.0 and 2.5,
% respectively. Then, bcdof(1)=2 and bcdof(2)=10; and
% bcval(1)=1.0 and bcval(2)=2.5.
%-----------------------------------------------------------
n=length(bcdof);
sdof=size(kk);
for i=1:n
c=bcdof(i);
for j=1:sdof
kk(c,j)=0;
end
kk(c,c)=1;
ff(c)=bcval(i);
end
________________________________________________________________________

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2013-3-21
The error was made on line 71 by Colonel Mustard with the Wrench.
No? How about a Clue ?
  1 个评论
Deviprakash Upadhyaya
Deviprakash Upadhyaya 2013-3-22
Sorry for the problem created.I have edited the question and made it proper.Will be getting a warning message in " disp=kk\ff; ".Problem with inversing the system matrix " kk ".

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 General Applications 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by