OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fixvel.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| fixvel ../engine/source/constraints/general/impvel/fixvel.F
25!||--- called by ------------------------------------------------------
26!|| resol ../engine/source/engine/resol.F
27!||--- calls -----------------------------------------------------
28!|| vinter_smooth ../engine/source/tools/curve/vinter_smooth.F
29!|| vinterdp ../engine/source/tools/curve/vinterdp.F
30!||--- uses -----------------------------------------------------
31!|| nodal_arrays_mod ../common_source/modules/nodal_arrays.F90
32!|| python_call_funct_cload_mod ../engine/source/loads/general/python_call_funct_cload.F90
33!|| python_funct_mod ../common_source/modules/python_mod.F90
34!|| sensor_mod ../common_source/modules/sensor_mod.F90
35!||====================================================================
36 SUBROUTINE fixvel(IBFV ,A ,V ,NPC ,TF ,
37 2 VEL ,MS ,X ,SKEW ,AR ,
38 3 VR ,IN ,NSENSOR,SENSOR_TAB,
39 4 WEIGHT,DEPLA ,RBY ,IFRAME,
40 5 XFRAME,DR ,NODNX_SMS, NODES,
41 6 TT_DOUBLE,DEPLA_DOUBLE, PYTHON, WFEXT)
42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE python_funct_mod
46 USE python_call_funct_cload_mod, ONLY: python_call_funct_cload
47 USE nodal_arrays_mod
48 USE sensor_mod
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C G l o b a l P a r a m e t e r s
55C-----------------------------------------------
56#include "mvsiz_p.inc"
57C-----------------------------------------------
58C C o m m o n B l o c k s
59C-----------------------------------------------
60#include "com01_c.inc"
61#include "com04_c.inc"
62#include "com08_c.inc"
63#include "param_c.inc"
64#include "parit_c.inc"
65#include "sms_c.inc"
66C-----------------------------------------------
67C D u m m y A r g u m e n t s
68C-----------------------------------------------
69 TYPE(nodal_arrays_) ,INTENT(IN) :: NODES
70 INTEGER ,INTENT(IN) :: NSENSOR
71 INTEGER NPC(*)
72 INTEGER IBFV(NIFV,*),WEIGHT(*),IFRAME(LISKN,*),NODNX_SMS(*)
73 my_real, DIMENSION(3,*), TARGET :: DEPLA
74 REAL(KIND=8), dimension(3,*), TARGET :: depla_double
75 REAL(KIND=8) :: tt_double
76 my_real
77 . a(3,*), v(3,*), tf(*), vel(lfxvelr,*), ms(*), x(3,*),
78 . skew(lskew,*), ar(3,*), vr(3,*), in(*),
79 . rby(nrby,*),xframe(nxframe,*), dr(3,*)
80 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
81 TYPE(python_), INTENT(INOUT) :: PYTHON
82 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
83C-----------------------------------------------
84C L o c a l V a r i a b l e s
85C-----------------------------------------------
86 INTEGER N, I, ISK, J, L, K1, K2, K3, ISENS,K,
87 . ii, ic, nn, ideb, nr, nsk, nfk, ifm, n0,
88 . ilenc(mvsiz), iposc(mvsiz), iadc(mvsiz),
89 . lc(mvsiz), index(mvsiz), icoor, ismooth,jj
90 my_real ::
91 . ax, axi, aold, vv, a0, aa, startt, stopt,
92 . dydx,dw,dw2,dd,rx,ry,rz,vf,vfx,vfy,vfz,afx,afy,afz,rw_sms(mvsiz),
93 . r, cth, sth, skdir(3), fmdir(3), hx, hy, hz, coef,
94 . er(3),et(3),alpha, nor1, nor2, nor3, yimp, skewglob,
95 . aold0(6,mvsiz),dfj,qf0(9),ej(3),ddf,ddr,rf
96
97 REAL(KIND=8) :: tstmp,fac,facx,yc0(mvsiz),yc(mvsiz),tsc0(mvsiz),tsc(mvsiz),dydxc(mvsiz)
98 INTEGER, DIMENSION(MVSIZ) :: ILENC_SMOOTH,IPOSC_SMOOTH,IADC_SMOOTH
99 my_real, DIMENSION(MVSIZ) :: YC0_SMOOTH,YC_SMOOTH,TSC0_SMOOTH,TSC_SMOOTH,DYDXC_SMOOTH
100 INTEGER :: NINDEX_SMOOTH,NINDEX_SMOOTH_WSENS
101 INTEGER, DIMENSION(MVSIZ) ::INDEX_SMOOTH
102
103 INTEGER, DIMENSION(MVSIZ) :: ILENC_WO_SMOOTH,IPOSC_WO_SMOOTH,IADC_WO_SMOOTH
104 REAL(KIND=8), dimension(mvsiz) :: yc0_wo_smooth,yc_wo_smooth,tsc0_wo_smooth,tsc_wo_smooth,dydxc_wo_smooth
105 INTEGER :: NINDEX_WO_SMOOTH,NINDEX_WO_SMOOTH_WSENS
106 INTEGER, DIMENSION(MVSIZ) ::INDEX_WO_SMOOTH
107
108 REAL(KIND=8), dimension(:,:), POINTER :: d
109 REAL(KIND=8) :: dt2_double
110
111 INTEGER, DIMENSION(MVSIZ) :: NSENS
112 REAL(KIND=8), DIMENSION(MVSIZ) :: TS
113C=======================================================================
114#if MYREAL4
115! Simple precision
116 d => depla_double(1:3,1:numnod)
117#elif 1
118! Double precision
119 d => depla(1:3,1:numnod)
120#endif
121 dt2_double = dt2
122 IF(n2d == 1)THEN
123 ax=one
124 ELSE
125 ax=zero
126 ENDIF
127C
128 nsens(1:mvsiz) = 0
129 ts(1:mvsiz) = zero
130C
131 ideb = 0
132 ismooth = 0
133 DO nn=1,nfxvel,nvsiz
134 IF (ibfv(8,nn) == 1) cycle
135 ic = 0
136 IF (nsensor > 0) THEN
137 DO ii = 1, min(nfxvel-ideb,nvsiz)
138 n = ii+ideb
139 IF (ibfv(13,n) > 0 ) cycle
140 startt = vel(2,n)
141 stopt = vel(3,n)
142 facx = vel(5,n)
143 IF (tt_double < startt) cycle
144 IF (tt_double > stopt) cycle
145 i = iabs(ibfv(1,n))
146 isens = ibfv(4,n)
147 IF (isens > 0) THEN
148 tstmp = tt_double - sensor_tab(isens)%TSTART
149 IF (tstmp < zero) cycle
150 ELSE
151 tstmp = tt_double
152 ENDIF
153 ic = ic + 1
154 index(ic) = n
155 nsens(ic) = isens
156 ts(ic) = tstmp
157 IF (ibfv(7,n) == 1) THEN
158 tsc(ic) = (ts(ic)+half*dt2_double)*facx
159 ELSEIF (ibfv(7,n) == 2) THEN
160 tsc0(ic) = ts(ic)*facx
161 tsc(ic) = (ts(ic)+dt2_double)*facx
162 ELSE
163 tsc(ic) = (ts(ic)+dt2_double)*facx
164 ENDIF
165 IF(idtmins==0.AND.idtmins_int==0)THEN
166 rw_sms(ic)=one
167 ELSE
168 isk = ibfv(2,n)/10
169 ifm = ibfv(9,n)
170 icoor = ibfv(10,n)
171 j = ibfv(2,n)
172 IF (ifm<=1) j=j-10*isk
173 IF(nodnx_sms(i)==0)THEN
174 rw_sms(ic)=one
175 ELSE
176 rw_sms(ic)=zero
177 END IF
178 END IF
179 ENDDO
180 ELSE
181C sans sensor
182 ts = tt_double
183 DO ii = 1, min(nfxvel-ideb,nvsiz)
184 n = ii+ideb
185 IF(ibfv(13,n) > 0 ) cycle
186 startt = vel(2,n)
187 stopt = vel(3,n)
188 facx = vel(5,n)
189 IF (tt_double < startt) cycle
190 IF (tt_double > stopt) cycle
191 i=iabs(ibfv(1,n))
192 ic = ic + 1
193 index(ic) = n
194 IF(ibfv(7,n) == 1) THEN
195 tsc(ic) = (tt_double + half*dt2_double)*facx
196 ELSE
197 tsc(ic) = (tt_double+dt2_double)*facx
198 ENDIF
199 IF(idtmins==0.AND.idtmins_int==0)THEN
200 rw_sms(ic)=one
201 ELSE
202 isk=ibfv(2,n)/10
203 ifm = ibfv(9,n)
204 icoor = ibfv(10,n)
205 j=ibfv(2,n)
206 IF (ifm<=1) j=j-10*isk
207 IF(nodnx_sms(i)==0)THEN
208 rw_sms(ic)=one
209 ELSE
210 rw_sms(ic)=zero
211 END IF
212 END IF
213 ENDDO
214 ENDIF
215C
216 ideb = ideb + min(nfxvel-ideb,nvsiz)
217C-------
218 nindex_smooth = 0
219 nindex_wo_smooth = 0
220C------- First pass for impdisp with sensor - in the top of the list
221 DO ii=1,ic
222 n = index(ii)
223 l = ibfv(3,n)
224 ismooth = 0
225 IF (l > 0) ismooth = npc(2*nfunct+l+1)
226 IF(nsens(ii) > 0 .and. ts(ii) >= zero .and. ibfv(7,n) == 2) THEN
227 IF(ismooth>0) THEN
228 nindex_smooth = nindex_smooth + 1
229 index_smooth(nindex_smooth) = ii
230 tsc_smooth(nindex_smooth) = tsc(ii)
231 tsc0_smooth(nindex_smooth) = tsc0(ii)
232 ELSE IF(ismooth==0) THEN
233 nindex_wo_smooth = nindex_wo_smooth + 1
234 index_wo_smooth(nindex_wo_smooth) = ii
235 tsc_wo_smooth(nindex_wo_smooth) = tsc(ii)
236 tsc0_wo_smooth(nindex_wo_smooth) = tsc0(ii)
237 ELSE IF(ismooth<0) THEN
238 l = -ismooth
239 i = iabs(ibfv(1,n))
240 CALL python_call_funct_cload(python, l, tsc(ii),yc(ii),i,nodes)
241 CALL python_call_funct_cload(python, l, tsc0(ii),yc0(ii),i,nodes)
242 ENDIF
243 ENDIF
244 ENDDO
245 nindex_smooth_wsens = nindex_smooth
246 nindex_wo_smooth_wsens = nindex_wo_smooth
247C------- Second pass
248 DO ii=1,ic
249 n = index(ii)
250 l = ibfv(3,n)
251 ismooth = 0
252 IF (l > 0) ismooth = npc(2*nfunct+l+1)
253 IF(nsens(ii) == 0 .or. ts(ii) < zero .or. ibfv(7,n) /= 2) THEN
254 IF(ismooth>0) THEN
255 nindex_smooth = nindex_smooth + 1
256 index_smooth(nindex_smooth) = ii
257 tsc_smooth(nindex_smooth) = tsc(ii)
258 ELSE IF(ismooth==0) THEN
259 nindex_wo_smooth = nindex_wo_smooth + 1
260 index_wo_smooth(nindex_wo_smooth) = ii
261 tsc_wo_smooth(nindex_wo_smooth) = tsc(ii)
262 ELSE IF(ismooth<0) THEN
263 l = -ismooth
264 i = iabs(ibfv(1,n))
265 CALL python_call_funct_cload(python, l, tsc(ii),yc(ii),i,nodes)
266 ENDIF
267 ENDIF
268 ENDDO
269
270 DO ii=1,nindex_wo_smooth
271 jj = index_wo_smooth(ii)
272 n = index(jj)
273 l = ibfv(3,n)
274 IF (l > 0) ismooth = npc(2*nfunct+l+1)
275 IF(ncycle == 0)THEN
276 iposc_wo_smooth(ii) = 0
277 ELSE
278 iposc_wo_smooth(ii) = ibfv(5,n)
279 ENDIF
280 iadc_wo_smooth(ii) = npc(l) / 2 + 1
281 ilenc_wo_smooth(ii) = npc(l+1) / 2 - iadc_wo_smooth(ii) - iposc_wo_smooth(ii)
282 ENDDO
283
284 DO ii=1,nindex_smooth
285 jj = index_smooth(ii)
286 n = index(jj)
287 l = ibfv(3,n)
288 IF (l > 0) ismooth = npc(2*nfunct+l+1)
289 IF(ncycle == 0)THEN
290 iposc_smooth(ii) = 0
291 ELSE
292 iposc_smooth(ii) = ibfv(5,n)
293 ENDIF
294 iadc_smooth(ii) = npc(l) / 2 + 1
295 ilenc_smooth(ii) = npc(l+1) / 2 - iadc_smooth(ii) - iposc_smooth(ii)
296 ENDDO
297
298 IF (nindex_wo_smooth > 0) THEN
299 IF (nindex_wo_smooth_wsens > 0)
300 . CALL vinterdp(tf,iadc_wo_smooth,iposc_wo_smooth,ilenc_wo_smooth,nindex_wo_smooth_wsens,
301 1 tsc0_wo_smooth,dydxc_wo_smooth,yc0_wo_smooth)
302 CALL vinterdp(tf,iadc_wo_smooth,iposc_wo_smooth,ilenc_wo_smooth,nindex_wo_smooth,
303 1 tsc_wo_smooth,dydxc_wo_smooth,yc_wo_smooth)
304
305 DO ii=1,nindex_wo_smooth
306 jj = index_wo_smooth(ii)
307 yc(jj) = yc_wo_smooth(ii)
308 yc0(jj) = yc0_wo_smooth(ii)
309 iposc(jj) = iposc_wo_smooth(ii)
310 ENDDO
311 ENDIF
312 IF (nindex_smooth > 0) THEN
313 IF (nindex_smooth_wsens > 0)
314 . CALL vinter_smooth(tf,iadc_smooth,iposc_smooth,ilenc_smooth,nindex_smooth_wsens,
315 1 tsc0_smooth,dydxc_smooth,yc0_smooth)
316 CALL vinter_smooth(tf,iadc_smooth,iposc_smooth,ilenc_smooth,nindex_smooth,
317 1 tsc_smooth,dydxc_smooth,yc_smooth)
318 DO ii=1,nindex_smooth
319 jj = index_smooth(ii)
320 yc(jj) = yc_smooth(ii)
321 yc0(jj) = yc0_smooth(ii)
322 iposc(jj) = iposc_smooth(ii)
323 ENDDO
324 ENDIF
325
326
327
328 IF (ivector == 0) THEN
329 DO ii=1,ic
330 n = index(ii)
331 ibfv(5,n) = iposc(ii)
332 fac = vel(1,n)
333 yimp = vel(6,n)
334 vel(6,n) = yc(ii)
335C IF sms on the degree of freedom,
336C dw was already counted in WFEXT, in sms_fixvel...
337 dw = vel(4,n)
338 wfext = wfext + rw_sms(ii)*dw
339 vel(4,n) = (one-rw_sms(ii))*vel(4,n)
340 IF (nsens(ii) > 0 .and. ts(ii) > zero .and. ibfv(7,n) == 2) yc0(ii) = yc0(ii) * fac
341 yc(ii) = yc(ii) * fac
342 yimp = yimp * fac
343 i = iabs(ibfv(1,n))
344 isk = ibfv(2,n)/10
345 skewglob = 0
346 IF(isk/=0)THEN
347 IF (skew(1,isk) == one .AND. skew(2,isk) == zero .AND.
348 . skew(3,isk) == zero .AND. skew(4,isk) == zero .AND.
349 . skew(5,isk) == one .AND. skew(6,isk) == zero .AND.
350 . skew(7,isk) == zero .AND. skew(8,isk) == zero .AND.
351 . skew(9,isk) == one ) skewglob = 1
352 ENDIF
353 ifm = ibfv(9,n)
354 j = ibfv(2,n)
355 icoor = ibfv(10,n)
356 skdir = zero
357 fmdir = zero
358 r = one
359 IF (ifm<=1) j=j-10*isk
360 IF (j<=3) THEN
361 axi=one-ax+ax*x(2,i)
362 IF (isk<=1.AND.ifm<=1.AND.icoor == 0) THEN
363 IF (ibfv(7,n) == 2) THEN
364 IF (nsens(ii) > 0 .and. ts(ii) > zero) THEN
365 yc(ii) = (yc(ii)-yc0(ii))/dt2_double
366 ELSE
367 yc(ii)=(yc(ii)-d(j,i))/dt2_double
368 ENDIF
369 ENDIF
370 IF(ibfv(7,n)>=1) yc(ii)=(yc(ii)-v(j,i))/dt12
371 aold=a(j,i)
372 a(j,i)=yc(ii)
373 dw = fourth*ms(i)*axi*weight(i)
374 . * (a(j,i)*dt12 + two*v(j,i))*(a(j,i)-aold)
375
376 ELSEIF(isk > 1.AND.ifm<=1.AND.icoor == 0.AND.
377 . skewglob == 1)THEN
378 IF (ibfv(7,n) == 2) THEN
379 IF (nsens(ii) > 0 .and. ts(ii) > zero) THEN
380 yc(ii) = (yc(ii)-yc0(ii))/dt2_double
381 ELSE
382 yc(ii)=(yc(ii)-d(j,i))/dt2_double
383 ENDIF
384 ENDIF
385 IF (ibfv(7,n)>=1) yc(ii)=(yc(ii)-v(j,i))/dt12
386 aold=a(j,i)
387 a(j,i)=yc(ii)
388 dw = fourth*ms(i)*axi*weight(i)
389 . * (a(j,i)*dt12 + two*v(j,i))*(a(j,i)-aold)
390 ELSEIF (isk > 1 .AND. icoor == 0) THEN
391 k1=3*j-2
392 k2=3*j-1
393 k3=3*j
394 IF (ibfv(7,n) == 2) THEN
395 dd = skew(k1,isk)*d(1,i) +
396 . skew(k2,isk)*d(2,i) +
397 . skew(k3,isk)*d(3,i)
398 ENDIF
399 vv = skew(k1,isk)*v(1,i) +
400 . skew(k2,isk)*v(2,i) +
401 . skew(k3,isk)*v(3,i)
402 a0 = skew(k1,isk)*a(1,i) +
403 . skew(k2,isk)*a(2,i) +
404 . skew(k3,isk)*a(3,i)
405
406 IF (ibfv(7,n) == 2) THEN
407 IF (nsens(ii) > 0 .and. ts(ii) > zero) THEN
408 yc(ii) = (yc(ii)-yc0(ii))/dt2_double
409 ELSE
410 yc(ii) = (yc(ii)-dd)/dt2_double
411 ENDIF
412 ENDIF
413 IF (ibfv(7,n) >= 1) yc(ii) = (yc(ii)-vv)/dt12
414 aa = yc(ii) - a0
415 a(1,i)=a(1,i)+skew(k1,isk)*aa
416 a(2,i)=a(2,i)+skew(k2,isk)*aa
417 a(3,i)=a(3,i)+skew(k3,isk)*aa
418 dw = fourth*ms(i)*(yc(ii)*dt12+two*vv)*aa*axi*weight(i)
419 ELSEIF ((isk<=1.AND.ifm<=1.AND.icoor == 1) .OR.
420 & (isk > 1 .AND. icoor == 1)) THEN
421 coef = (skew(10,isk)-x(1,i))*skew(7,isk) +
422 & (skew(11,isk)-x(2,i))*skew(8,isk) +
423 & (skew(12,isk)-x(3,i))*skew(9,isk)
424 hx = coef*skew(7,isk)+x(1,i)-skew(10,isk)
425 hy = coef*skew(8,isk)+x(2,i)-skew(11,isk)
426 hz = coef*skew(9,isk)+x(3,i)-skew(12,isk)
427 r = sqrt(hx*hx+hy*hy+hz*hz)
428 IF(r<=em20) THEN
429 cth = zero
430 sth = zero
431 ELSE
432 cth = (skew(1,isk)*hx+skew(2,isk)*hy+skew(3,isk)*hz)/r
433 sth = (skew(4,isk)*hx+skew(5,isk)*hy+skew(6,isk)*hz)/r
434 ENDIF
435 IF (j == 1) THEN
436 skdir(1) =cth*skew(1,isk)+ sth*skew(4,isk)
437 skdir(2) =cth*skew(2,isk)+ sth*skew(5,isk)
438 skdir(3) =cth*skew(3,isk)+ sth*skew(6,isk)
439 ELSEIF(j == 2) THEN
440 er(1) =cth*skew(1,isk)+ sth*skew(4,isk)
441 er(2) =cth*skew(2,isk)+ sth*skew(5,isk)
442 er(3) =cth*skew(3,isk)+ sth*skew(6,isk)
443 nor1=sqrt(er(1)*er(1)+er(2)*er(2)+er(3)*er(3))
444 IF (nor1 > em20) er=er/nor1
445 er=er/max(nor1,em20)
446 et(1) =cth*skew(4,isk)- sth*skew(1,isk)
447 et(2) =cth*skew(5,isk)- sth*skew(2,isk)
448 et(3) =cth*skew(6,isk)- sth*skew(3,isk)
449 nor2=sqrt(et(1)*et(1)+et(2)*et(2)+et(3)*et(3))
450 et=et/max(nor2,em20)
451c-----------------vitesse angulaire imposee---------
452 alpha=-yc(ii)*dt2_double/two
453c-----------------sinon-----------------------------
454c ALPHA=-YC(II)*DT2_DOUBLE/(TWO*R)
455 skdir=alpha*er+et
456 ELSEIF(j == 3) THEN
457 skdir(1)=skew(7,isk)
458 skdir(2)=skew(8,isk)
459 skdir(3)=skew(9,isk)
460 IF(ibfv(7,n) == 2) yimp= skdir(1)*d(1,i) +
461 & skdir(2)*d(2,i) +
462 & skdir(3)*d(3,i)
463 ENDIF
464 nor3 = sqrt(skdir(1)*skdir(1)
465 & +skdir(2)*skdir(2)
466 & +skdir(3)*skdir(3))
467 skdir=skdir/max(nor3,em20)
468 vv = v(1,i)*skdir(1) + v(2,i)*skdir(2) + v(3,i)*skdir(3)
469 a0 = a(1,i)*skdir(1) + a(2,i)*skdir(2) + a(3,i)*skdir(3)
470 IF(ibfv(7,n) == 2) yc(ii)=(yc(ii)-yimp)/dt2_double
471 IF(ibfv(7,n)>=1) THEN
472 IF (j == 2) THEN
473c YC(II)=(YC(II)-VV)/DT12
474c-----------------vitesse angulaire imposee--------
475 yc(ii)=(r*yc(ii)-vv)/dt12
476c----------------------------------------------------
477 ELSE
478 yc(ii)=(yc(ii)-vv)/dt12
479 ENDIF
480 ENDIF
481 aa = yc(ii) - a0
482 IF(r<=em20) aa=zero
483 a(1,i)=a(1,i)+skdir(1)*aa
484 a(2,i)=a(2,i)+skdir(2)*aa
485 a(3,i)=a(3,i)+skdir(3)*aa
486 dw = fourth*ms(i)*(yc(ii)*dt12+two*vv)*aa*axi*weight(i)
487 ELSEIF (ifm > 1.AND.icoor == 0) THEN
488C Fixed velocity in moving frame (transl)
489 k1=3*j-2
490 k2=3*j-1
491 k3=3*j
492 rx = x(1,i) - xframe(10,ifm)
493 ry = x(2,i) - xframe(11,ifm)
494 rz = x(3,i) - xframe(12,ifm)
495 IF (ibfv(7,n) == 2) THEN
496 qf0(1:9) = xframe(19:27,ifm)
497 ej(1) = qf0(k1)
498 ej(2) = qf0(k2)
499 ej(3) = qf0(k3)
500 dd = ej(1)*d(1,i) +ej(2)*d(2,i) + ej(3)*d(3,i)
501 ddf= ej(1)*xframe(34,ifm)+ej(2)*xframe(35,ifm)+ej(3)*xframe(36,ifm)
502 dfj= (ej(1)-xframe(k1,ifm))*rx+(ej(2)-xframe(k2,ifm))*ry+
503 . (ej(3)-xframe(k3,ifm))*rz + ddf
504 IF (nsens(ii) > 0 .and. ts(ii) > zero) THEN
505 yc(ii) = (yc(ii)-yc0(ii))/dt2_double
506 ELSE
507 yc(ii) = (yc(ii)-dd+dfj)/dt2_double
508 ENDIF
509 vf = zero
510 ELSE
511 vfx = xframe(31,ifm)+xframe(14,ifm)*rz-xframe(15,ifm)*ry
512 vfy = xframe(32,ifm)+xframe(15,ifm)*rx-xframe(13,ifm)*rz
513 vfz = xframe(33,ifm)+xframe(13,ifm)*ry-xframe(14,ifm)*rx
514 vf = xframe(k1,ifm)*vfx+xframe(k2,ifm)*vfy+xframe(k3,ifm)*vfz
515 ENDIF
516 vv = xframe(k1,ifm)*v(1,i)
517 . + xframe(k2,ifm)*v(2,i)
518 . + xframe(k3,ifm)*v(3,i)
519 a0 = xframe(k1,ifm)*a(1,i)
520 . + xframe(k2,ifm)*a(2,i)
521 . + xframe(k3,ifm)*a(3,i)
522 yc(ii)=(yc(ii)-vv+vf)/dt12
523 aa = yc(ii) - a0
524 a(1,i)=a(1,i)+xframe(k1,ifm)*aa
525 a(2,i)=a(2,i)+xframe(k2,ifm)*aa
526 a(3,i)=a(3,i)+xframe(k3,ifm)*aa
527 dw = fourth*ms(i)*(yc(ii)*dt12+two*vv)*aa*axi*weight(i)
528 ELSEIF (ifm > 1.AND.icoor == 1) THEN
529 coef = (xframe(10,ifm)-x(1,i))*xframe(7,ifm) +
530 & (xframe(11,ifm)-x(2,i))*xframe(8,ifm) +
531 & (xframe(12,ifm)-x(3,i))*xframe(9,ifm)
532 hx = coef*xframe(7,ifm)+x(1,i)-xframe(10,ifm)
533 hy = coef*xframe(8,ifm)+x(2,i)-xframe(11,ifm)
534 hz = coef*xframe(9,ifm)+x(3,i)-xframe(12,ifm)
535 r = sqrt(hx*hx+hy*hy+hz*hz)
536 IF(r<=em20) THEN
537 cth = zero
538 sth = zero
539 ELSE
540 cth = (xframe(1,ifm)*hx +
541 & xframe(2,ifm)*hy +
542 & xframe(3,ifm)*hz)/r
543 sth = (xframe(4,ifm)*hx +
544 & xframe(5,ifm)*hy +
545 & xframe(6,ifm)*hz)/r
546 ENDIF
547 IF (j == 1) THEN
548 fmdir(1) =cth*xframe(1,ifm)+ sth*xframe(4,ifm)
549 fmdir(2) =cth*xframe(2,ifm)+ sth*xframe(5,ifm)
550 fmdir(3) =cth*xframe(3,ifm)+ sth*xframe(6,ifm)
551 ELSEIF (j == 2) THEN
552 er(1) =cth*xframe(1,ifm)+ sth*xframe(4,ifm)
553 er(2) =cth*xframe(2,ifm)+ sth*xframe(5,ifm)
554 er(3) =cth*xframe(3,ifm)+ sth*xframe(6,ifm)
555 nor1=sqrt(er(1)*er(1)+er(2)*er(2)+er(3)*er(3))
556 er=er/max(nor1,em20)
557 et(1) =cth*xframe(4,ifm)- sth*xframe(1,ifm)
558 et(2) =cth*xframe(5,ifm)- sth*xframe(2,ifm)
559 et(3) =cth*xframe(6,ifm)- sth*xframe(3,ifm)
560 nor2=sqrt(et(1)*et(1)+et(2)*et(2)+et(3)*et(3))
561 et=et/max(nor2,em20)
562c---------------vitesse angulaire imposee--------
563 alpha=-yc(ii)*dt2_double/two
564c-----------------------------------------------
565c ALPHA=-YC(II)*DT2_DOUBLE/(TWO*R)
566 fmdir=alpha*er+et
567 ELSEIF (j == 3) THEN
568 fmdir(1)=xframe(7,ifm)
569 fmdir(2)=xframe(8,ifm)
570 fmdir(3)=xframe(9,ifm)
571 ENDIF
572 nor3 = sqrt(fmdir(1)*fmdir(1)
573 & +fmdir(2)*fmdir(2)
574 & +fmdir(3)*fmdir(3))
575 fmdir=fmdir/max(nor3,em20) ! bug
576 rx = x(1,i) - xframe(10,ifm)
577 ry = x(2,i) - xframe(11,ifm)
578 rz = x(3,i) - xframe(12,ifm)
579 vfx = xframe(31,ifm)+xframe(14,ifm)*rz-xframe(15,ifm)*ry
580 vfy = xframe(32,ifm)+xframe(15,ifm)*rx-xframe(13,ifm)*rz
581 vfz = xframe(33,ifm)+xframe(13,ifm)*ry-xframe(14,ifm)*rx
582 vv = fmdir(1)*v(1,i) + fmdir(2)*v(2,i) + fmdir(3)*v(3,i)
583 vf = fmdir(1)*vfx + fmdir(2)*vfy + fmdir(3)*vfz
584 a0 = fmdir(1)*a(1,i) + fmdir(2)*a(2,i) + fmdir(3)*a(3,i)
585 IF (j == 2) THEN
586c YC(II)=(YC(II)-VV+VF)/DT12
587c---------------vitesse angulaire imposee--------
588 yc(ii)=(r*yc(ii)-vv+vf)/dt12
589c-------------------------------------------------
590 ELSE
591 yc(ii)=(yc(ii)-vv+vf)/dt12
592 ENDIF
593 aa = yc(ii) - a0
594 IF(r<=em20) aa=zero
595 a(1,i)=a(1,i)+fmdir(1)*aa
596 a(2,i)=a(2,i)+fmdir(2)*aa
597 a(3,i)=a(3,i)+fmdir(3)*aa
598 dw = fourth*ms(i)*(yc(ii)*dt12+two*vv)*aa*axi*weight(i)
599 ENDIF
600 ELSEIF(j<=6)THEN
601 j = j - 3
602 IF(isk<=1.AND.ifm<=1.AND.icoor == 0)THEN
603 IF (ibfv(7,n) == 2) THEN
604 IF (nsens(ii) > 0 .and. ts(ii) > zero) THEN
605 yc(ii) = (yc(ii)-yc0(ii))/dt2_double
606 ELSE
607 yc(ii)=(yc(ii)-dr(j,i))/dt2_double
608 ENDIF
609 ENDIF
610 IF (ibfv(7,n)>=1) yc(ii)=(yc(ii)-vr(j,i))/dt12
611 aold = ar(j,i)
612 ar(j,i)= yc(ii)
613 IF(ibfv(6,n) == 0)THEN
614 dw=fourth *
615 . (ar(j,i)*dt12 + two*vr(j,i))*in(i)*(ar(j,i)-aold) *
616 . weight(i)
617 ELSE
618 nr = ibfv(6,n)
619 dw= fourth *
620 . (ar(j,i)*dt12 + two*vr(j,i))*weight(i)*(ar(j,i)-aold)*
621 . (rby(16+j,nr) + rby(19+j,nr) + rby(22+j,nr))
622 ENDIF
623 ELSEIF (isk > 1.AND.icoor == 0) THEN
624 k1=3*j-2
625 k2=3*j-1
626 k3=3*j
627 IF(ibfv(7,n) == 2) THEN
628 dd = skew(k1,isk)*dr(1,i) +
629 . skew(k2,isk)*dr(2,i) +
630 . skew(k3,isk)*dr(3,i)
631 END IF
632 vv = skew(k1,isk)*vr(1,i) +
633 . skew(k2,isk)*vr(2,i) +
634 . skew(k3,isk)*vr(3,i)
635 a0 = skew(k1,isk)*ar(1,i) +
636 . skew(k2,isk)*ar(2,i) +
637 . skew(k3,isk)*ar(3,i)
638 IF (ibfv(7,n) == 2) yc(ii)=(yc(ii)-dd)/dt2_double
639 IF (ibfv(7,n)>=1) yc(ii)=(yc(ii)-vv)/dt12
640 aa = yc(ii) - a0
641 ar(1,i)=ar(1,i)+skew(k1,isk)*aa
642 ar(2,i)=ar(2,i)+skew(k2,isk)*aa
643 ar(3,i)=ar(3,i)+skew(k3,isk)*aa
644 IF(ibfv(6,n) == 0)THEN
645 dw= fourth*(yc(ii)*dt12+two*vv)*in(i)*aa*weight(i)
646 ELSE
647 nr = ibfv(6,n)
648 dw= fourth*(yc(ii)*dt12+two*vv)*weight(i)*aa*
649 . ((rby(17,nr)*skew(k1,isk)
650 . +rby(18,nr)*skew(k2,isk)
651 . +rby(19,nr)*skew(k3,isk))*skew(k1,isk) +
652 . (rby(20,nr)*skew(k1,isk)
653 . +rby(21,nr)*skew(k2,isk)
654 . +rby(22,nr)*skew(k3,isk))*skew(k2,isk) +
655 . (rby(23,nr)*skew(k1,isk)
656 . +rby(24,nr)*skew(k2,isk)
657 . +rby(25,nr)*skew(k3,isk))*skew(k3,isk))
658 ENDIF
659 ELSEIF ((isk<=1.AND.ifm<=1.AND.icoor == 1) .OR.
660 & (isk > 1.AND.icoor == 1)) THEN
661 coef = (skew(10,isk)-x(1,i))*skew(7,isk) +
662 & (skew(11,isk)-x(2,i))*skew(8,isk) +
663 & (skew(12,isk)-x(3,i))*skew(9,isk)
664 hx = coef*skew(7,isk)+x(1,i)-skew(10,isk)
665 hy = coef*skew(8,isk)+x(2,i)-skew(11,isk)
666 hz = coef*skew(9,isk)+x(3,i)-skew(12,isk)
667 r = sqrt(hx*hx+hy*hy+hz*hz)
668 IF(r<=em20) THEN
669 cth = zero
670 sth = zero
671 ELSE
672 cth = (skew(1,isk)*hx+skew(2,isk)*hy+skew(3,isk)*hz)/r
673 sth = (skew(4,isk)*hx+skew(5,isk)*hy+skew(6,isk)*hz)/r
674 ENDIF
675 IF (j == 1) THEN
676 skdir(1) =cth*skew(1,isk)+ sth*skew(4,isk)
677 skdir(2) =cth*skew(2,isk)+ sth*skew(5,isk)
678 skdir(3) =cth*skew(3,isk)+ sth*skew(6,isk)
679 ELSEIF(j == 2) THEN
680 skdir(1) =cth*skew(4,isk)- sth*skew(1,isk)
681 skdir(2) =cth*skew(5,isk)- sth*skew(2,isk)
682 skdir(3) =cth*skew(6,isk)- sth*skew(3,isk)
683 ELSEIF(j == 3) THEN
684 skdir(1)=skew(7,isk)
685 skdir(2)=skew(8,isk)
686 skdir(3)=skew(9,isk)
687 IF(ibfv(7,n) == 2) yimp= skdir(1)*dr(1,i) +
688 & skdir(2)*dr(2,i) +
689 & skdir(3)*dr(3,i)
690 ENDIF
691 nor3 = sqrt(skdir(1)*skdir(1)
692 & +skdir(2)*skdir(2)
693 & +skdir(3)*skdir(3))
694 skdir=skdir/max(nor3,em20)
695 IF(ibfv(7,n) == 2) yc(ii)=(yc(ii)-yimp)/dt2_double
696 vv =vr(1,i)*skdir(1)+ vr(2,i)*skdir(2)+ vr(3,i)*skdir(3)
697 a0 =ar(1,i)*skdir(1)+ ar(2,i)*skdir(2)+ ar(3,i)*skdir(3)
698 IF(ibfv(7,n)>=1) yc(ii)=(yc(ii)-vv)/dt12
699 aa = yc(ii) - a0
700 IF(r<=em20) aa=zero
701 ar(1,i)=ar(1,i)+skdir(1)*aa
702 ar(2,i)=ar(2,i)+skdir(2)*aa
703 ar(3,i)=ar(3,i)+skdir(3)*aa
704 IF(ibfv(6,n) == 0)THEN
705 dw= fourth*(yc(ii)*dt12+two*vv)*in(i)*aa*weight(i)
706 ELSE
707 nr = ibfv(6,n)
708 dw= fourth*(yc(ii)*dt12+two*vv)*weight(i)*aa*
709 & ((rby(17,nr)*skdir(1)
710 & +rby(18,nr)*skdir(2)
711 & +rby(19,nr)*skdir(3))*skdir(1) +
712 & (rby(20,nr)*skdir(1)
713 & +rby(21,nr)*skdir(2)
714 & +rby(22,nr)*skdir(3))*skdir(2) +
715 & (rby(23,nr)*skdir(1)
716 & +rby(24,nr)*skdir(2)
717 & +rby(25,nr)*skdir(3))*skdir(3))
718 ENDIF
719 ELSEIF (ifm > 1.AND.icoor == 0) THEN
720C Fixed velocity in moving frame (rota)
721 k1=3*j-2
722 k2=3*j-1
723 k3=3*j
724 IF (ibfv(7,n) == 2) THEN
725 qf0(1:9) = xframe(19:27,ifm)
726 ej(1) = qf0(k1)
727 ej(2) = qf0(k2)
728 ej(3) = qf0(k3)
729 ddr = ej(1)*dr(1,i) +ej(2)*dr(2,i) + ej(3)*dr(3,i)
730 ddf= ej(1)*xframe(13,ifm)+ej(2)*xframe(14,ifm)+ej(3)*xframe(15,ifm)
731 IF (nsens(ii) > 0 .and. ts(ii) > zero) THEN
732 yc(ii) = (yc(ii)-yc0(ii))/dt2_double
733 ELSE
734 yc(ii) = (yc(ii)-ddr)/dt2_double + ddf
735 ENDIF
736 vf = zero
737 ELSE
738 vf = xframe(k1,ifm)*xframe(13,ifm)
739 . + xframe(k2,ifm)*xframe(14,ifm)
740 . + xframe(k3,ifm)*xframe(15,ifm)
741 ENDIF
742 vv = xframe(k1,ifm)*vr(1,i)
743 . + xframe(k2,ifm)*vr(2,i)
744 . + xframe(k3,ifm)*vr(3,i)
745 a0 = xframe(k1,ifm)*ar(1,i)
746 . + xframe(k2,ifm)*ar(2,i)
747 . + xframe(k3,ifm)*ar(3,i)
748 aa = (yc(ii)-vv+vf)/dt12 - a0
749 ar(1,i)=ar(1,i)+xframe(k1,ifm)*aa
750 ar(2,i)=ar(2,i)+xframe(k2,ifm)*aa
751 ar(3,i)=ar(3,i)+xframe(k3,ifm)*aa
752 IF(ibfv(6,n) == 0)THEN
753 dw= fourth*(yc(ii)+vv)*in(i)*aa*weight(i)
754 ELSE
755 nr = ibfv(6,n)
756 dw= fourth*(yc(ii)+vv)*weight(i)*aa*
757 . ((rby(17,nr)*xframe(k1,ifm)
758 . +rby(18,nr)*xframe(k2,ifm)
759 . +rby(19,nr)*xframe(k3,ifm))*xframe(k1,ifm) +
760 . (rby(20,nr)*xframe(k1,ifm)
761 . +rby(21,nr)*xframe(k2,ifm)
762 . +rby(22,nr)*xframe(k3,ifm))*xframe(k2,ifm) +
763 . (rby(23,nr)*xframe(k1,ifm)
764 . +rby(24,nr)*xframe(k2,ifm)
765 . +rby(25,nr)*xframe(k3,ifm))*xframe(k3,ifm))
766 ENDIF
767 ELSEIF (ifm > 1.AND.icoor == 1) THEN
768 coef = (xframe(10,ifm)-x(1,i))*xframe(7,ifm) +
769 & (xframe(11,ifm)-x(2,i))*xframe(8,ifm) +
770 & (xframe(12,ifm)-x(3,i))*xframe(9,ifm)
771 hx = coef*xframe(7,ifm)+x(1,i)-xframe(10,ifm)
772 hy = coef*xframe(8,ifm)+x(2,i)-xframe(11,ifm)
773 hz = coef*xframe(9,ifm)+x(3,i)-xframe(12,ifm)
774 r = sqrt(hx*hx+hy*hy+hz*hz)
775 IF(r<=em20) THEN
776 cth = zero
777 sth = zero
778 ELSE
779 cth = (xframe(1,ifm)*hx +
780 & xframe(2,ifm)*hy +
781 & xframe(3,ifm)*hz)/r
782 sth = (xframe(4,ifm)*hx +
783 & xframe(5,ifm)*hy +
784 & xframe(6,ifm)*hz)/r
785 ENDIF
786 IF (j == 1) THEN
787 fmdir(1) =cth*xframe(1,ifm)+ sth*xframe(4,ifm)
788 fmdir(2) =cth*xframe(2,ifm)+ sth*xframe(5,ifm)
789 fmdir(3) =cth*xframe(3,ifm)+ sth*xframe(6,ifm)
790 ELSEIF (j == 2) THEN
791 fmdir(1) =cth*xframe(4,ifm)- sth*xframe(1,ifm)
792 fmdir(2) =cth*xframe(5,ifm)- sth*xframe(2,ifm)
793 fmdir(3) =cth*xframe(6,ifm)- sth*xframe(3,ifm)
794 ELSEIF (j == 3) THEN
795 fmdir(1)=xframe(7,ifm)
796 fmdir(2)=xframe(8,ifm)
797 fmdir(3)=xframe(9,ifm)
798 ENDIF
799 nor3 = sqrt(fmdir(1)*fmdir(1)
800 & +fmdir(2)*fmdir(2)
801 & +fmdir(3)*fmdir(3))
802 fmdir=fmdir/max(nor3,em20)
803 vv = fmdir(1)*vr(1,i) +
804 & fmdir(2)*vr(2,i) +
805 & fmdir(3)*vr(3,i)
806 vf = fmdir(1)*xframe(13,ifm) +
807 & fmdir(2)*xframe(14,ifm) +
808 & fmdir(3)*xframe(15,ifm)
809 a0 = fmdir(1)*ar(1,i) +
810 & fmdir(2)*ar(2,i) +
811 & fmdir(3)*ar(3,i)
812 aa = (yc(ii)-vv+vf)/dt12 - a0
813 IF(r<=em20) aa=zero
814 ar(1,i)=ar(1,i)+fmdir(1)*aa
815 ar(2,i)=ar(2,i)+fmdir(2)*aa
816 ar(3,i)=ar(3,i)+fmdir(3)*aa
817 IF(ibfv(6,n) == 0)THEN
818 dw= fourth*(yc(ii)+vv)*in(i)*aa*weight(i)
819 ELSE
820 nr = ibfv(6,n)
821 dw= fourth*(yc(ii)+vv)*weight(i)*aa*
822 & ((rby(17,nr)*fmdir(1)
823 & +rby(18,nr)*fmdir(2)
824 & +rby(19,nr)*fmdir(3))*fmdir(1) +
825 & (rby(20,nr)*fmdir(1)
826 & +rby(21,nr)*fmdir(2)
827 & +rby(22,nr)*fmdir(3))*fmdir(2) +
828 & (rby(23,nr)*fmdir(1)
829 & +rby(24,nr)*fmdir(2)
830 & +rby(25,nr)*fmdir(3))*fmdir(3))
831 ENDIF
832 ENDIF
833 ENDIF
834 wfext=wfext + dt1*dw
835C DT2_DOUBLE*dw into memory ; if sms on the degree of freedom,
836C part of DT2_DOUBLE*dw was already computed and stored in sms_fixvel...
837 vel(4,n) = vel(4,n)+dt2_double*dw
838 ENDDO
839c----------------------------------------------------------------------------
840 ELSE ! partie vectorielle
841c----------------------------------------------------------------------------
842 ENDIF
843 ENDDO ! NN=1,NFXVEL
844c-----------
845 RETURN
846 END
alpha
#define alpha
Definition eval.h:35
fixvel
subroutine fixvel(ibfv, a, v, npc, tf, vel, ms, x, skew, ar, vr, in, nsensor, sensor_tab, weight, depla, rby, iframe, xframe, dr, nodnx_sms, nodes, tt_double, depla_double, python, wfext)
Definition fixvel.F:42
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
vinter_smooth
subroutine vinter_smooth(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter_smooth.F:30
vinterdp
subroutine vinterdp(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinterdp.F:29