32
33
34
35#include "implicit_f.inc"
36
37
38
39#include "mvsiz_p.inc"
40
41
42
43 INTEGER JFT, JLT,NVS,IVS(*),IXTG1(4,*)
45 . x2(*), y2(*), x3(*), y3(*),area2(*),
46 . x4(*), y4(*), x5(*), y5(*),area4(*),area5(*),
47 . x6(*), y6(*), z4(*), z5(*),z6(*),area6(*),
48 . n4x(*),n4y(*),n4z(*),n5x(*),n5y(*),n5z(*),n6x(*),n6y(*),
49 . n6z(*),px2(*), py2(*), px3(*), py3(*), nu(*),thk2(*),
50 . aldt(*),alpe(*),pb1(mvsiz,3,3),pb2(mvsiz,3,3),pb3(mvsiz,3,6)
51
52
53
54 INTEGER I, J,EP
56 . areai(mvsiz),
57 . al1(mvsiz), al2(mvsiz), al3(mvsiz), almax, almin
59 . x32(mvsiz), y32(mvsiz),al4, al5, al6,ah(mvsiz),
60 . h1(mvsiz),h2(mvsiz),h3(mvsiz),h4(mvsiz),h5(mvsiz),h6(mvsiz),
61 . c1(mvsiz),c2(mvsiz),c3(mvsiz),s1(mvsiz),s2(mvsiz),s3(mvsiz),
62 . h14(mvsiz),h24(mvsiz),h25(mvsiz),h35(mvsiz),h4d,h5d,h6d,
63 . beta1(mvsiz),beta2(mvsiz),beta3(mvsiz),x42,y42,z2,x53,y53,
64 . h36(mvsiz),h16(mvsiz),al41,al42,al52,al53,ai1,ai2,ai3,
65 . beta14(mvsiz),beta24(mvsiz),beta25(mvsiz),beta35(mvsiz),
66 . beta36(mvsiz),beta16(mvsiz),al61,al63,x61,y61
68 . c11(mvsiz,3,3), c21(mvsiz,3,3),c22(mvsiz,3),b11,b12,b13,
69 . b21,b22,b23,b31,b32,b33,b1,b2,b3,h15,h26,h34,fac
70
71
72
73
74 DO i=jft,jlt
75 areai(i)=one/area2(i)
76 px2(i)=y3(i)*areai(i)
77 py2(i)=-x3(i)*areai(i)
78 px3(i)=-y2(i)*areai(i)
79 py3(i)=x2(i)*areai(i)
80 x32(i) =x3(i)-x2(i)
81 y32(i) =y3(i)-y2(i)
82 ENDDO
83
84 DO i=jft,jlt
85 al1(i) = sqrt(x32(i)*x32(i) + y32(i)*y32(i))
86 al2(i) = sqrt(x3(i) * x3(i) + y3(i) * y3(i))
87 al3(i) = sqrt(x2(i) * x2(i) + y2(i) * y2(i))
88
89 h1(i) = al1(i)*areai(i)
90 h2(i) = al2(i)*areai(i)
91 h3(i) = al3(i)*areai(i)
92 ai1 = one/al1(i)
93 ai2 = one/al2(i)
94 ai3 = one/al3(i)
95 c1(i) = y32(i)*ai1
96 s1(i) = -x32(i)*ai1
97 c2(i) = -y3(i)*ai2
98 s2(i) = x3(i)*ai2
99 c3(i) = y2(i)*ai3
100 s3(i) = -x2(i)*ai3
101 beta1(i) = s3(i)*s2(i)+c3(i)*c2(i)
102 beta2(i) = s3(i)*s1(i)+c3(i)*c1(i)
103 beta3(i) = s1(i)*s2(i)+c1(i)*c2(i)
104 z2 =z4(i)*z4(i)
105 al41 = x4(i) * x4(i) + y4(i) * y4(i) +z2
106 x42 =x4(i)-x2(i)
107 y42 =y4(i)-y2(i)
108 al42 = x42 * x42 + y42 * y42+z2
109 h4(i) = al3(i)*area4(i)
110 h14(i) = area4(i)*sqrt(al42)
111 h24(i) = area4(i)*sqrt(al41)
112 h4d= (al42-al41)/al3(i)
113 beta14(i) = half*(h4d-al3(i))/sqrt(al41)
114 beta24(i) = -half*(h4d+al3(i))/sqrt(al42)
115 z2 =z5(i)*z5(i)
116 al52 = x5(i) * x5(i) + y5(i) * y5(i)+z2
117 x53 =x5(i)-x32(i)
118 y53 =y5(i)-y32(i)
119 al53 = x53 * x53 + y53 * y53+z2
120 h5(i)= al1(i)*area5(i)
121 h25(i) = area5(i)*sqrt(al53)
122 h35(i) = area5(i)*sqrt(al52)
123 h5d= (al53-al52)/al1(i)
124 beta25(i) = half*(h5d-al1(i))/sqrt(al52)
125 beta35(i) = -half*(h5d+al1(i))/sqrt(al53)
126 z2 =z6(i)*z6(i)
127 al63 = x6(i) * x6(i) + y6(i) * y6(i)+z2
128 x61 =x6(i)+x3(i)
129 y61 =y6(i)+y3(i)
130 al61 = x61 * x61 + y61 * y61+z2
131 h6(i) = al2(i)*area6(i)
132 h36(i) = area6(i)*sqrt(al61)
133 h16(i) = area6(i)*sqrt(al63)
134 h6d= (al61-al63)/al2(i)
135 beta36(i) = half*(h6d-al2(i))/sqrt(al63)
136 beta16(i) = -half*(h6d+al2(i))/sqrt(al61)
137 ENDDO
138 DO i=jft,jlt
139 c11(i,1,1) = h1(i)
140 c11(i,2,1) = h1(i)*beta3(i)
141 c11(i,3,1) = h1(i)*beta2(i)
142 c11(i,1,2) = h2(i)*beta3(i)
143 c11(i,2,2) = h2
144 c11(i,3,2) = h2(i)*beta1(i)
145 c11(i,1,3) = h3(i)*beta2(i)
146 c11(i,2,3) = h3(i)*beta1(i)
147 c11(i,3,3) = h3(i)
148 ENDDO
149 DO i=jft,jlt
150 c21(i,1,1) = beta24(i)*h14(i)
151 c21(i,2,1) = 0.0
152 c21(i,3,1) = beta36(i)*h16(i)
153 c21(i,1,2) = beta14(i)*h24(i)
154 c21(i,2,2) = beta35(i)*h25(i)
155 c21(i,3,2) = 0.0
156 c21(i,1,3) = 0.0
157 c21(i,2,3) = beta25(i)*h35(i)
158 c21(i,3,3) = beta16(i)*h36(i)
159 ENDDO
160 DO i=jft,jlt
161 c22(i,1) = h4(i)
162 c22(i,2) = h5(i)
163 c22(i,3) = h6(i)
164 ENDDO
165
166 DO ep=nvs+1,jlt
167 i =ivs(ep)
168
169
170 IF (ixtg1(1,i)==0) THEN
171 c21(i,1,1) = -c11(i,3,1)
172 c21(i,1,2) = -c11(i,3,2)
173 c21(i,1,3) = -c11(i,3,3)
174 c22(i,1) = 0.0
175 h4(i)=ep20
176
177 ELSEIF (ixtg1(1,i)<0) THEN
178 c21(i,1,1) = c11(i,3,1)
179 c21(i,1,2) = c11(i,3,2)
180 c21(i,1,3) = c11(i,3,3)
181 c22(i,1) = 0.0
182 h4(i)=ep20
183 ENDIF
184
185
186 IF (ixtg1(2,i)==0) THEN
187 c21(i,2,1) = -c11(i,1,1)
188 c21(i,2,2) = -c11(i,1,2)
189 c21(i,2,3) = -c11(i,1,3)
190 c22(i,2) = 0.0
191 h5(i)=ep20
192
193 ELSEIF (ixtg1(2,i)<0) THEN
194 c21(i,2,1) = c11(i,1,1)
195 c21(i,2,2) = c11(i,1,2)
196 c21(i,2,3) = c11(i,1,3)
197 c22(i,2) = 0.0
198 h5(i)=ep20
199 ENDIF
200
201
202 IF (ixtg1(3,i)==0) THEN
203 c21(i,3,1) = -c11(i,2,1)
204 c21(i,3,2) = -c11(i,2,2)
205 c21(i,3,3) = -c11(i,2,3)
206 c22(i,3) = 0.0
207 h6(i)=ep20
208
209 ELSEIF (ixtg1(3,i)<0) THEN
210 c21(i,3,1) = c11(i,2,1)
211 c21(i,3,2) = c11(i,2,2)
212 c21(i,3,3) = c11(i,2,3)
213 c22(i,3) = 0.0
214 h6(i)=ep20
215 ENDIF
216 ENDDO
217
218 DO i=jft,jlt
219 h15 =-two*h1(i)*h5(i)/(h1(i)+h5(i))
220 h26 =-two*h2(i)*h6(i)/(h2(i)+h6(i))
221 h34 =-two*h3(i)*h4(i)/(h3(i)+h4(i))
222 almin = -half*
min(h15,h26,h34)
223 ah(i) = almin*almin*thk2(i)
224
225 b11=c1(i)*c1(i)*h15
226 b12=c2(i)*c2(i)*h26
227 b13=c3(i)*c3(i)*h34
228 b21= h15 - b11
229 b22= h26 - b12
230 b23= h34 - b13
231 b31=two*c1(i)*s1(i)*h15
232 b32=two*c2(i)*s2(i)*h26
233 b33=two*c3(i)*s3(i)*h34
234
235 pb1(i,1,1) =b11*c11(i,1,1)+b12*c11(i,2,1)+b13*c11(i,3,1)
236 pb1(i,1,2) =b11*c11(i,1,2)+b12*c11(i,2,2)+b13*c11(i,3,2)
237 pb1(i,1,3) =b11*c11(i,1,3)+b12*c11(i,2,3)+b13*c11(i,3,3)
238 pb1(i,2,1) =b21*c11(i,1,1)+b22*c11(i,2,1)+b23*c11(i,3,1)
239 pb1(i,2,2) =b21*c11(i,1,2)+b22*c11(i,2,2)+b23*c11(i,3,2)
240 pb1(i,2,3) =b21*c11(i,1,3)+b22*c11(i,2,3)+b23*c11(i,3,3)
241 pb1(i,3,1) =b31*c11(i,1,1)+b32*c11(i,2,1)+b33*c11(i,3,1)
242 pb1(i,3,2) =b31*c11(i,1,2)+b32*c11(i,2,2)+b33*c11(i,3,2)
243 pb1(i,3,3) =b31*c11(i,1,3)+b32*c11(i,2,3)+b33*c11(i,3,3)
244
245 pb2(i,1,1) =b13*c22(i,1)
246 pb2(i,1,2) =b11*c22(i,2)
247 pb2(i,1,3) =b12*c22(i,3)
248 pb2(i,2,1) =b23*c22(i,1)
249 pb2(i,2,2) =b21*c22(i,2)
250 pb2(i,2,3) =b22*c22(i,3)
251 pb2(i,3,1) =b33*c22(i,1)
252 pb2(i,3,2) =b31*c22(i,2)
253 pb2(i,3,3) =b32*c22(i,3)
254
255 pb3(i,1,1) =b11
256 pb3(i,1,2) =b12
257 pb3(i,1,3) =b13
258 pb3(i,2,1) =b21
259 pb3(i,2,2) =b22
260 pb3(i,2,3) =b23
261 pb3(i,3,1) =b31
262 pb3(i,3,2) =b32
263 pb3(i,3,3) =b33
264 pb3(i,1,4) =c21(i,2,1)
265 pb3(i,1,5) =c21(i,2,2)
266 pb3(i,1,6) =c21(i,2,3)
267 pb3(i,2,4) =c21(i,3,1)
268 pb3(i,2,5) =c21(i,3,2)
269 pb3(i,2,6) =c21(i,3,3)
270 pb3(i,3,4) =c21(i,1,1)
271 pb3(i,3,5) =c21(i,1,2)
272 pb3(i,3,6) =c21(i,1,3)
273 ENDDO
274
275 DO i=jft,jlt
276 fac =one+zep6*(1+nu(i))*ah(i)
277 almax = sqrt(fac)*
max(al1(i),al2(i),al3(i))
278 aldt(i)= area2(i) /almax
279 alpe(i)=one
280 ENDDO
281
282
283
284 RETURN
285