38
39
40
41#include "implicit_f.inc"
42
43
44
45#include "mvsiz_p.inc"
46
47
48
49#include "scr17_c.inc"
50
51
52
53 INTEGER, INTENT(IN) :: NEL,ICS
54 my_real,
DIMENSION(MVSIZ) ,
INTENT(OUT) ::
area,llsh
55 my_real,
DIMENSION(MVSIZ) ,
INTENT(IN) :: voln,
56 . x1, x2, x3, x4, x5, x6, x7, x8,
57 . y1, y2, y3, y4, y5, y6, y7, y8,
58 . z1, z2, z3, z4, z5, z6, z7, z8
59
60
61
62 INTEGER I, J, N
64 . rx(mvsiz),ry(mvsiz),rz(mvsiz),sx(mvsiz),sy(mvsiz),sz(mvsiz),
65 . vq(3,3,mvsiz), lxyz0(3),deta1(mvsiz),xx,yy,zz,
66 . xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),yl2(mvsiz),
67 . yl3(mvsiz),yl4(mvsiz),zl1(mvsiz),
68 . xn(mvsiz,4) , yn(mvsiz,4) , zn(mvsiz,4)
70 . al1,al2,ll(mvsiz),corel(2,4)
72 . x13,x24,y13,y24,l13,l24,c1,c2,thkly,posly,
73 . fac,visce,rx1,ry1,sx1,sy1,s1,fac1,fac2,faci,fac11,facdt
74
75 SELECT CASE(ics)
76 CASE (1)
77 DO i=1,nel
78 xn(i,1) = half*(x1(i)+x4(i))
79 yn(i,1) = half*(y1(i)+y4(i))
80 zn(i,1) = half*(z1(i)+z4(i))
81 xn(i,2) = half*(x2(i)+x3(i))
82 yn(i,2) = half*(y2(i)+y3(i))
83 zn(i,2) = half*(z2(i)+z3(i))
84 xn(i,3) = half*(x6(i)+x7(i))
85 yn(i,3) = half*(y6(i)+y7(i))
86 zn(i,3) = half*(z6(i)+z7(i))
87 xn(i,4) = half*(x5(i)+x8(i))
88 yn(i,4) = half*(y5(i)+y8(i))
89 zn(i,4) = half*(z5(i)+z8(i))
90 ENDDO
91 CASE (10)
92 DO i=1,nel
93 xn(i,1) = half*(x1(i)+x5(i))
94 yn(i,1) = half*(y1(i)+y5(i))
95 zn(i,1) = half*(z1(i)+z5(i))
96 xn(i,2) = half*(x2(i)+x6(i))
97 yn(i,2) = half*(y2(i)+y6(i))
98 zn(i,2) = half*(z2(i)+z6(i))
99 xn(i,3) = half*(x3(i)+x7(i))
100 yn(i,3) = half*(y3(i)+y7(i))
101 zn(i,3) = half*(z3(i)+z7(i))
102 xn(i,4) = half*(x4(i)+x8(i))
103 yn(i,4) = half*(y4(i)+y8(i))
104 zn(i,4) = half*(z4(i)+z8(i))
105 ENDDO
106 CASE (100)
107 DO i=1,nel
108 xn(i,1) = half*(x1(i)+x2(i))
109 yn(i,1) = half*(y1(i)+y2(i))
110 zn(i,1) = half*(z1(i)+z2(i))
111 xn(i,2) = half*(x5(i)+x6(i))
112 yn(i,2) = half*(y5(i)+y6(i))
113 zn(i,2) = half*(z5(i)+z6(i))
114 xn(i,3) = half*(x8(i)+x7(i))
115 yn(i,3) = half*(y8(i)+y7(i))
116 zn(i,3) = half*(z8(i)+z7(i))
117 xn(i,4) = half*(x4(i)+x3(i))
118 yn(i,4) = half*(y4(i)+y3(i))
119 zn(i,4) = half*(z4(i)+z3(i))
120 ENDDO
121 END SELECT
122
123 DO i=1,nel
124 rx(i)=xn(i,2)+xn(i,3)-xn(i,1)-xn(i,4)
125 ry(i)=yn(i,2)+yn(i,3)-yn(i,1)-yn(i,4)
126 rz(i)=zn(i,2)+zn(i,3)-zn(i,1)-zn(i,4)
127 sx(i)=xn(i,3)+xn(i,4)-xn(i,1)-xn(i,2)
128 sy(i)=yn(i,3)+yn(i,4)-yn(i,1)-yn(i,2)
129 sz(i)=zn(i,3)+zn(i,4)-zn(i,1)-zn(i,2)
130 ENDDO
131
132 CALL clsys3(rx, ry, rz, sx, sy, sz,
133 . vq, deta1,nel)
134
135 DO i=1,nel
136 lxyz0(1)=fourth*(xn(i,1)+xn(i,2)+xn(i,3)+xn(i,4))
137 lxyz0(2)=fourth*(yn(i,1)+yn(i,2)+yn(i,3)+yn(i,4))
138 lxyz0(3)=fourth*(zn(i,1)+zn(i,2)+zn(i,3)+zn(i,4))
139 xx=xn(i,2)-xn(i,1)
140 yy=yn(i,2)-yn(i,1)
141 zz=zn(i,2)-zn(i
142 xl2(i)=vq(1,1,i)*xx+vq(2,1,i)*yy+vq(3,1,i)*zz
143 yl2(i)=vq(1,2,i)*xx+vq(2,2,i)*yy+vq(3,2,i)*zz
144 xx=xn(i,2)-lxyz0(1)
145 yy=yn(i,2)-lxyz0(2)
146 zz=zn(i,2)-lxyz0(3)
147 zl1(i)=vq(1,3,i)*xx+vq(2,3,i)*yy+vq(3,3,i)*zz
148
149 xx=xn(i,3)-xn(i,1)
150 yy=yn(i,3)-yn(i,1)
151 zz=zn(i,3)-zn(i,1)
152 xl3(i)=vq(1,1,i)*xx+vq(2,1,i)*yy+vq(3,1,i)*zz
153 yl3(i)=vq(1,2,i)*xx+vq(2,2,i)*yy+vq(3,2,i)*zz
154
155 xx=xn(i,4)-xn(i,1)
156 yy=yn(i,4)-yn(i,1)
157 zz=zn(i,4)-zn(i,1)
158 xl4(i)=vq(1,1,i)*xx+vq(2,1,i)*yy+vq(3,1,i)*zz
159 yl4(i)=vq(1,2,i)*xx+vq(2,2,i)*yy+vq(3,2,i)*zz
160 area(i)=fourth*deta1(i)
161 ENDDO
162 fac = two
163 facdt = five_over_4
164
165 IF (idt1sol>0) facdt =four_over_3
166
167 DO i=1,nel
168 lxyz0(1)=fourth*(xl2(i)+xl3(i)+xl4(i))
169 lxyz0(2)=fourth*(yl2(i)+yl3(i)+yl4(i))
170 corel(1,1)=-lxyz0(1)
171 corel(1,2)=xl2(i)-lxyz0(1)
172 corel(1,3)=xl3(i)-lxyz0(1)
173 corel(1,4)=xl4(i)-lxyz0(1)
174 corel(2,1)=-lxyz0(2)
175 corel(2,2)=yl2(i)-lxyz0(2)
176 corel(2,3)=yl3(i)-lxyz0(2)
177 corel(2,4)=yl4(i)-lxyz0(2)
178 x13=(corel(1,1)-corel(1,3))*half
179 x24=(corel(1,2)-corel(1,4))*half
180 y13=(corel(2,1)-corel(2,3))*half
181 y24=(corel(2,2)-corel(2,4))*half
182
183 l13=x13*x13+y13*y13
184 l24=x24*x24+y24*y24
186 c1 =corel(1,2)*corel(2,4)-corel(2,2)*corel(1,4)
187 c2 =corel(1,1)*corel(2,3)-corel(2,1)*corel(1,3)
188 al2 =
max(abs(c1),abs(c2))/
area(i)
189 rx1=x24-x13
190 ry1=y24-y13
191 sx1=-x24-x13
192 sy1=-y24-y13
193 c1=sqrt(rx1*rx1+ry1*ry1)
194 c2=sqrt(sx1*sx1+sy1*sy1)
195 s1=fourth*(
max(c1,c2)/
min(c1,c2)-one)
196 fac1=
min(half,s1)+one
198 fac2=3.413*
max(zero,fac2-0.7071)
199 fac2=0.78+0.22*fac2*fac2*fac2
200 faci=two*fac1*fac2
201 s1 = sqrt(faci*(facdt+al2)*al1)
204 ENDDO
205
206 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine clsys3(rx, ry, rz, sx, sy, sz, vq, det, nel, mvsiz)