31
32
33
34#include "implicit_f.inc"
35
36
37
38 INTEGER IEL, JEL
40 . x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, xp, yp, zp,
41 . nrx, nry, nrz,
area, rvlh, rvlg
42
43
44
45 INTEGER I, J
47 . vx1, vy1, vz1, vx2, vy2, vz2, s1, s12, s2, nr1, nr2,
48 . x0, y0, z0, ksi(5), eta(5), dksi(4), deta(4), r(5),
49 . xls, yls, zls, s(4), v, fln, arg,
50 . d2, l12, l22, l32, l42, lm2
52
53 x0=fourth*(x1+x2+x3+x4)
54 y0=fourth*(y1+y2+y3+y4)
55 z0=fourth*(z1+z2+z3+z4)
56
57
58 d2 =(x0-xp)**2+(y0-yp)**2+(z0-zp)**2
59 l12=(x2-x1)**2+(y2-y1)**2+(z2-z1)**2
60 l22=(x3-x2)**2+(y3-y2)**2+(z3-z2)**2
61 l32=(x4-x3)**2+(y4-y3)**2+(z4-z3)**2
62 l42=(x1-x4)**2+(y1-y4)**2+(z1-z4)**2
63 lm2=
max(l12,l22,l32,l42)
64
65 IF(d2>twenty5*lm2) THEN
67 rvlh=
area*(nrx*(x0-xp)+nry*(y0-yp)+nrz*(z0-zp))/(d2**three_half)
68 ELSE
69
70 vx1=x2-x1
71 vy1=y2-y1
72 vz1=z2-z1
73 vx2=x4-x1
74 vy2=y4-y1
75 vz2=z4-z1
76
77 s1 =vx1*vx1+vy1*vy1+vz1*vz1
78 s12=vx1*vx2+vy1*vy2+vz1*vz2
79 nr1=sqrt(s1)
80
81 vx2=vx2-s12/s1*vx1
82 vy2=vy2-s12/s1*vy1
83 vz2=vz2-s12/s1*vz1
84
85 s2=vx2*vx2+vy2*vy2+vz2*vz2
86 nr2=sqrt(s2)
87 vx1=vx1/nr1
88 vy1=vy1/nr1
89 vz1=vz1/nr1
90 vx2=vx2/nr2
91 vy2=vy2/nr2
92 vz2=vz2/nr2
93
94 xls=(xp-x0)*vx1+(yp-y0)*vy1+(zp-z0)*vz1
95 yls=(xp-x0)*vx2+(yp-y0)*vy2+(zp-z0)*vz2
96 zls=(xp-x0)*nrx+(yp-y0)*nry+(zp-z0)*nrz
97
98 ksi(1)=(x1-x0)*vx1+(y1-y0)*vy1+(z1-z0)*vz1
99 eta(1)=(x1-x0)*vx2+(y1-y0)*vy2+(z1-z0)*vz2
100 ksi(2)=(x2-x0)*vx1+(y2-y0)*vy1+(z2-z0)*vz1
101 eta(2)=(x2-x0)*vx2+(y2-y0)*vy2+(z2-z0)*vz2
102 ksi(3)=(x3-x0)*vx1+(y3-y0)*vy1+(z3-z0)*vz1
103 eta(3)=(x3-x0)*vx2+(y3-y0)*vy2+(z3-z0)*vz2
104 ksi(4)=(x4-x0)*vx1+(y4-y0)*vy1+(z4-z0)*vz1
105 eta(4)=(x4-x0)*vx2+(y4-y0)*vy2+(z4-z0)*vz2
106 ksi(5)=ksi(1)
107 eta(5)=eta(1)
108
109 dksi(1)=ksi(2)-ksi(1)
110 dksi(2)=ksi(3)-ksi(2)
111 dksi(3)=ksi(4)-ksi(3)
112 dksi(4)=ksi(1)-ksi(4)
113 deta(1)=eta(2)-eta(1)
114 deta(2)=eta(3)-eta(2)
115 deta(3)=eta(4)-eta(3)
116 deta(4)=eta(1)-eta(4)
117 r(1)=sqrt((xp-x1)**2+(yp-y1)**2+(zp-z1)**2)
118 r(2)=sqrt((xp-x2)**2+(yp-y2)**2+(zp-z2)**2)
119 r(3)=sqrt((xp-x3)**2+(yp-y3)**2+(zp-z3)**2)
120 r(4)=sqrt((xp-x4)**2+(yp-y4)**2+(zp-z4)**2)
121 s(1)=sqrt(l12)
122 s(2)=sqrt(l22)
123 s(3)=sqrt(l32)
124 s(4)=sqrt(l42)
125 r(5)=r(1)
126
127 DO i=1,4
128 cs(i)=dksi(i)/s(i)
129 sn(i)=deta(i)/s(i)
130 ENDDO
131
132 rvlh=zero
133
134 IF (zls/=zero) THEN
135 DO i=1,4
136 j=i+1
137 rvlh=rvlh
138 . +atan((deta(i)*((xls-ksi(i))**2+zls**2)-dksi(i)*(xls-ksi(i))*(yls-eta(i)))/(r
139 . -atan((deta(i)*((xls-ksi(j))**2+zls**2)-dksi(i)*(xls-ksi(j))*(yls-eta(j)))/(r(j)*zls*dksi(i)))
140 ENDDO
141 ENDIF
142
143
144 rvlg=zero
145 DO i=1,4
146 j=i+1
147 v=(xls-ksi(i))*sn(i)-(yls-eta(i))*cs(i)
148 arg=(r(i)+r(j)-s(i))/(r(i)+r(j)+s(i))
149 IF (arg>zero) THEN
150 fln=-log(arg)
151 rvlg=rvlg+v*fln
152 ELSE
153 WRITE(6,'(A,E13.5)') '? ARG=',arg
154 ENDIF
155 ENDDO
156 rvlg=-rvlg+zls*rvlh
157
158 ENDIF
159
160 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)