39
40
41
43
44
45
46#include "implicit_f.inc"
47
48
49
50#include "mvsiz_p.inc"
51
52
53
54#include "param_c.inc"
55
56
57
58 INTEGER NEL,NGL(*)
60 . vol(*), geo(npropg,*),vzl(*), deltax(*), det(*)
62 . x1(*), x2(*), x3(*), x4(*), x5(*), x6(*),
63 . y1(*), y2(*), y3(*), y4(*), y5(*), y6(*),
64 . z1(*), z2(*), z3(*), z4(*), z5(*), z6(*)
65
66
67
68 INTEGER I, NFAC
70 . x21(mvsiz) , x31(mvsiz) , x41(mvsiz) , x54(mvsiz),x64(mvsiz),
71 . y21(mvsiz) , y31(mvsiz) , y41(mvsiz) , y54(mvsiz),y64(mvsiz),
72 . z21(mvsiz) , z31(mvsiz) , z41(mvsiz) , z54(mvsiz),z64(mvsiz),
73 . jac1(mvsiz), jac2(mvsiz) ,jac3(mvsiz),
74 . jac4(mvsiz), jac5(mvsiz) ,jac6(mvsiz),
75 . jac7(mvsiz), jac8(mvsiz) ,jac9(mvsiz),
76 . jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz)
77
79 . xioff(mvsiz), aream(mvsiz),
80 . atest(mvsiz),
area(6,mvsiz)
81
82
83
84 DO i=1,nel
85 x21(i)=x2(i)-x1(i)
86 x31(i)=x3(i)-x1(i)
87 x41(i)=x4(i)-x1(i)
88 x54(i)=x5(i)-x4(i)
89 x64(i)=x6(i)-x4(i)
90
91 y21(i)=y2(i)-y1(i)
92 y31(i)=y3(i)-y1(i)
93 y41(i)=y4(i)-y1(i)
94 y54(i)=y5(i)-y4(i)
95 y64(i)=y6(i)-y4(i)
96
97 z21(i)=z2(i)-z1(i)
98 z31(i)=z3(i)-z1(i)
99 z41(i)=z4(i)-z1(i)
100 z54(i)=z5(i)-z4(i)
101 z64(i)=z6(i)-z4(i)
102 ENDDO
103
104 DO i=1,nel
105
106 jac1(i)=x21(i)+x54(i)
107 jac2(i)=y21(i)+y54(i)
108 jac3(i)=z21(i)+z54(i)
109 ENDDO
110
111 DO i=1,nel
112
113 jac4(i)=x31(i)+x64(i)
114 jac5(i)=y31(i)+y64(i)
115 jac6(i)=z31(i)+z64(i)
116
117 jac7(i)=third*(x41(i)+x5(i)-x2(i)+x6(i)-x3(i))
118 jac8(i)=third*(y41(i)+y5(i)-y2(i)+y6(i)-y3(i))
119 jac9(i)=third*(z41(i)+z5(i)-z2(i)+z6(i)-z3(i))
120 ENDDO
121
122 DO i=1,nel
123 jac_59_68(i)=jac5(i)*jac9(i)-jac6(i)*jac8(i)
124 jac_67_49(i)=jac6(i)*jac7(i)-jac4(i)*jac9(i)
125 jac_48_57(i)=jac4(i)*jac8(i)-jac5(i)*jac7(i)
126 ENDDO
127
128 DO i=1,nel
129 det(i)=one_over_8*(jac1(i)*jac_59_68(i)+jac2(i)*jac_67_49(i)+jac3(i)*jac_48_57(i))
130 vol(i)=det(i)
131 ENDDO
132
133 DO i=1,nel
134 IF(det(i)>zero) cycle
136 . msgtype=msgerror,
137 . anmode=aninfo,
138 . i1=ngl(i))
139 ENDDO
140
141 DO i=1,nel
142 vzl(i) = fourth*(jac9(i)*(x54(i)*y64(i)-x21(i)*y31(i)-x64(i)*y54(i)+x31(i)*y21(i))
143 . -jac8(i)*(x54(i)*z64(i)+x31(i)*z21(i)-x21(i)*z31(i)-x64(i)*z54(i))
144 . +jac7(i)*(y54(i)*z64(i)+y31(i)*z21(i)-y21(i)*z31(i)-y64(i)*z54(i)))
145
146 ENDDO
147 DO i=1,nel
148 xioff(i)=one
149 aream(i)=zero
150 ENDDO
151
152 CALL slen(x1,x2,x5,x4,y1,y2,y4,y5,z1,z2,z5,z4,1,
area, aream)
153 CALL slen(x2,x5,x6,x3,y2,y5,y6,y3,z2,z5,z6,z3,2,
area, aream)
154 CALL slen(x1,x4,x6,x3,y1,y4,y6,y3,z1,z4,z6,z3,3,
area, aream)
155 CALL slen(x1,x2,x3,x3,y1,y2,y3,y3,z1,z2,z3,z3,4,
area, aream)
156 CALL slen(x4,x5,x6,x6,y4,y5,y6,y6,z4,z5,z6,z6,5,
area, aream)
157
158 DO i=1,nel
159 atest(i)=em4*aream(i)
160 nfac=0
161 IF(
area(1,i)<atest(i)) nfac=nfac+1
162 IF(
area(2,i)<atest(i)) nfac=nfac+1
163 IF(
area(3,i)<atest(i)) nfac=nfac+1
164 IF(
area(4,i)<atest(i)) nfac=nfac+1
165 IF(
area(5,i)<atest(i)) nfac=nfac+1
166 IF(nfac>=2) xioff(i)=ep03
167 ENDDO
168 DO i=1,nel
169 deltax(i)=four*det(i)*xioff(i)/sqrt(aream(i))
170 ENDDO
171
172 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine slen(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, j, area, aream)
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)