32
33
34
35#include "implicit_f.inc"
36
37
38
39#include "mvsiz_p.inc"
40
41
42
43#include "vect01_c.inc"
44#include "com04_c.inc"
45#include "com08_c.inc"
46#include "param_c.inc"
47
48
49
50 INTEGER :: IMS(NUMNOD,*), NC1(*), NC2(*), NC3(*), NC4(*)
51 my_real :: v(3,*), w(3,*), fill(numnod,*), dfill(numnod,*), x(3,*), dalph1(*), dalph2(*)
52
53
54
55 INTEGER I, N1, N2, N3, N4, NP
56 my_real fi1(mvsiz), fi2(mvsiz), fi3(mvsiz), fi4(mvsiz),
57 . fa(mvsiz), vdy1(mvsiz), vdy2(mvsiz), vdy3(mvsiz), vdy4(mvsiz
58 . vdz1(mvsiz), vdz2(mvsiz), vdz3(mvsiz), vdz4(mvsiz), vdy(mvsiz), vdz(mvsiz
59 . abf, dn, p1, p2, p3, p4, pt, psy, psz, pty,
60 . ptz, ps, pst, pts, ds0, dt0, ds, dt,
61 . df1(mvsiz), df2(mvsiz), df3(mvsiz), df4(mvsiz)
62
63 DO i=lft,llt
64 fi1(i)=fill(nc1(i),1)
65 fi2(i)=fill(nc2(i),1)
66 fi3(i)=fill(nc3(i),1)
67 fi4(i)=fill(nc4(i),1)
68 abf=abs(fi1(i))+abs(fi2(i))+abs(fi3(i))+abs(fi4(i))
69 n1=nint(sign(one,fi1(i)))
70 n2=nint(sign(one,fi2(i)))
71 n3=nint(sign(one,fi3(i)))
72 n4=nint(sign(one,fi4(i)))
74 dn=dt1*np
75 IF(dn/=zero)THEN
76 fa(i)=-dalph1(i)*abf/dn
77 ELSE
78 fa(i)=zero
79 ENDIF
80 ENDDO
81
82
83
84 DO i=lft,llt
85 vdy1(i)=v(2,nc1(i)) - w(2,nc1(i))
86 vdz1(i)=v(3,nc1(i)) - w(3,nc1(i))
87
88 vdy2(i)=v(2,nc2(i)) - w(2,nc2(i))
89 vdz2(i)=v(3,nc2(i)) - w(3,nc2(i))
90
91 vdy3(i)=v(2,nc3(i)) - w(2,nc3(i))
92 vdz3(i)=v(3,nc3(i)) - w(3,nc3(i))
93
94 vdy4(i)=v(2,nc4(i)) - w(2,nc4(i))
95 vdz4(i)=v(3,nc4(i)) - w(3,nc4(i))
96 ENDDO
97
98
99
100 DO i=lft,llt
101 p1=fi1(i) + one
102 p2=fi2(i) + one
103 p3=fi3(i) + one
104 p4=fi4(i) + one
105 pt=(p1+p2+p3+p4)
107 vdy(i)=(vdy1(i)*p1+vdy2(i)*p2+vdy3(i)*p3+vdy4(i)*p4)/pt
108 vdz(i)=(vdz1(i)*p1+vdz2(i)*p2+vdz3(i)*p3+vdz4(i)*p4)/pt
109 ENDDO
110
111 DO i=lft,llt
112 psy=-x(2,nc1(i))+x(2,nc2(i))+x(2,nc3(i))-x(2,nc4(i))
113 psz=-x(3,nc1(i))+x(3,nc2(i))+x(3,nc3(i))-x(3,nc4(i))
114 pty=-x(2,nc1(i))-x(2,nc2(i))+x(2,nc3(i))+x(2,nc4(i))
115 ptz=-x(3,nc1(i))-x(3,nc2(i))+x(3,nc3(i))+x(3,nc4(i))
116 ps=sqrt(psy**2+psz**2)
117 pt=sqrt(pty**2+ptz**2)
118 pst=psy*ptz-psz*pty
119 pts=-pst
120 ds0=-four*(pty*vdz(i)-ptz*vdy(i))/pts
121 dt0=-four*(psy*vdz(i)-psz*vdy(i))/pst
122 IF(fi1(i)>=zero)THEN
123 ds=-four*(pty*vdz1(i)-ptz*vdy1(i))/pts
124 dt=-four*(psy*vdz1(i)-psz*vdy1(i))/pst
125 ELSE
126 ds=ds0
127 dt=dt0
128 ENDIF
131 df1(i)=fourth*((-two*ds-two*dt+ds*dt*dt1)*fi1(i)
132 . + ( two*ds-ds*dt*dt1)*fi2(i)
133 . + ( ds*dt*dt1)*fi3(i)
134 . + ( two*dt-ds*dt*dt1)*fi4(i) )
135 IF(fi2(i)>=zero)THEN
136 ds=-four*(pty*vdz2(i)-ptz*vdy2(i))/pts
137 dt=-four*(psy*vdz2(i)-psz*vdy2(i))/pst
138 ELSE
139 ds=ds0
140 dt=dt0
141 ENDIF
144 df2(i)=fourth*(( -two*ds+ds*dt*dt1)*fi1(i)
145 . + ( two*ds-two*dt-ds*dt*dt1)*fi2(i)
146 . + ( +two*dt+ds*dt*dt1)*fi3(i)
147 . + ( -ds*dt*dt1)*fi4(i) )
148 IF(fi3(i)>=zero)THEN
149 ds=-four*(pty*vdz3(i)-ptz*vdy3(i))/pts
150 dt=-four*(psy*vdz3(i)-psz*vdy3(i))/pst
151 ELSE
152 ds=ds0
153 dt=dt0
154 ENDIF
157 df3(i)=fourth*(( +ds*dt*dt1)*fi1(i)
158 . + ( -two*dt-ds*dt*dt1)*fi2(i)
159 . + (+two*ds+two*dt+ds*dt*dt1)*fi3(i)
160 . + (-two*ds -ds*dt*dt1)*fi4(i) )
161 IF(fi4(i)>=zero)THEN
162 ds=-four*(pty*vdz4(i)-ptz*vdy4(i))/pts
163 dt=-four*(psy*vdz4(i)-psz*vdy4(i))/pst
164 ELSE
165 ds=ds0
166 dt=dt0
167 ENDIF
170 df4(i)=fourth*(( -two*dt+ds*dt*dt1)*fi1(i)
171 . + ( -ds*dt*dt1)*fi2(i)
172 . + (+two*ds +ds*dt*dt1)*fi3(i)
173 . + (-two*ds+two*dt-ds*dt*dt1)*fi4(i) )
174 ENDDO
175
177
178 DO i = lft,llt
179 dfill(nc1(i),1)=dfill(nc1(i),1)+df1(i)-fa(i)
180 dfill(nc2(i),1)=dfill(nc2(i),1)+df2(i)-fa(i)
181 dfill(nc3(i),1)=dfill(nc3(i),1)+df3(i)-fa(i)
182 dfill(nc4(i),1)=dfill(nc4(i),1)+df4(i)-fa(i)
183 ims(nc1(i),1)=ims(nc1(i),1)+1
184 ims(nc2(i),1)=ims(nc2(i),1)+1
185 ims(nc3(i),1)=ims(nc3(i),1)+1
186 ims(nc4(i),1)=ims(nc4(i),1)+1
187 ENDDO
188
190
191 IF(jmult>1)THEN
192
193 DO i=lft,llt
194 fi1(i)=fill(nc1(i),2)
195 fi2(i)=fill(nc2(i),2)
196 fi3(i)=fill(nc3(i),2)
197 fi4(i)=fill(nc4(i),2)
198 abf=abs(fi1(i))+abs(fi2(i))+abs(fi3(i))+abs(fi4(i))
199 n1=nint(sign(one,fi1(i)))
200 n2=nint(sign(one,fi2(i)))
201 n3=nint(sign(one,fi3(i)))
202 n4=nint(sign(one,fi4(i)))
204 dn=dt1*np
205 IF(dn/=zero)THEN
206 fa(i)=-dalph2(i)*abf/dn
207 ELSE
208 fa(i)=zero
209 ENDIF
210 ENDDO
211
212
213
214 DO i=lft,llt
215 p1=fi1(i) + one
216 p2=fi2(i) + one
217 p3=fi3(i) + one
218 p4=fi4(i) + one
219 pt=(p1+p2+p3+p4)
221 vdy(i)=(vdy1(i)*p1+vdy2(i)*p2+vdy3(i)*p3+vdy4(i)*p4)/pt
222 vdz(i)=(vdz1(i)*p1+vdz2(i)*p2+vdz3
223 ENDDO
224
225 DO i=lft,llt
226 psy=-x(2,nc1(i))+x(2,nc2(i))+x(2,nc3(i))-x(2,nc4(i))
227 psz=-x(3,nc1
228 pty=-x(2,nc1
229 ptz=-x(3,nc1(i))-x(3,nc2(i))+x(3,nc3(i))+x
230 ps=sqrt(psy**2+psz**2)
231 pt=sqrt(pty**2+ptz**2)
232 pst=psy*ptz-psz*pty
233 pts=-pst
234 ds0=-four*(pty*vdz(i)-ptz*vdy(i))/pts
235 dt0=-four*(psy*vdz(i)-psz*vdy(i))/pst
236 IF(fi1(i)>=zero)THEN
237 ds=-four*(pty*vdz1(i)-ptz*vdy1(i))/pts
238 dt=-four*(psy*vdz1(i)-psz*vdy1(i))/pst
239 ELSE
240 ds=ds0
241 dt=dt0
242 ENDIF
245 df1(i)=fourth*((-two*ds-two*dt+ds*dt*dt1)*fi1(i)
246 . +( two*ds -ds*dt*dt1)*fi2(i)
247 . +( ds*dt*dt1)*fi3(i)
248 . +( two*dt-ds*dt*dt1)*fi4(i) )
249 IF(fi2(i)>=zero)THEN
250 ds=-four*(pty*vdz2(i)-ptz*vdy2(i))/pts
251 dt=-four*(psy*vdz2(i)-psz*vdy2(i))/pst
252 ELSE
253 ds=ds0
254 dt=dt0
255 ENDIF
258 df2(i)=fourth*((-two*ds +ds*dt*dt1)*fi1(i)
259 . +( two*ds-two*dt-ds*dt*dt1)*fi2(i)
260 . +( +two*dt+ds*dt*dt1)*fi3(i)
261 . +( -ds*dt*dt1)*fi4(i) )
262 IF(fi3(i)>=zero)THEN
263 ds=-four*(pty*vdz3(i)-ptz*vdy3(i))/pts
264 dt=-four*(psy*vdz3(i)-psz*vdy3(i))/pst
265 ELSE
266 ds=ds0
267 dt=dt0
268 ENDIF
271 df3(i)=fourth*(( +ds*dt*dt1)*fi1(i)
272 . +( -two*dt-ds*dt*dt1)*fi2(i)
273 . +(+two*ds+two*dt+ds*dt*dt1)*fi3(i)
274 . +(-two*ds -ds*dt*dt1)*fi4(i) )
275 IF(fi4(i)>=zero)THEN
276 ds=-four*(pty*vdz4(i)-ptz*vdy4(i))/pts
277 dt=-four*(psy*vdz4(i)-psz*vdy4(i))/pst
278 ELSE
279 ds=ds0
280 dt=dt0
281 ENDIF
284 df4(i)=fourth*(( -two*dt+ds*dt*dt1)*fi1(i)
285 . +( -ds*dt*dt1)*fi2(i)
286 . +(+two*ds +ds*dt*dt1)*fi3(i)
287 . +(-two*ds+two*dt-ds*dt*dt1)*fi4(i) )
288 ENDDO
289
291
292 DO i = lft,llt
293 dfill(nc1(i),2)=dfill(nc1(i),2)+df1(i)-fa(i)
294 dfill(nc2(i),2)=dfill(nc2(i),2)+df2(i)-fa(i)
295 dfill(nc3(i),2)=dfill(nc3(i),2)+df3(i)-fa(i)
296 dfill(nc4(i),2)=dfill(nc4(i),2)+df4(i)-fa(i)
297 ims(nc1(i),2)=ims(nc1(i),2)+1
298 ims(nc2(i),2)=ims(nc2(i),2)+1
299 ims(nc3(i),2)=ims(nc3(i),2)+1
300 ims(nc4(i),2)=ims(nc4(i),2)+1
301 ENDDO
302
304
305 ENDIF
306
307 RETURN