51
52
53
56 USE pblast_mod
57 USE sensor_mod
59 USE loads_mod
62 USE pblast_mod
63
64
65
66#include "implicit_f.inc"
67#include "param_c.inc"
68
69
70
71#include "com04_c.inc"
72#include "com08_c.inc"
73#include "tabsiz_c.inc"
74
75
76
77 INTEGER GET_U_NUMSENS,GET_U_SENS_FPAR,GET_U_SENS_IPAR,
78 . GET_U_SENS_VALUE,SET_U_SENS_VALUE
81
82
83
84 INTEGER ,INTENT(IN) :: NSENSOR
85 CHARACTER(LEN=NCHARLINE100):: KEYWORD
87 . skin_vector(3,*),tf(*),x(3,*),v(3,*),dis(3,numnod)
88 TYPE (H3D_DATABASE) :: H3D_DATA
89 INTEGER , DIMENSION(NUMSKINP0), INTENT(IN) :: IMAPSKP
90 integer
91 . h3d_part(*),is_written_skin(*),info1,npc(*)
92 INTEGER LLOADP(SLLOADP)
93 INTEGER ILOADP(SIZLOADP,*),IB(NIBCLD,*)
94 INTEGER TAGNCONT(NLOADP_HYD_INTER,NUMNOD),
95 . LOADP_HYD_INTER(NLOADP_HYD),NODAL_IPART(*)
97 . fac(lfacload,nloadp),xframe(nxframe,*),forc(lfaccld,*)
98 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
99 TYPE (LOADS_) , INTENT(IN) :: LOADS
100 INTEGER , DIMENSION(LISKN,NUMFRAM+1) ,INTENT(IN) :: IFRAME
101 TYPE (TTABLE) ,DIMENSION(NTABLE) ,INTENT(IN) :: TABLE
102 TYPE(PBLAST_),INTENT(IN) :: PBLAST
103
104
105
106 INTEGER NL, N1, ISK, N2, N3, N4, ,NSKIN,
107 . IAD ,NP ,IFUNC ,NPRES,NSKIN0,NSKIN1,N1FRAM,DIR_HSP,I,N
108 INTEGER K1, K2, K3, ISENS,K,LL,IERR,
109 . N_OLD, ISMOOTH,IDEL,NINTERP ,NPL,TAGN1,TAGN2,TAGN3,TAGN4,
110 . FUN_CX,FUN_VEL,DIR_VEL,IFRA2, IANIM,IJK,UP_BOUND,
111 . IZ_UPDATE,ABAC_ID,ISIZ_SEG,IERR1,
112 . Phi_I, ID, USER_ID, ITA_SHIFT,NDT,NDT0,
113 . NITER,ITER,IMODEL,IL,IS,SEGCONT,FUN_HSP,IFRA1,NP0
115 . nx, ny, nz, axi, aa, a0, vv, fx, fy, fz, ax, dydx, ts,
116 . sixth,x_old, f1, f2,xsens,fcx,fcy,fcypinch,fp,
117 . fcx1,fcy1,fcx2,fcy2,vx,vy,vz,vel,vseg,
norm
118 my_real finter, ps, zx,zy,zz,finter_smooth
120 . rxi,ryi,rzi,sxi,syi,szi
121 my_real coormean,ymean,zmean,pvel,nsign,dnorm,
122 . xdet,ydet,zdet,tdet,wtnt,pmin,dx,dy,dz,normd, p,
123 . t0inf_loc, ta_shift, tt_star
124 INTEGER :: IFUN,IFRA,M1,M2,NDIM,NPOINT, IIOUT,SHIFT,FUNCTYPE
125 my_real :: a11,a12,a21,a22,b1,b2,det,len,dirx,diry,dirz,
126 . beta,gamma,r,s,rmax,xfacr,xfact,yfac,segp,press,disp
127 my_real,
DIMENSION(3) :: p0,dir,a,b,c,d,m
128
129 EXTERNAL finter,finter_smooth
130
131 IF (keyword /= 'VECT/PEXT') RETURN
133 is_written_skin(nskin+1:numskin) = 0
134 skin_vector(1:3,nskin+1:numskin)=zero
135 np0= 0
136 nskin0 = nskin
137
138
139 n_old = 0
140 x_old = zero
151 IF (n1==0.OR.n2==0.OR.n3==0.OR.n4==-1) cycle
152
153 np0 = np0 + 1
154 nskin = nskin0+ imapskp(np0)
155 IF (nodal_ipart(n1)>0) THEN
156 IF (h3d_part(nodal_ipart(n1))==1) is_written_skin(nskin)=1
157 END IF
158 isens = 0
159 xsens = one
160 DO k=1,nsensor
161 IF(ib(6,
nl)== sensor_tab(k)%SENS_ID) isens=k
162 ENDDO
163 IF(isens==0)THEN
164 ts=tt
165 ELSE
166 ts = tt- sensor_tab(isens)%TSTART
167 IF(ts < zero) cycle
168 ENDIF
169 IF(idel > 0 ) cycle
170 IF(functype == 1)THEN
171 IF(n_old/=n5.OR.x_old/=ts) THEN
172 ismooth = 0
173 IF (n5 > 0) ismooth = npc(2*nfunct+n5+1)
174 IF (ismooth == 0) THEN
175 f1 = finter(n5,ts*fcx,npc,tf,dydx)
176 ELSE IF(ismooth > 0) THEN
177 f1 = finter_smooth(n5,ts*fcx,npc,tf,dydx)
178 ELSE
179 ENDIF
180 ENDIF
181 n_old = n5
182 x_old = ts
183 ELSEIF(functype == 2)THEN
184 IF(n_old/=n5) THEN
185 disp = (dis(3,n1)+dis(3,n2)+dis(3,n3)+dis(3,n4))/4.0
186 ismooth = 0
187 IF (n5 > 0) ismooth = npc(2*nfunct+n5+1)
188 IF (ismooth == 0) THEN
189 f1 = finter(n5,disp*fcx,npc,tf,dydx)
190 ELSE IF(ismooth > 0) THEN
191 f1 = finter_smooth(n5,disp*fcx,npc,tf,dydx)
192 ELSE
193 ENDIF ! IF (ismooth == 0)
194 ENDIF
195 n_old = n5
196 x_old = disp
197 ELSEIF(functype == 3)THEN
198 IF(n_old/=n5) THEN
199 vel = (v(3,n1)+v(3,n2)+v(3,n3)+v(3,n4))/4.0
200 ismooth = 0
201 IF (n5 > 0) ismooth = npc(2*nfunct+n5+1)
202 IF (ismooth == 0) THEN
203 f1 = finter(n5,vel*fcx,npc,tf,dydx)
204 ELSE IF(ismooth > 0) THEN
205 f1 = finter_smooth(n5,vel*fcx,npc,tf,dydx)
206 ELSE
207 ENDIF
208 ENDIF
209 n_old = n5
210 x_old = vel
211 ENDIF
212 aa = fcy*f1*xsens
213 skin_vector(3,nskin)=skin_vector(3,nskin) + aa
214 END DO
215
216 shift = nloadp_f+pblast%NLOADP_B
217 DO np=1+shift,nloadp_hyd+shift
218 isiz_seg = iloadp(1,np)/4
219 ifunc = iloadp(3,np)
220 iad = iloadp(4,np)
222 isens = iloadp(7,np)
223 fcy = fac(1,np)
224 fcx = fac(2,np)
225
226 DO n=1, isiz_seg
227 n1 = lloadp(iad+4*(n-1))
228 n2 = lloadp(iad+4*(n-1)+1)
229 n3 = lloadp(iad+4*(n-1)+2)
230 n4 = lloadp(iad+4*(n-1)+3)
231 IF (n1==0.OR.n2==0.OR.n3==0) cycle
232 nskin = nskin0+ imapskp(np0+n)
233 IF (nodal_ipart(n1)>0) THEN
234 IF (h3d_part(nodal_ipart(n1))==1) is_written_skin(nskin)=1
235 END IF
236 ENDDO
237
238 IF(isens==0)THEN
239 ts=tt
240 ELSE
241 ts = tt-sensor_tab(isens)%TSTART
242 ENDIF
243 DO n=1, isiz_seg
244 n1 = lloadp(iad+4*(n-1))
245 n2 = lloadp(iad+4*(n-1)+1)
246 n3 = lloadp(iad+4*(n-1)+2)
247 n4 = lloadp(iad+4*(n-1)+3)
248 IF (n1==0.OR.n2==0.OR.n3==0) cycle
249
250 np0 = np0 + 1
251 IF(ts<zero) cycle
252 nskin = nskin0+ imapskp(np0)
253 f1 = finter(ifunc,ts*fcx,npc,tf,dydx)
254 aa = fcy*f1
255
256
257
258 segcont = 0
259
260 tagn1 = 0
261 tagn2 = 0
262 tagn3 = 0
263 tagn4 = 0
264 fp = one
266 npl = loadp_hyd_inter(np)
267 IF(n4/=0) THEN
268 IF(tagncont(npl,n1)==1.AND.tagncont(npl,n2)==1.AND.
269 . tagncont(npl,n3)==1.AND.tagncont(npl,n4)==1) THEN
270 segcont = 1
271 ELSE
272 tagn1 = tagncont(npl,n1)
273 tagn2 = tagncont(npl,n2)
274 tagn3 = tagncont(npl,n3)
275 tagn4 = tagncont(npl,n4)
276 fp = (tagn1+tagn2+tagn3+tagn4)/4
277 ENDIF
278 ELSE
279 IF(tagncont(npl,n1)==1.AND.tagncont(npl,n2)==1.AND.
280 . tagncontTHEN
281 segcont = 1
282 ELSE
283 tagn1 = tagncont(npl,n1)
284 tagn2 = tagncont(npl,n2)
285 tagn3 = tagncont(npl,n3)
286 fp = (tagn1+tagn2+tagn3)/3
287 ENDIF
288 ENDIF
289 IF (fp==zero) fp = one
290 ENDIF
291 IF (segcont==1) aa = zero
292 aa=aa*fp
293 skin_vector(3,nskin)=skin_vector(3,nskin) + aa
294 END DO
295 END DO
296
298
299 isiz_seg = iloadp(1,
nl)/4
301 DO n=1, isiz_seg
302 n1 = lloadp(iad+4*(n-1))
303 n2 = lloadp(iad+4*(n-1)+1)
304 n3 = lloadp(iad+4*(n-1)+2)
305 n4 = lloadp(iad+4*(n-1)+3)
306 IF (n1==0.OR.n2==0.OR.n3==0) cycle
307 nskin = nskin0+ imapskp(np0+n)
308 IF (nodal_ipart(n1)>0) THEN
309 IF (h3d_part(nodal_ipart(n1))==1) is_written_skin(nskin)=1
310 END IF
311 ENDDO
320 fun_vel=iloadp(11,
nl)
323
324
325 dir_vel=
max(iloadp(12,
nl),1)
327 isens=0
328 xsens = one
329 DO k=1,nsensor
330 IF(iloadp(6,
nl)== sensor_tab(k)%SENS_ID) isens=k
331 ENDDO
332 IF(isens==0)THEN
333 ts=tt
334 ELSE
335 ts = tt-sensor_tab(isens)%TSTART
336 ENDIF
337 DO i = 1,isiz_seg
338 n1=lloadp(iloadp(4,
nl)+4*(i-1))
339 n2=lloadp(iloadp(4,
nl)+4*(i-1)+1)
340 n3=lloadp(iloadp(4,
nl)+4*(i-1)+2)
341 n4=lloadp(iloadp(4,
nl)+4*(i-1)+3)
342 IF (n1==0.OR.n2==0.OR.n3==0) cycle
343 np0 = np0 + 1
344 IF(ts < zero) cycle
345 nskin = nskin0+ imapskp(np0)
346
347 aa = zero
348 vel = zero
349 pvel=zero
350
351
352 IF(n4/=0 .AND. n1/=n2 .AND. n1/=n3 .AND. n1/=n4 .AND.
353 . n2/=n3 .AND. n2/=n4 .AND. n3/=n4 )THEN
354
355 k1=3*dir_hsp-2
356 k2=3*dir_hsp-1
357 k3=3*dir_hsp
358
359 IF(fun_hsp /=0)THEN
360 coormean = (xframe(k1,ifra1)*(x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4))/four)+
361 . (xframe(k2,ifra1)*(x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4))/four)+
362 . (xframe(k3,ifra1)*(x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4))/four)
363 aa = fcy*finter(fun_hsp,(coormean-xframe(9+dir_hsp,ifra1))*fcx,npc,tf,dydx)
364 ENDIF
365 nx= (x(2,n3)-x(2,n1))*(x(3,n4)-x(3,n2)) - (x(3,n3)-x(3,n1))*(x(2,n4)-x(2,n2))
366 ny= (x(3,n3)-x(3,n1))*(x(1,n4)-x(1,n2)) - (x(1,n3)-x(1,n1))*(x(3,n4)-x(3,n2))
367 nz= (x(1,n3)-x(1,n1))*(x(2,n4)-x(2,n2)) - (x(2,n3)-x(2,n1))*(x(1,n4)-x(1,n2))
368 norm = sqrt(nx*nx+ny*ny+nz*nz)
369 aa = aa * half *
norm
370
371 k1=3*dir_vel-2
372 k2=3*dir_vel-1
373 k3=3*dir_vel
374
375 nsign = (nx * xframe(k1,ifra2) +
376 . ny * xframe(k2,ifra2) +
377 . nz * xframe(k3,ifra2))
378 IF(nsign/=zero) nsign = sign(one,nsign)
379
380 vseg= (xframe(k1,ifra2)*
381 . (v(1,n1) + v(1,n2) + v(1,n3) + v(1,n4)) /four)+
382 . (xframe(k2,ifra2)*
383 . (v(2,n1) + v(2,n2) + v(2,n3) + v(2,n4)) /four)+
384 . (xframe(k3,ifra2)*
385 . (v(3,n1) + v(3,n2) + v(3,n3) + v(3,n4)) /four)
386
387 IF(fun_vel /=0)THEN
388 vel = fcy2*finter(fun_vel,tt*fcx2,npc,tf,dydx)- vseg
389 ELSE
390 vel = - vseg
391 ENDIF
392 IF(fun_cx /=0)
393 . pvel = ( (-(nx/
norm)*vel*xframe(k1,ifra2)-
394 . (ny/
norm)*vel*xframe(k2,ifra2)-
395 . (nz/
norm)*vel*xframe(k3,ifra2))**2 )* fcy1*
396 . finter(fun_cx,tt*fcx1,npc,tf,dydx)/two
397
398 ELSE
399 IF(n1 == n2)THEN
400 n2 = n3
401 n3 = n4
402 n4 = 0
403 ELSEIF(n1 == n3)THEN
404 n3 = n4
405 n4 = 0
406 ELSEIF(n1 == n4)THEN
407 n4 = 0
408 ELSEIF(n2 == n3)THEN
409 n3 = n4
410 n4 = 0
411 ELSEIF(n2 == n4)THEN
412 n2 = n3
413 n3 = n4
414 n4 = 0
415 ELSEIF(n3 == n4)THEN
416 n4 = 0
417 ENDIF
418 IF (n4==0) n4=n3
419
420 IF(fun_hsp /=0)THEN
421 k1=3*dir_hsp-2
422 k2=3*dir_hsp-1
423 k3=3*dir_hsp
424
425 coormean = (xframe(k1,ifra1)*(x(1,n1)+x(1,n2)+x(1,n3))/three)+
426 . (xframe(k2,ifra1)*(x(2,n1)+x(2,n2)+x(2,n3))/three)+
427 . (xframe(k3,ifra1)*(x(3,n1)+x(3,n2)+x(3,n3))/three)
428 aa = fcy*finter(fun_hsp,(coormean-xframe(9+dir_hsp,ifra1))*fcx,npc,tf,dydx)
429 ENDIF
430 nx= (x(2,n3)-x(2,n1))*(x(3,n4)-x(3,n2)) - (x(3,n3)-x(3,n1))*(x(2,n4)-x(2,n2))
431 ny= (x(3,n3)-x(3,n1))*(x(1,n4)-x(1,n2)) - (x(1,n3)-x(1,n1))*(x(3,n4)-x(3,n2))
432 nz= (x(1,n3)-x(1,n1))*(x(2,n4)-x(2,n2)) - (x(2,n3)-x(2,n1))*(x(1,n4)-x(1,n2))
433 norm = sqrt(nx*nx+ny*ny+nz*nz)
434 aa = aa * half *
norm
435
436 k1=3*dir_vel-2
437 k2=3*dir_vel-1
438 k3=3*dir_vel
439
440 nsign = (nx * xframe(k1,ifra2) +
441 . ny * xframe(k2,ifra2) +
442 . nz * xframe(k3,ifra2))
443 IF(nsign/=zero) nsign = sign(one,nsign)
444
445 vseg= (xframe(k1,ifra2)*
446 . (v(1,n1) + v(1,n2) + v(1,n3)) /three)+
447 . (xframe(k2,ifra2)*
448 . (v(2,n1) + v(2,n2) + v(2,n3)) /three)+
449 . (xframe(k3,ifra2)*
450 . (v(3,n1) + v(3,n2) + v(3,n3)) /three)
451
452 IF(fun_vel /=0)THEN
453 vel = fcy2*finter(fun_vel,tt*fcx2,npc,tf,dydx)- vseg
454 ELSE
455 vel = - vseg
456 ENDIF
457 IF(fun_cx /=0)
458 . pvel = ( (-(nx/
norm)*vel*xframe(k1,ifra2)-
459 . (ny/
norm)*vel*xframe(k2,ifra2)-
460 . (nz/
norm)*vel*xframe(k3,ifra2))**2 )* fcy1*
461 . finter(fun_cx,tt*fcx1,npc,tf,dydx)/two
462 ENDIF
463 aa=-aa+pvel*nsign
464 skin_vector(3,nskin)=skin_vector(3,nskin) + aa
465 END DO
466 END DO
467
468 DO nl=1+nloadp_f,nloadp_f+pblast%NLOADP_B
469
470 isiz_seg = iloadp(1,
nl)/4
472 DO n=1, isiz_seg
473 n1 = lloadp(iad+4*(n-1))
474 n2 = lloadp(iad+4*(n-1)+1)
475 n3 = lloadp(iad+4*(n-1)+2)
476 n4 = lloadp(iad+4*(n-1)+3)
477 IF (n1==0.OR.n2==0.OR.n3==0) cycle
478 nskin = nskin0+ imapskp(np0+n)
479 IF (nodal_ipart(n1)>0) THEN
480 IF (h3d_part(nodal_ipart(n1))==1) is_written_skin(nskin)=1
481 END IF
482 ENDDO
486
487
488
489 DO i = 1,isiz_seg
490 n1=lloadp(iloadp(4,
nl)+4*(i-1))
491 n2=lloadp(iloadp(4,
nl)+4*(i-1)+1)
492 n3=lloadp(iloadp(4,
nl)+4*(i-1)+2)
493 n4=lloadp(iloadp(4,
nl)+4*(i-1)+3)
494 IF (n1==0.OR.n2==0.OR.n3==0) cycle
495 np0 = np0 + 1
496 IF(tt<tdet)cycle
497 nskin = nskin0+ imapskp(np0)
498 p = pblast%PBLAST_TAB(il)%PRES(i)
499 aa= -p
500 skin_vector(3,nskin)=skin_vector(3,nskin) + aa
501 ENDDO
502 END DO
503
504 DO nl=1,loads%NLOAD_CYL
505
506 isiz_seg = loads%LOAD_CYL(
nl)%NSEG
507 DO n=1, isiz_seg
508 n1 = loads%LOAD_CYL(
nl)%SEGNOD(n,1)
509 n2 = loads%LOAD_CYL(
nl)%SEGNOD(n,2)
510 n3 = loads%LOAD_CYL(
nl)%SEGNOD(n,3)
511 n4 = loads%LOAD_CYL(
nl)%SEGNOD(n,4)
512 IF (n1==0.OR.n2==0.OR.n3==0) cycle
513 nskin = nskin0+ imapskp(np0+n)
514 IF (nodal_ipart(n1)>0) THEN
515 IF (h3d_part(nodal_ipart(n1))==1) is_written_skin(nskin)=1
516 END IF
517 ENDDO
518 isens = loads%LOAD_CYL(
nl)%ISENS
519 iiout = 0
520 IF (isens > 0) THEN
521 IF (sensor_tab(isens)%STATUS == 0) THEN
522 np0 = np0 + isiz_seg
523 cycle
524 END IF
525 END IF
526 ifra = loads%LOAD_CYL(
nl)%IFRAME + 1
527 xfacr= loads%LOAD_CYL(
nl)%XSCALE_R
528 xfact= loads%LOAD_CYL(
nl)%XSCALE_T
529 yfac = loads%LOAD_CYL(
nl)%YSCALE
530 ifun = loads%LOAD_CYL(
nl)%ITABLE
531 ndim = table(ifun)%NDIM
532 npoint = SIZE(table(ifun)%X(1)%VALUES)
533 rmax = table(ifun)%X(1)%VALUES(npoint)
534 m1 = iframe(1,ifra)
535 m2 = iframe(2,ifra)
536 dirx = x(1,m1) - x(1,m2)
537 diry = x(2,m1) - x(2,m2)
538 dirz = x(3,m1) - x(3,m2)
539 len = sqrt(dirx**2 + diry**2 + dirz**2)
540
541 dir(1) = dirx / len
542 dir(2) = diry / len
543 dir(3) = dirz / len
544 p0(1) = x(1,m2)
545 p0(2)
546 p0(3) = x(3,m2)
547
548
549
550 DO n = 1,isiz_seg
551 n1 = loads%LOAD_CYL(
nl)%SEGNOD(n,1)
552 n2 = loads%LOAD_CYL(
nl)%SEGNOD(n,2)
553 n3 = loads%LOAD_CYL(
nl)%SEGNOD(n,3)
554 n4 = loads%LOAD_CYL(
nl)%SEGNOD(n,4)
555 press = zero
556 a(1) = x(1,n1)
557 a(2) = x(2,n1)
558 a(3) = x(3,n1)
559 b(1) = x(1,n2)
560 b(2) = x(2,n2)
561 b(3) = x(3,n2)
562 c(1) = x(1,n3)
563 c(2) = x(2,n3)
564 c(3) = x(3,n3)
565 np0 = np0 + 1
566 IF (n4 == 0) THEN
568 . ifun ,table ,xfacr ,xfact
569 press = abs(segp) * yfac
570
571 ELSE
572 d(1) = x
573 d(2) = x(2,n4)
574 d(3) = x(3,n4)
575 m(1) = (x(1,n1) + x(1,n2) + x(1,n3) + x(1,n4)) * fourth
576 m(2) = (x(2,n1) + x(2,n2) + x(2,n3
577 m(3) = (x(3,n1) + x(3,n2) + x(3,n3) + x(3,n4)) * fourth
578
580 . ifun ,table ,xfacr ,xfact ,segp )
581 press = press + segp * fourth
582
584 . ifun ,table ,xfacr ,xfact ,segp )
585 press = press + segp * fourth
586
588 . ifun ,table ,xfacr ,xfact ,segp )
589 press = press + segp * fourth
590
592 . ifun ,table ,xfacr ,xfact ,segp )
593 press = press
594 press = abs(press) * yfac
595 END IF
596 nskin = nskin0+ imapskp(np0)
597 skin_vector(3,nskin)=skin_vector(3,nskin) +press
598 enddo
599 END DO
600
601 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
integer, parameter ncharline100
subroutine ninterp(ifunc, npc, pld, npoint, xv, yv)
subroutine press_seg3(a, b, c, n1, dir, ifunc, table, xfacr, xfact, press)
character *2 function nl()
integer function set_u_sens_value(nsens, ivar, var)
integer function get_u_sens_ipar(nsens, ivar, var)
integer function get_u_sens_value(nsens, ivar, var)
integer function get_u_sens_fpar(nsens, ivar, var)
integer function get_u_numsens(idsens)