CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
shape8hr.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine shape8hr (xl, xsj, shp, gs, a)
 

Function/Subroutine Documentation

◆ shape8hr()

subroutine shape8hr ( real*8, dimension(3,20), intent(in)  xl,
real*8, intent(out)  xsj,
real*8, dimension(4,20), intent(out)  shp,
real*8, dimension(8,4), intent(out)  gs,
real*8, intent(out)  a 
)
20 !
21 ! shape functions and derivatives for a 8-node linear isoparametric
22 ! mean strain solid element with hourglass control
23 !
24 ! Reference: Flanagan, D.P., Belytschko, T.; "Uniform strain hexahedron
25 ! and quadrilateral with orthogonal Hourglass control". Int. J. Num.
26 ! Meth. Engg., Vol. 17, 679-706, 1981.
27 !
28 ! author: Otto-Ernst Bernhardi
29 !
30  implicit none
31 !
32  integer i,j,k,l
33  real*8 shp(4,20),xl(3,20),xsj,vol
34  real*8 x1,x2,x3,x4,x5,x6,x7,x8
35  real*8 y1,y2,y3,y4,y5,y6,y7,y8
36  real*8 z1,z2,z3,z4,z5,z6,z7,z8
37  real*8 gb(8,4),gs(8,4),s0,a
38 !
39  intent(in) xl
40 !
41  intent(out) shp,gs,a,xsj
42 !
43  gb = reshape((
44  & / 1.0d0, 1.0d0,-1.0d0,-1.0d0,-1.0d0,-1.0d0, 1.0d0, 1.0d0,
45  & 1.0d0,-1.0d0,-1.0d0, 1.0d0,-1.0d0, 1.0d0, 1.0d0,-1.0d0,
46  & 1.0d0,-1.0d0, 1.0d0,-1.0d0, 1.0d0,-1.0d0, 1.0d0,-1.0d0,
47  & -1.0d0, 1.0d0,-1.0d0, 1.0d0, 1.0d0,-1.0d0, 1.0d0,-1.0d0 /),
48  & (/8,4/))
49 !
50 ! shape functions and their global derivatives
51 !
52 ! shape functions
53 !
54  shp(4, 1)=1.0d0/8.0d0
55  shp(4, 2)=1.0d0/8.0d0
56  shp(4, 3)=1.0d0/8.0d0
57  shp(4, 4)=1.0d0/8.0d0
58  shp(4, 5)=1.0d0/8.0d0
59  shp(4, 6)=1.0d0/8.0d0
60  shp(4, 7)=1.0d0/8.0d0
61  shp(4, 8)=1.0d0/8.0d0
62 !
63 ! extract node point coordinates of element
64  x1=xl(1,1)
65  x2=xl(1,2)
66  x3=xl(1,3)
67  x4=xl(1,4)
68  x5=xl(1,5)
69  x6=xl(1,6)
70  x7=xl(1,7)
71  x8=xl(1,8)
72  y1=xl(2,1)
73  y2=xl(2,2)
74  y3=xl(2,3)
75  y4=xl(2,4)
76  y5=xl(2,5)
77  y6=xl(2,6)
78  y7=xl(2,7)
79  y8=xl(2,8)
80  z1=xl(3,1)
81  z2=xl(3,2)
82  z3=xl(3,3)
83  z4=xl(3,4)
84  z5=xl(3,5)
85  z6=xl(3,6)
86  z7=xl(3,7)
87  z8=xl(3,8)
88 !
89 ! Average displacement derivative operator matrix,
90 ! using eqn 16 in the reference above.
91 ! Generated using maxima/wxmaxima and the following input lines.
92 ! Note that shp array must be divided by the element volume.
93 ! h1:(1-r)*(1-s)*(1-t)/8;
94 ! h2:(1+r)*(1-s)*(1-t)/8;
95 ! h3:(1+r)*(1+s)*(1-t)/8;
96 ! h4:(1-r)*(1+s)*(1-t)/8;
97 ! h5:(1-r)*(1-s)*(1+t)/8;
98 ! h6:(1+r)*(1-s)*(1+t)/8;
99 ! h7:(1+r)*(1+s)*(1+t)/8;
100 ! h8:(1-r)*(1+s)*(1+t)/8;
101 ! H:matrix([h1,h2,h3,h4,h5,h6,h7,h8])$
102 ! Bx:diff(H,r);
103 ! By:diff(H,s);
104 ! Bz:diff(H,t);
105 ! B:matrix([Bx[1,1],Bx[1,2],Bx[1,3],Bx[1,4],Bx[1,5],Bx[1,6],Bx[1,7],Bx[1,8]],
106 ! [By[1,1],By[1,2],By[1,3],By[1,4],By[1,5],By[1,6],By[1,7],By[1,8]],
107 ! [Bz[1,1],Bz[1,2],Bz[1,3],Bz[1,4],Bz[1,5],Bz[1,6],Bz[1,7],Bz[1,8]]);
108 ! x:matrix([x1,x2,x3,x4,x5,x6,x7,x8],
109 ! [y1,y2,y3,y4,y5,y6,y7,y8],
110 ! [z1,z2,z3,z4,z5,z6,z7,z8]);
111 ! xr:B.transpose(x);
112 ! det:determinant(xr);
113 ! vol:factor(integrate(integrate(integrate(det,r,-1,1),s,-1,1),t,-1,1));
114 ! shp: matrix( [0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0])$
115 ! shp[1][1]:ratsimp(diff(vol,x1));
116 ! shp[1][2]:ratsimp(diff(vol,x2));
117 ! shp[1][3]:ratsimp(diff(vol,x3));
118 ! ... (more lines)
119 ! shp[3][8]:ratsimp(diff(vol,z8));
120 ! fortran(shp);
121 !
122 ! local derivatives of the shape functions: xi-derivative
123 !
124  shp(1,1) = ((y5-y4)*z8+(y2-y5)*z6+(-y8+y6-y4+y2)*z5+(y8+y5-y3-y2)
125  1 *z4+(y4-y2)*z3+(-y6-y5+y4+y3)*z2)/1.2d+1
126  shp(1,2) = -((y6-y3)*z7+(-y7+y5-y3+y1)*z6+(y1-y6)*z5+(y3-y1)*z4+(
127  1 y7+y6-y4-y1)*z3+(-y6-y5+y4+y3)*z1)/1.2d+1
128  shp(1,3) = -((y7-y4)*z8+(-y8+y6-y4+y2)*z7+(y2-y7)*z6+(y8+y7-y2-y1
129  1 )*z4+(-y7-y6+y4+y1)*z2+(y4-y2)*z1)/1.2d+1
130  shp(1,4) = -((y7-y5+y3-y1)*z8+(y3-y8)*z7+(y8-y1)*z5+(-y8-y7+y2+y1
131  1 )*z3+(y1-y3)*z2+(y8+y5-y3-y2)*z1)/1.2d+1
132  shp(1,5) = ((y7+y6-y4-y1)*z8+(y6-y8)*z7+(-y8-y7+y2+y1)*z6+(y8-y1)
133  1 *z4+(y1-y6)*z2+(y8-y6+y4-y2)*z1)/1.2d+1
134  shp(1,6) = ((y7-y5)*z8+(-y8-y5+y3+y2)*z7+(y8+y7-y2-y1)*z5+(y2-y7)
135  1 *z3+(-y7+y5-y3+y1)*z2+(y5-y2)*z1)/1.2d+1
136  shp(1,7) = -((y6+y5-y4-y3)*z8+(-y8-y5+y3+y2)*z6+(y6-y8)*z5+(y8-y3
137  1 )*z4+(y8-y6+y4-y2)*z3+(y3-y6)*z2)/1.2d+1
138  shp(1,8) = ((y6+y5-y4-y3)*z7+(y5-y7)*z6+(-y7-y6+y4+y1)*z5+(y7-y5+
139  1 y3-y1)*z4+(y7-y4)*z3+(y4-y5)*z1)/1.2d+1
140 !
141 ! local derivatives of the shape functions: eta-derivative
142 !
143  shp(2,1) = -((x5-x4)*z8+(x2-x5)*z6+(-x8+x6-x4+x2)*z5+(x8+x5-x3-x2
144  1 )*z4+(x4-x2)*z3+(-x6-x5+x4+x3)*z2)/1.2d+1
145  shp(2,2) = ((x6-x3)*z7+(-x7+x5-x3+x1)*z6+(x1-x6)*z5+(x3-x1)*z4+(x
146  1 7+x6-x4-x1)*z3+(-x6-x5+x4+x3)*z1)/1.2d+1
147  shp(2,3) = ((x7-x4)*z8+(-x8+x6-x4+x2)*z7+(x2-x7)*z6+(x8+x7-x2-x1)
148  1 *z4+(-x7-x6+x4+x1)*z2+(x4-x2)*z1)/1.2d+1
149  shp(2,4) = ((x7-x5+x3-x1)*z8+(x3-x8)*z7+(x8-x1)*z5+(-x8-x7+x2+x1)
150  1 *z3+(x1-x3)*z2+(x8+x5-x3-x2)*z1)/1.2d+1
151  shp(2,5) = -((x7+x6-x4-x1)*z8+(x6-x8)*z7+(-x8-x7+x2+x1)*z6+(x8-x1
152  1 )*z4+(x1-x6)*z2+(x8-x6+x4-x2)*z1)/1.2d+1
153  shp(2,6) = -((x7-x5)*z8+(-x8-x5+x3+x2)*z7+(x8+x7-x2-x1)*z5+(x2-x7
154  1 )*z3+(-x7+x5-x3+x1)*z2+(x5-x2)*z1)/1.2d+1
155  shp(2,7) = ((x6+x5-x4-x3)*z8+(-x8-x5+x3+x2)*z6+(x6-x8)*z5+(x8-x3)
156  1 *z4+(x8-x6+x4-x2)*z3+(x3-x6)*z2)/1.2d+1
157  shp(2,8) = -((x6+x5-x4-x3)*z7+(x5-x7)*z6+(-x7-x6+x4+x1)*z5+(x7-x5
158  1 +x3-x1)*z4+(x7-x4)*z3+(x4-x5)*z1)/1.2d+1
159 !
160 ! local derivatives of the shape functions: zeta-derivative
161 !
162  shp(3,1) = ((x5-x4)*y8+(x2-x5)*y6+(-x8+x6-x4+x2)*y5+(x8+x5-x3-x2)
163  1 *y4+(x4-x2)*y3+(-x6-x5+x4+x3)*y2)/1.2d+1
164  shp(3,2) = -((x6-x3)*y7+(-x7+x5-x3+x1)*y6+(x1-x6)*y5+(x3-x1)*y4+(
165  1 x7+x6-x4-x1)*y3+(-x6-x5+x4+x3)*y1)/1.2d+1
166  shp(3,3) = -((x7-x4)*y8+(-x8+x6-x4+x2)*y7+(x2-x7)*y6+(x8+x7-x2-x1
167  1 )*y4+(-x7-x6+x4+x1)*y2+(x4-x2)*y1)/1.2d+1
168  shp(3,4) = -((x7-x5+x3-x1)*y8+(x3-x8)*y7+(x8-x1)*y5+(-x8-x7+x2+x1
169  1 )*y3+(x1-x3)*y2+(x8+x5-x3-x2)*y1)/1.2d+1
170  shp(3,5) = ((x7+x6-x4-x1)*y8+(x6-x8)*y7+(-x8-x7+x2+x1)*y6+(x8-x1)
171  1 *y4+(x1-x6)*y2+(x8-x6+x4-x2)*y1)/1.2d+1
172  shp(3,6) = ((x7-x5)*y8+(-x8-x5+x3+x2)*y7+(x8+x7-x2-x1)*y5+(x2-x7)
173  1 *y3+(-x7+x5-x3+x1)*y2+(x5-x2)*y1)/1.2d+1
174  shp(3,7) = -((x6+x5-x4-x3)*y8+(-x8-x5+x3+x2)*y6+(x6-x8)*y5+(x8-x3
175  1 )*y4+(x8-x6+x4-x2)*y3+(x3-x6)*y2)/1.2d+1
176  shp(3,8) = ((x6+x5-x4-x3)*y7+(x5-x7)*y6+(-x7-x6+x4+x1)*y5+(x7-x5+
177  1 x3-x1)*y4+(x7-x4)*y3+(x4-x5)*y1)/1.2d+1
178 !
179 ! computation of element volume (eqn 15)
180 !
181  vol=0.0d0
182  do k=1,8
183  vol=vol+xl(1,k)*shp(1,k)
184  enddo
185 !
186 ! computation of the average jacobian determinant
187 !
188  xsj=vol/8.0d0
189 !
190 ! hourglass control vectors(see appendix 2 from the reference above).
191 ! divide shp array by element volume
192  a=0.0d0
193  do i=1,8
194  do j=1,3
195  a=a+shp(j,i)*shp(j,i)
196  s0=shp(j,i)/vol
197  shp(j,i)=s0
198  do k=1,4
199  gs(i,k)=gb(i,k)
200  do l=1,8
201  gs(i,k)=gs(i,k)-s0*xl(j,l)*gb(l,k)
202  enddo
203  enddo
204  enddo
205  enddo
206 c
207 c calculate hourglass control stiffness factor a
208 c (to be used in hgstiffness() and hgforce())
209 c
210 c in ABAQUS, a 0.005 is used as default value.
211  a=0.0005*a/vol
212 c
213  return
z1
static double * z1
Definition: filtermain.c:48
x1
static double * x1
Definition: filtermain.c:48
Hosted by OpenAircraft.com, (Michigan UAV, LLC)