44
45
46
48
49
50
51#include "implicit_f.inc"
52#include "comlock.inc"
53
54
55
56#include "mvsiz_p.inc"
57
58
59
60
61
62
63 INTEGER, INTENT(IN) :: NEL
64
66 . x(mvsiz,8), y(mvsiz,8), z(mvsiz,8),
67 . px1(mvsiz,8),px2(mvsiz,8),px3(mvsiz,8),px4(mvsiz,8),
68 . px5(mvsiz,8),px6(mvsiz,8),px7(mvsiz,8),px8(mvsiz,8),
69 . py1(mvsiz,8),py2(mvsiz,8),py3(mvsiz,8),py4(mvsiz,8),
70 . py5(mvsiz,8),py6(mvsiz,8),py7(mvsiz,8),py8(mvsiz,8),
71 . pz1(mvsiz,8),pz2(mvsiz,8),pz3(mvsiz,8),pz4(mvsiz,8),
72 . pz5(mvsiz,8),pz6(mvsiz,8),pz7(mvsiz,8),pz8(mvsiz,8),
73 . vlinc(mvsiz,8),vol(mvsiz),deltax(mvsiz) ,off(mvsiz)
74 INTEGER NGL(MVSIZ)
75
76
77
78 INTEGER I, IPT
79
81 . aj11(mvsiz), aj12(mvsiz), aj13(mvsiz), aj21(mvsiz),
82 . aj22(mvsiz), aj23(mvsiz), aj31(mvsiz), aj32(mvsiz),
83 . aj33(mvsiz), ai11(mvsiz), ai12, ai13,
84 . ai21(mvsiz), ai22, ai23, ai31(mvsiz),
85 . ai32, ai33,
86 . a1111, a1113,
87 . a1115, a1117, a2111, a2113,
88 . a2115, a2117, a3111, a3113,
89 . a3115, a3117, a1221, a1222,
90 . a1225, a1226, a2221, a2222,
91 . a2225, a2226, a3221, a3222
92
94 . a3225, a3226, a1331, a1332,
95 . a1333, a1334, a2331, a2332,
96 . a2333, a2334, a3331, a3332,
97 . a3333, a3334, x12(mvsiz), x34(mvsiz), x56(mvsiz),
98 . x78(mvsiz), y12(mvsiz), y34(mvsiz), y56(mvsiz), y78(mvsiz),
99 . z12(mvsiz), z34(mvsiz), z56(mvsiz), z78(mvsiz), x14(mvsiz),
100 . x23(mvsiz), x58(mvsiz), x67(mvsiz), y14(mvsiz), y23(mvsiz),
101 . y58(mvsiz
102 . z67(mvsiz), x15(mvsiz), x26(mvsiz), x37(mvsiz), x48(mvsiz),
103 . y15(mvsiz), y26(mvsiz), y37(mvsiz), y48(mvsiz), z15(mvsiz),
104 . z26(mvsiz), z37(mvsiz), z48(mvsiz), h(8), vlinv,
105 . delta1,delta2,delta3,
106 . spx1,spx2,spx3,spx4,spx5,spx6,spx7,spx8,
107 . spy1,spy2,spy3,spy4,spy5,spy6,spy7,spy8,
108 . spz1,spz2,spz3,spz4,spz5,spz6,spz7,spz8,
109 . xx1,yy1,zz1,xx2,yy2,zz2,xx3,yy3,zz3,smax,vmin,
110 . p1(8), p2(8), p3(8)
111 DO 5 i=1,nel
112 x12(i)=x(i,1)-x(i,2)
113 y12(i)=y(i,1)-y(i,2)
114 z12(i)=z(i,1)-z(i,2)
115 x34(i)=x(i,3)-x(i,4)
116 y34(i)=y(i,3)-y(i,4)
117 z34(i)=z(i,3)-z(i,4)
118 x56(i)=x(i,5)-x(i,6)
119 y56(i)=y(i,5)-y(i,6)
120 z56(i)=z(i,5)-z(i,6)
121 x78(i)=x(i,7)-x(i,8)
122 y78(i)=y(i,7)-y(i,8)
123 z78(i)=z(i,7)-z(i,8)
124 x14(i)=x(i,1)-x(i,4)
125 y14(i)=y(i,1)-y(i,4)
126 z14(i)=z(i,1)-z(i,4)
127 x23(i)=x(i,2)-x(i,3)
128 y23(i)=y(i,2)-y(i,3)
129 z23(i)=z(i,2)-z(i,3)
130 x58(i)=x(i,5)-x(i,8)
131 y58(i)=y(i,5)-y(i,8)
132 z58(i)=z(i,5)-z(i,8)
133 x67(i)=x(i,6)-x(i,7)
134 y67(i)=y(i,6)-y(i,7)
135 z67(i)=z(i,6)-z(i,7)
136 x15(i)=x(i,1)-x(i,5)
137 y15(i)=y(i,1)-y(i,5)
138 z15(i)=z(i,1)-z(i,5)
139 x26(i)=x(i,2)-x(i,6)
140 y26(i)=y(i,2)-y(i,6)
141 z26(i)=z(i,2)-z(i,6)
142 x37(i)=x(i,3)-x(i,7)
143 y37(i)=y(i,3)-y(i,7)
144 z37(i)=z(i,3)-z(i,7)
145 x48(i)=x(i,4)-x(i,8)
146 y48(i)=y(i,4)-y(i,8)
147 z48(i)=z(i,4)-z(i,8)
148
149 vol(i)=0.
150 5 CONTINUE
151
152 DO 100 ipt=1,8
153
154
155
156 CALL basisf (h,p1,p2,p3,ipt)
157
158 DO i=1,nel
159 aj11(i)=p1(1)*x12(i)+p1(3)*x34(i)+p1(5)*x56(i)+p1(7)*x78(i)
160 aj12(i)=p1(1)*y12(i)+p1(3)*y34(i)+p1(5)*y56(i)+p1(7)*y78(i)
161 aj13(i)=p1(1)*z12(i)+p1(3)*z34(i)+p1(5)*z56(i)+p1(7)*z78(i)
162 aj21(i)=p2(1)*x14(i)+p2(2)*x23(i)+p2(5)*x58(i)+p2(6)*x67(i)
163 aj22(i)=p2(1)*y14(i)+p2(2)*y23(i)+p2(5)*y58(i)+p2(6)*y67(i)
164 aj23(i)=p2(1)*z14(i)+p2(2)*z23(i)+p2(5)*z58(i)+p2(6)*z67(i)
165 aj31(i)=p3(1)*x15(i)+p3(2)*x26(i)+p3(3)*x37(i)+p3(4)*x48(i)
166 aj32(i)=p3(1)*y15(i)+p3(2)*y26(i)+p3(3)*y37(i)+p3(4)*y48(i)
167 aj33(i)=p3(1)*z15(i)+p3(2)*z26(i)+p3(3)*z37(i)+p3(4)*z48(i)
168 ENDDO
169
170 DO 20 i=1,nel
171 ai11(i)= aj22(i)*aj33(i)-aj23(i)*aj32(i)
172 ai21(i)=-aj21(i)*aj33(i)+aj23(i)*aj31(i)
173 ai31(i)= aj21(i)*aj32(i)-aj22(i)*aj31(i)
174 vlinc(i,ipt)=aj11(i)*ai11(i)+aj12(i)*ai21(i)+aj13(i)*ai31(i)
175 20 CONTINUE
176
177 DO 130 i=1,nel
178
179 IF(off(i)==zero)vlinc(i,ipt)=one
180 vol(i)=vol(i)+vlinc(i,ipt)
181 IF(vlinc(i,ipt)>zero) GO TO 130
182 CALL ancmsg(msgid=172,anmode=aninfo,
183 . i1=ngl(i),i2=ipt)
185 130 CONTINUE
186
187 DO i=1,nel
188 vlinv=one/vlinc(i,ipt)
189 ai11(i)=vlinv*ai11(i)
190 ai21(i)=vlinv*ai21(i)
191 ai31(i)=vlinv*ai31(i)
192 ai12=vlinv*(-aj12(i)*aj33(i)+aj13(i)*aj32(i))
193 ai22=vlinv*( aj11(i)*aj33(i)-aj13(i)*aj31(i))
194 ai32=vlinv*(-aj11(i)*aj32(i)+aj12(i)*aj31(i))
195 ai13=vlinv*( aj12(i)*aj23(i)-aj13(i)*aj22(i))
196 ai23=vlinv*(-aj11(i)*aj23(i)+aj13(i)*aj21(i))
197 ai33=vlinv*( aj11(i)*aj22(i)-aj12(i)*aj21(i))
198
199 a1111=ai11(i)*p1(1)
200 a1113=ai11(i)*p1(3)
201 a1115=ai11(i)*p1(5)
202 a1117=ai11(i)*p1(7)
203 a2111=ai21(i)*p1(1)
204 a2113=ai21(i)*p1(3)
205 a2115=ai21(i)*p1(5)
206 a2117=ai21(i)*p1(7)
207 a3111=ai31(i)*p1(1)
208 a3113=ai31(i)*p1(3)
209 a3115=ai31(i)*p1(5)
210 a3117=ai31(i)*p1(7)
211 a1221=ai12*p2(1)
212 a1222=ai12*p2(2)
213 a1225=ai12*p2(5)
214 a1226=ai12*p2(6)
215 a2221=ai22*p2(1)
216 a2222=ai22*p2(2)
217
218 a2225=ai22*p2(5)
219 a2226=ai22*p2(6)
220 a3221=ai32*p2(1)
221 a3222=ai32*p2(2)
222 a3225=ai32*p2(5)
223 a3226=ai32*p2(6)
224 a1331=ai13*p3(1)
225 a1332=ai13*p3(2)
226 a1333=ai13*p3(3)
227 a1334=ai13*p3(4)
228 a2331=ai23*p3(1)
229 a2332=ai23*p3(2)
230 a2333=ai23*p3(3)
231 a2334=ai23*p3(4)
232 a3331=ai33*p3(1)
233 a3332=ai33*p3(2)
234 a3333=ai33*p3(3)
235 a3334=ai33*p3(4)
236 px1(i,ipt)= a1111+a1221+a1331
237 px2(i,ipt)=-a1111+a1222+a1332
238 px3(i,ipt)= a1113-a1222+a1333
239 px4(i,ipt)=-a1113-a1221+a1334
240 px5(i,ipt)= a1115+a1225-a1331
241 px6(i,ipt)=-a1115+a1226-a1332
242 px7(i,ipt)= a1117-a1226-a1333
243 px8(i,ipt)=-a1117-a1225-a1334
244 py1(i,ipt)= a2111+a2221+a2331
245 py2(i,ipt)=-a2111+a2222+a2332
246 py3(i,ipt)= a2113-a2222+a2333
247 py4(i,ipt)=-a2113-a2221+a2334
248 py5(i,ipt)= a2115+a2225-a2331
249 py6(i,ipt)=-a2115+a2226-a2332
250 py7(i,ipt)= a2117-a2226-a2333
251 py8(i,ipt)=-a2117-a2225-a2334
252 pz1(i,ipt)= a3111+a3221+a3331
253 pz2(i,ipt)=-a3111+a3222+a3332
254 pz3(i,ipt)= a3113-a3222+a3333
255 pz4(i,ipt)=-a3113-a3221+a3334
256 pz5(i,ipt)= a3115+a3225-a3331
257 pz6(i,ipt)=-a3115+a3226-a3332
258 pz7(i,ipt)= a3117-a3226-a3333
259 pz8(i,ipt)=-a3117-a3225-a3334
260 ENDDO
261
262 100 CONTINUE
263
264 DO i=1,nel
265
266 xx1 = x(i,1) + x(i,2) + x(i,3) + x(i,4)
267 . - x(i,5) - x(i,6) - x(i,7) - x(i,8)
268 yy1 = y(i,1) + y(i,2) + y(i,3) + y(i,4)
269 . - y(i,5) - y(i,6) - y(i,7) - y(i,8)
270 zz1 = z(i,1) + z(i,2) + z(i,3) + z(i,4)
271 . - z(i,5) - z(i,6) - z(i,7) - z(i,8)
272 xx2 = x(i,1) + x(i,2) + x(i,5) + x(i,6)
273 . - x(i,3) - x(i,4) - x(i,7) - x(i,8)
274 yy2 = y(i,1) + y(i,2) + y(i,5) + y(i,6)
275 . - y(i,3) - y(i,4) - y(i,7) - y(i,8)
276 zz2 = z(i,1) + z(i,2) + z(i,5) + z(i,6)
277 . - z(i,3) - z(i,4) - z(i,7) - z(i,8)
278 xx3 = x(i,1) + x(i,4) + x(i,5) + x(i,8)
279 . - x(i,3) - x(i,2) - x(i,7) - x(i,6)
280 yy3 = y(i,1) + y(i,4) + y(i,5) + y(i,8)
281 . - y(i,3) - y(i,2) - y(i,7) - y(i,6)
282 zz3 = z(i,1) + z(i,4) + z(i,5) + z(i,8)
283 . - z(i,3) - z(i,2) - z(i,7) - z(i,6)
284
285 smax = (yy1 * zz2 - yy2 * zz1)**2
286 . + (zz1 * xx2 - zz2 * xx1)**2
287 . + (xx1 * yy2 - xx2 * yy1)**2
288 smax =
max(smax,(yy1 * zz3 - yy3 * zz1)**2
289 . + (zz1 * xx3 - zz3 * xx1)**2
290 . + (xx1 * yy3 - xx3 * yy1)**2)
291 smax =
max(smax,(yy3 * zz2 - yy2 * zz3)**2
292 . + (zz3 * xx2 - zz2 * xx3)**2
293 . + (xx3 * yy2 - xx2 * yy3)**2)
294
295 vmin =
min(vlinc(i,1),vlinc(i,2),vlinc(i,3),vlinc(i,4),
296 . vlinc(i,5),vlinc(i,6),vlinc(i,7),vlinc(i,8))
297 deltax(i)=hundred28*vmin/sqrt(smax)
298 ENDDO
299 RETURN
subroutine basisf(h, p1, p2, p3, ipt)
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
1.15.0