Advanced Mathematics and Mechanics Applications Using MATLAB phần 4 doc

184 ADVANCED MATH AND MECHANICS APPLICATIONS USING MATLAB

−2

−1

0

1

2

0

1

2

3

4

0

0.5

1

1.5

2

2.5

3

3.5

4

x axis

Surface Plot of a General Polyhedron

y axis

z axis

Figure 5.6: Surface Plot of a General Polyhedron

© 2003 by CRC Press LLC

GAUSS INTEGRATION WITH GEOMETRIC PROPERTY APPLICATIONS 185

Program polhdrun

1: function polhdrun

2: % Example: polhdrun

3: % ~~~~~~~~~~~~~~~~~

4: %

5: % This program illustrates the use of routine

6: % polhedrn to calculate the geometrical

7: % properties of a polyhedron.

8: %

9: % User m functions called: 10: % crosmat, polyxy, cubrange, pyramid,

11: % polhdplt, polhedrn

12:

13: x=[2 2 2 2 2 2 0 0 0 0 0 0]-1;

14: y=[0 4 4 2 3 3 0 4 4 2 3 3];

15: z=[0 0 4 1 1 2 0 0 4 1 1 2]; 16: idface=[1 2 3 6 5 4 6 3; ...

17: 1 3 9 7 0 0 0 0; ... 18: 1 7 8 2 0 0 0 0; ...

19: 2 8 9 3 0 0 0 0; ...

20: 7 9 12 10 11 12 9 8; ... 21: 4 10 12 6 0 0 0 0; ...

22: 4 5 11 10 0 0 0 0; ...

23: 5 6 12 11 0 0 0 0]; 24: polhdplt(x,y,z,idface,[1,1,1]);

25: [v,rc,vrr,irr]=polhedrn(x,y,z,idface) 26:

27: %=============================================

28:

29: function [v,rc,vrr,irr]=polhedrn(x,y,z,idface)

30: %

31: % [v,rc,vrr,irr]=polhedrn(x,y,z,idface) 32: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

33: % 34: % This function determines the volume,

35: % centroidal coordinates and inertial moments

36: % for an arbitrary polyhedron. 37: %

38: % x,y,z - vectors containing the corner

39: % indices of the polyhedron 40: % idface - a matrix in which row j defines the

41: % corner indices of the j’th face.

© 2003 by CRC Press LLC

186 ADVANCED MATH AND MECHANICS APPLICATIONS USING MATLAB

42: % Each face is traversed in a

43: % counterclockwise sense relative to

44: % the outward normal. The column 45: % dimension equals the largest number

46: % of indices needed to define a face. 47: % Rows requiring fewer than the

48: % maximum number of corner indices are

49: % padded with zeros on the right. 50: %

51: % v - the volume of the polyhedron

52: % rc - the centroidal radius 53: % vrr - the integral of R*R’*d(vol)

54: % irr - the inertia tensor for a rigid body

55: % of unit mass obtained from vrr as 56: % eye(3,3)*sum(diag(vrr))-vrr

57: % 58: % User m functions called: pyramid

59: %----------------------------------------------

60:

61: r=[x(:),y(:),z(:)]; nf=size(idface,1);

62: v=0; vr=0; vrr=0;

63: for k=1:nf 64: i=idface(k,:); i=i(find(i>0));

65: [u,ur,urr]=pyramid(r(i,:)); 66: v=v+u; vr=vr+ur; vrr=vrr+urr;

67: end

68: rc=vr/v; irr=eye(3,3)*sum(diag(vrr))-vrr; 69:

70: %=============================================

71:

72: function [area,xbar,ybar,axx,axy,ayy]=polyxy(x,y)

73: % 74: % [area,xbar,ybar,axx,axy,ayy]=polyxy(x,y)

75: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

76: % 77: % This function computes the area, centroidal

78: % coordinates, and inertial moments of an

79: % arbitrary polygon. 80: %

81: % x,y - vectors containing the corner 82: % coordinates. The boundary is

83: % traversed in a counterclockwise

84: % direction 85: %

86: % area - the polygon area

© 2003 by CRC Press LLC

GAUSS INTEGRATION WITH GEOMETRIC PROPERTY APPLICATIONS 187

87: % xbar,ybar - the centroidal coordinates

88: % axx - integral of x^2*dxdy

89: % axy - integral of xy*dxdy 90: % ayy - integral of y^2*dxdy

91: % 92: % User m functions called: none

93: %----------------------------------------------

94:

95: n=1:length(x); n1=n+1;

96: x=[x(:);x(1)]; y=[y(:);y(1)];

97: a=(x(n).*y(n1)-y(n).*x(n1))’; 98: area=sum(a)/2; a6=6*area;

99: xbar=a*(x(n)+x(n1))/a6; ybar=a*(y(n)+y(n1))/a6;

100: ayy=a*(y(n).^2+y(n).*y(n1)+y(n1).^2)/12; 101: axy=a*(x(n).*(2*y(n)+y(n1))+x(n1).* ...

102: (2*y(n1)+y(n)))/24; 103: axx=a*(x(n).^2+x(n).*x(n1)+x(n1).^2)/12;

104:

105: %============================================= 106:

107: function [v,vr,vrr,h,area,n]=pyramid(r)

108: % 109: % [v,vr,vrr,h,area,n]=pyramid(r)

110: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 111: %

112: % This function determines geometrical

113: % properties of a pyramid with the apex at the 114: % origin and corner coordinates of the base

115: % stored in the rows of r.

116: % 117: % r - matrix containing the corner

118: % coordinates of a polygonal base stored 119: % in the rows of matrix r.

120: %

121: % v - the volume of the pyramid 122: % vr - the first moment of volume relative to

123: % the origin

124: % vrr - the second moment of volume relative 125: % to the origin

126: % h - the pyramid height 127: % area - the base area

128: % n - the outward directed unit normal to

129: % the base 130: %

131: % User m functions called: crosmat, polyxy

© 2003 by CRC Press LLC

188 ADVANCED MATH AND MECHANICS APPLICATIONS USING MATLAB

132: %----------------------------------------------

133:

134: ns=size(r,1); 135: na=sum(crosmat(r,r([2:ns,1],:)))’/2;

136: area=norm(na); n=na/area; p=null(n’); 137: i=p(:,1); j=p(:,2);

138: if det([p,n])<0, j=-j; end;

139: r1=r(1,:); rr=r-r1(ones(ns,1),:); 140: x=rr*i; y=rr*j;

141: [areat,xc,yc,axx,axy,ayy]=polyxy(x,y);

142: rc=r1’+xc*i+yc*j; h=r1*n; 143: v=h*area/3; vr=v*3/4*rc;

144: axx=axx-area*xc^2; ayy=ayy-area*yc^2;

145: axy=axy-area*xc*yc; 146: vrr=h/5*(area*rc*rc’+axx*i*i’+ayy*j*j’+ ...

147: axy*(i*j’+j*i’)); 148:

149: %=============================================

150:

151: function polhdplt(x,y,z,idface,colr)

152: %

153: % polhdplt(x,y,z,idface,colr) 154: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~

155: % 156: % This function makes a surface plot of an

157: % arbitrary polyhedron.

158: % 159: % x,y,z - vectors containing the corner

160: % indices of the polyhedron

161: % idface - a matrix in which row j defines the 162: % corner indices of the j’th face.

163: % Each face is traversed in a 164: % counterclockwise sense relative to

165: % the outward normal. The column

166: % dimension equals the largest number 167: % of indices needed to define a face.

168: % Rows requiring fewer than the

169: % maximum number of corner indices are 170: % padded with zeros on the right.

171: % colr - character string or a vector 172: % defining the surface color

173: %

174: % User m functions called: cubrange 175: %----------------------------------------------

176:

© 2003 by CRC Press LLC

GAUSS INTEGRATION WITH GEOMETRIC PROPERTY APPLICATIONS 189

177: if nargin<5, colr=[1 0 1]; end

178: hold off, close; nf=size(idface,1);

179: v=cubrange([x(:),y(:),z(:)],1.1); 180: for k=1:nf

181: i=idface(k,:); i=i(find(i>0)); 182: xi=x(i); yi=y(i); zi=z(i);

183: fill3(xi,yi,zi,colr); hold on;

184: end 185: axis(v); grid on;

186: xlabel(’x axis’); ylabel(’y axis’);

187: zlabel(’z axis’); 188: title(’Surface Plot of a General Polyhedron’);

189: figure(gcf); hold off;

190:

191: %=============================================

192:

193: function c=crosmat(a,b)

194: %

195: % c=crosmat(a,b) 196: % ~~~~~~~~~~~~~~

197: %

198: % This function computes the vector cross 199: % product for vectors stored in the rows

200: % of matrices a and b, and returns the 201: % results in the rows of c.

202: %

203: % User m functions called: none 204: %----------------------------------------------

205:

206: c=[a(:,2).*b(:,3)-a(:,3).*b(:,2),... 207: a(:,3).*b(:,1)-a(:,1).*b(:,3),...

208: a(:,1).*b(:,2)-a(:,2).*b(:,1)]; 209:

210: %=============================================

211:

212: % function range=cubrange(xyz,ovrsiz)

213: % See Appendix B

© 2003 by CRC Press LLC