46
47
48
49#include "implicit_f.inc"
50
51
52
53#include "mvsiz_p.inc"
54#include "com08_c.inc"
55#include "scr05_c.inc"
56
57
58
59 INTEGER NEL,IOUT,IPROP,NUVAR,IGTYP,ISENS
61 . uvar(nuvar,*), fx(*), fy(*), fz(*), eint(*),
62 . xmom(*), ymom(*), zmom(*), xkm(*), xkr(*),
63 . xcm(*) ,xcr(*), mass(*) ,iner(*),
64 . off(*), rot1(3,mvsiz), rot2(3,mvsiz),
65 . dx(*), dy(*), dz(*), rx(*), ry(*), rz(*), x0_err(3,*)
66 DOUBLE PRECISION XL(MVSIZ,3)
67
68
69
70 INTEGER I, K, JTYP, ISK1, ISK2,
71 . IFUN_XX,IFUN_YY,IFUN_ZZ,IFUN_RX,IFUN_RY,IFUN_RZ,
72 . IFUN_CXX,IFUN_CYY,IFUN_CZZ,IFUN_CRX,IFUN_CRY,IFUN_CRZ,
73 . IFUN_FMX,IFUN_FMY,IFUN_FMZ,IFUN_FMRX,IFUN_FMRY,IFUN_FMRZ,
74 . KFUNC,KMAT,KPROP,GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU,
75 . FDOF(6),ICOMBT,ICOMBR
77 . lx2,ly2,lz2,cxx,cyy,czz,crx,cry,crz,
78 . cr1,cr2,cr3,cr4,cr5,cr6,cx, ms,in,fac_ctx,fac_crx,
79 . my1(mvsiz),my2(mvsiz),mz1(mvsiz),mz2(mvsiz),mx1(mvsiz),
80 . mx2(mvsiz),fold(mvsiz,6),dxold(mvsiz,6),l2(mvsiz),
81 . knx(mvsiz),kny(mvsiz),knz(mvsiz),
82 . krx(mvsiz),kry(mvsiz),krz(mvsiz),
83 . vx(mvsiz),vy(mvsiz),vz(mvsiz),
84 . vx2(mvsiz),vy2(mvsiz),vz2(mvsiz),
85 . th1(mvsiz),th2(mvsiz),th3(mvsiz),
86 . vrx(mvsiz),vry(mvsiz),vrz(mvsiz),
87 . get_u_mat, get_u_geo, get_u_func,
88 . sma(6),smi(6),knn(mvsiz),krr(mvsiz),cr,
89 . dxs(mvsiz),dys(mvsiz),dzs(mvsiz),drxs(mvsiz),
90 . drys(mvsiz),drzs(mvsiz),kf(6),dxsk(6),fm(6),
91 . fcomb(6),deq(mvsiz),req(mvsiz),xcent(6),smeqt,smeqr,
92 . fac_loc_l,fac_loc_t,fac_x,fac_r
93 DOUBLE PRECISION X0DP(3)
94
95 EXTERNAL get_u_mat,get_u_geo, get_u_func,
97
98 parameter(kfunc=29)
99 parameter(kmat=31)
100 parameter(kprop=33)
101
102 dt = dt1
103 IF (dt==zero) dt = ep30
104 jtyp = nint(get_u_geo(1,iprop))
105
106 cr = get_u_geo(12,iprop)
107 cr1 = get_u_geo(15,iprop)
108 cr2 = get_u_geo(16,iprop)
109 cr3 = get_u_geo(17,iprop)
110 cr4 = get_u_geo(18,iprop)
111 cr5 = get_u_geo(19,iprop)
112 cr6 = get_u_geo(20,iprop)
113 cxx = get_u_geo(21,iprop)
114 cyy = get_u_geo(22,iprop)
115 czz = get_u_geo(23,iprop)
116 crx = get_u_geo(24,iprop)
117 cry = get_u_geo(25,iprop)
118 crz = get_u_geo(26,iprop)
119 fac_loc_l = get_u_geo(27,iprop)
120 fac_loc_t = get_u_geo(28,iprop)
121
122 fac_x = one / fac_loc_l
123 fac_r = one
124 fac_ctx = fac_loc_t / fac_loc_l
125 fac_crx = fac_loc_t
126
127 IF (igtyp==45) THEN
135 k = 28
136
137 xcent(1:6) = zero
138 DO i=1,6
139 kf(i) = get_u_geo(40+i,iprop)
140 fm(i) = get_u_geo(46+i,iprop)
141 smi(i) = get_u_geo(k+1,iprop)
142 sma(i) = get_u_geo(k+2,iprop)
143 fcomb(i) = get_u_geo(54+i,iprop)
144 xcent(i) = half*(sma(i)+smi(i))
145 k = k+2
146 END DO
147
148 icombt = nint(fcomb(1)+fcomb(2)+fcomb(3))
149 icombr = nint(fcomb(4)+fcomb(5)+fcomb(6))
150 ELSE
151
152 icombt = 0
153 icombr = 0
154 ENDIF
155
162
169
170 DO i=1,nel
171 knx(i) = uvar(19,i)
172 kny(i) = uvar(20,i)
173 knz(i) = uvar(21,i)
174 krx(i) = uvar(31,i)
175 kry(i) = uvar(32,i)
176 krz(i) = uvar(33,i)
177 knn(i) = uvar(17,i)
178 krr(i) = uvar(18,i)
179 IF ((igtyp==45).AND.(isens==1)) THEN
180 DO k=1,6
181 dxsk(k) = uvar(4+k-1,i)
182 IF (abs(sma(k))>em20) dxsk(k) =
min(dxsk(k),sma(k))
183 IF (abs(smi(k))>em20) dxsk(k) =
max(dxsk(k),smi(k))
184 END DO
185 dxs(i) = dxsk(1)
186 dys(i) = dxsk(2)
187 dzs(i) = dxsk(3)
188 drxs(i) = dxsk(4)
189 drys(i) = dxsk(5)
190 drzs(i) = dxsk(6)
191 ENDIF
192 dxold(i,1) = dx(i)
193 dxold(i,2) = dy(i)
194 dxold(i,3) = dz(i)
195 dxold(i,4) = rx(i)
196 dxold(i,5) = ry(i)
197 dxold(i,6) = rz(i)
198 fold(i,1) = fx(i)
199 fold(i,2) = fy(i)
200 fold(i,3) = fz(i)
201 fold(i,4) = xmom(i)
202 fold(i,5) = ymom(i)
203 fold(i,6) = zmom(i)
204 IF (iresp == 1) THEN
205
206 IF (tt == zero) THEN
207 uvar(1,i)=xl(i,1)
208 uvar(2,i)=xl(i,2)
209 uvar(3,i)=xl(i,3)
210 x0_err(1,i)= xl(i,1)-uvar(1,i)
211 x0_err(2,i)= xl(i,2)-uvar(2,i)
212 x0_err(3,i)= xl(i,3)-uvar(3,i)
213 ENDIF
214 x0dp(1) = uvar(1,i)
215 x0dp(2) = uvar(2,i)
216 x0dp(3) = uvar(3,i)
217 x0dp(1) = x0dp(1) + x0_err(1,i)
218 x0dp(2) = x0dp(2) + x0_err(2,i)
219 x0dp(3) = x0dp(3) + x0_err(3,i)
220 ELSE
221
222 x0dp(1) = uvar(1,i)
223 x0dp(2) = uvar(2,i)
224 x0dp(3) = uvar(3,i)
225 ENDIF
226
227 dx(i) = xl(i,1)-x0dp(1)
228 dy(i) = xl(i,2)-x0dp(2)
229 dz(i) = xl(i,3)-x0dp(3)
230 IF (icombt > 1) THEN
231 deq(i) = sqrt(fcomb(1)*(dx(i)-xcent(1))*(dx(i)-xcent(1))
232 . +fcomb(2)*(dy(i)-xcent(2))*(dy(i)-xcent(2))
233 . +fcomb(3)*(dz(i)-xcent(3))*(dz(i)-xcent(3)))
234 ELSE
235 deq(i) = zero
236 ENDIF
237
238 rx(i) = rot2(1,i)-rot1(1,i)
239 ry(i) = rot2(2,i)-rot1(2,i)
240 rz(i) = rot2(3,i)-rot1(3,i)
241 IF (icombr > 1) THEN
242 req(i) = sqrt(fcomb(4)*(rx(i)-xcent(4))*(rx(i)-xcent(4))
243 . +fcomb(5)*(ry(i)-xcent(5))*(ry(i)-xcent(5))
244 . +fcomb(6)*(rz(i)-xcent(6))*(rz(i)-xcent(6)))
245 ELSE
246 req(i) = zero
247 ENDIF
248
249 vx(i) = (dx(i) - dxold(i,1)) / dt
250 vy(i) = (dy(i) - dxold(i,2)) / dt
251 vz(i) = (dz(i) - dxold(i,3)) / dt
252 vrx(i) = (rx(i) - dxold(i,4)) / dt
253 vry(i) = (ry(i) - dxold(i,5)) / dt
254 vrz(i) = (rz(i) - dxold(i,6)) / dt
255 xkm(i) = zero
256 xkr(i) = zero
257 xcm(i) = zero
258 xcr(i) = zero
259 ENDDO
260
261 CALL xddl33(nel, dx, fx, knx, ifun_xx, xkm,fac_x)
262 CALL xddl33(nel, dy, fy, kny, ifun_yy, xkm,fac_x)
263 CALL xddl33(nel, dz, fz, knz, ifun_zz, xkm,fac_x)
264 CALL xddl33(nel, rx, xmom, krx, ifun_rx, xkr,fac_r)
265 CALL xddl33(nel, ry, ymom, kry, ifun_ry, xkr,fac_r)
266 CALL xddl33(nel, rz, zmom, krz, ifun_rz, xkr,fac_r)
267
268
269 IF (igtyp==45) THEN
270 CALL stdpl(nel,dx,vx,fx,knx,knn,kf(1),cr,cr1,smi(1),sma(1),
271 . xkm,fm(1),uvar(10,1),ifun_fmx,icombt,deq,fcomb(1),xcent(1),fac_x)
272 CALL stdpl(nel,dy,vy,fy,kny,knn,kf(2),cr,cr2,smi(2),sma(2),
273 . xkm,fm(2),uvar(11,1),ifun_fmy,icombt,deq,fcomb(2),xcent(2),fac_x)
274 CALL stdpl(nel,dz,vz,fz,knz,knn,kf(3),cr,cr3,smi(3),sma(3),
275 . xkm,fm(3),uvar(12,1),ifun_fmz,icombt,deq,fcomb(3),xcent(3),fac_x)
276
277 CALL stdpl(nel,rx,vrx,xmom,krx,krr,kf(4),cr,cr4,smi(4),sma(4),
278 . xkr,fm(4),uvar(13,1),ifun_fmrx,icombr,req,fcomb(4),xcent(4),fac_r)
279 CALL stdpl(nel,ry,vry,ymom,kry,krr,kf(5),cr,cr5,smi(5),sma(5),
280 . xkr,fm(5),uvar(14,1),ifun_fmry,icombr,req,fcomb(5),xcent(5),fac_r)
281 CALL stdpl(nel,rz,vrz,zmom,krz,krr,kf(6),cr,cr6,smi(6),sma(6),
282 . xkr,fm(6),uvar(15,1),ifun_fmrz,icombr,req,fcomb(6),xcent(6),fac_r)
283 ENDIF
284
285 IF ((igtyp==45).AND.(isens==1)) THEN
286 CALL sens_block(nel,dx,fx,knx,knn,cr,cr1,dxs,fdof(1),xkm)
287 CALL sens_block(nel,dy,fy,kny,knn,cr,cr2,dys,fdof(2),xkm)
288 CALL sens_block(nel,dz,fz,knz,knn,cr,cr3,dzs,fdof(3),xkm)
289 CALL sens_block(nel,rx,xmom,krx,krr,cr,cr4,drxs,fdof(4),xkr)
290 CALL sens_block(nel,ry,ymom,kry,krr,cr,cr5,drys,fdof(5),xkr)
291 CALL sens_block(nel,rz,zmom,krz,krr,cr,cr6,drzs,fdof(6),xkr)
292 ENDIF
293
294 DO i=1,nel
295 ms = zero
296 in = zero
297 IF((uvar(34,i)+uvar(35,i))/=zero)
298 . ms = (uvar(34,i)*uvar(35,i))/(uvar(34,i)+uvar(35,i))
299 IF((uvar(36,i)+uvar(37,i))/=zero)
300 . in = (uvar(36,i)*uvar(37,i))/(uvar(36,i)+uvar(37,i))
301 cx =
max(cr1,cr2,cr3,cr4,cr5,cr6)
302 xcm(i) = cx*sqrt(xkm(i)*ms)
303 xcr(i) = cx*sqrt(xkr(i)*in)
304
305 fx(i)= fx(i) + cr1*sqrt(knx(i)*ms)*vx(i)
306 fy(i)= fy(i) + cr2*sqrt(kny(i)*ms)*vy(i)
307 fz(i)= fz(i) + cr3*sqrt(knz(i)*ms)*vz(i)
308 xmom(i)= xmom(i) + cr4*sqrt(krx(i)*in)*vrx(i)
309 ymom(i)= ymom(i) + cr5*sqrt(kry(i)*in)*vry(i)
310 zmom(i)= zmom(i) + cr6*sqrt(krz(i)*in)*vrz(i)
311 ENDDO
312
313 CALL xddl33i(nel, vx, fx, cxx, ifun_cxx, xcm, fac_ctx)
314 CALL xddl33i(nel, vy, fy, cyy, ifun_cyy, xcm, fac_ctx)
315 CALL xddl33i(nel, vz, fz, czz, ifun_czz, xcm, fac_ctx)
316 CALL xddl33i(nel, vrx, xmom, crx, ifun_crx, xcr, fac_crx)
317 CALL xddl33i(nel, vry, ymom, cry, ifun_cry, xcr, fac_crx)
318 CALL xddl33i(nel, vrz, zmom, crz, ifun_crz, xcr, fac_crx)
319
320
321 DO i=1,nel
322 eint(i) = eint(i) + half*(
323 . (dx(i)-dxold(i,1)) * (fx(i)+fold(i,1))
324 . + (dy(i)-dxold(i,2)) * (fy(i)+fold(i,2))
325 . + (dz(i)-dxold(i,3)) * (fz(i)+fold(i,3))
326 . + (rx(i)-dxold(i,4)) * (xmom(i)+fold(i,4))
327 . + (ry(i)-dxold(i,5)) * (ymom(i)+fold(i,5))
328 . + (rz(i)-dxold(i,6)) * (zmom(i)+fold(i,6)))
329 mass(i) = zero
330 iner(i) = zero
331 ENDDO
332
333 RETURN
subroutine xddl33i(nel, dx, fx, kx, ifun, kmx, fac_x)
subroutine xddl33(nel, dx, fx, kx, ifun, kmx, fac_x)
subroutine sens_block(nel, dx, fx, kx, knn, cr, crx, dxs, flag, kmx)
subroutine stdpl(nel, dx, vx, fx, kx, knn, kf, cr, crx, dxmi, dxma, kmx, fm, frp, ifun, icomb, deq, fcomb, xcent, fac_x)
subroutine def_fdof(jtyp, fdof)
integer function get_u_pid(ip)
integer function get_u_pnu(ivar, ip, k)
integer function get_u_mid(im)
integer function get_u_mnu(ivar, im, k)