OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
h3d_skin_vector.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!|| h3d_skin_vector ../engine/source/output/h3d/h3d_results/h3d_skin_vector.f
25!||--- called by ------------------------------------------------------
26!|| genh3d ../engine/source/output/h3d/h3d_results/genh3d.F
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| finter_smooth ../engine/source/tools/curve/finter_smooth.F
30!|| get_u_numsens ../engine/source/user_interface/usensor.F
31!|| get_u_sens_fpar ../engine/source/user_interface/usensor.F
32!|| get_u_sens_ipar ../engine/source/user_interface/usensor.F
33!|| get_u_sens_value ../engine/source/user_interface/usensor.F
34!|| press_seg3 ../engine/source/loads/general/load_pcyl/press_seg3.F
35!|| set_u_sens_value ../engine/source/user_interface/usensor.F
36!||--- uses -----------------------------------------------------
37!|| h3d_inc_mod ../engine/share/modules/h3d_inc_mod.F
38!|| h3d_mod ../engine/share/modules/h3d_mod.F
39!|| loads_mod ../common_source/modules/loads/loads_mod.F90
40!|| names_and_titles_mod ../common_source/modules/names_and_titles_mod.F
41!|| pblast_mod ../common_source/modules/loads/pblast_mod.F90
42!|| pinchtype_mod ../common_source/modules/pinchtype_mod.F
43!|| sensor_mod ../common_source/modules/sensor_mod.F90
44!|| table_mod ../engine/share/modules/table_mod.F
45!||====================================================================
46 SUBROUTINE h3d_skin_vector(SKIN_VECTOR,NODAL_IPART,NSENSOR,
47 . IS_WRITTEN_SKIN ,H3D_PART,INFO1 ,KEYWORD ,
48 . IB ,ILOADP,LLOADP,FAC ,NPC,TF ,SENSOR_TAB,
49 . TAGNCONT,LOADP_HYD_INTER,FORC,XFRAME,X ,V ,
50 . IMAPSKP,LOADS ,TABLE,IFRAME,DIS,PBLAST)
51C-----------------------------------------------
52C M o d u l e s
53C-----------------------------------------------
54 USE h3d_mod
55 USE pinchtype_mod
56 USE pblast_mod
57 USE sensor_mod
58 USE h3d_inc_mod
59 USE loads_mod
60 USE table_mod
62 USE pblast_mod
63C-----------------------------------------------
64C I m p l i c i t T y p e s
65C-----------------------------------------------
66#include "implicit_f.inc"
67#include "param_c.inc"
68C-----------------------------------------------
69C C o m m o n B l o c k s
70C-----------------------------------------------
71#include "com04_c.inc"
72#include "com08_c.inc"
73#include "tabsiz_c.inc"
74C-----------------------------------------------
75C E x t e r n a l F u n c t i o n s
76C-----------------------------------------------
77 INTEGER GET_U_NUMSENS,GET_U_SENS_FPAR,GET_U_SENS_IPAR,
78 . GET_U_SENS_VALUE,SET_U_SENS_VALUE
79 EXTERNAL GET_U_NUMSENS,GET_U_SENS_FPAR,GET_U_SENS_IPAR,
80 . GET_U_SENS_VALUE,SET_U_SENS_VALUE
81C-----------------------------------------------,
82C D u m m y A r g u m e n t s
83C-----------------------------------------------
84 INTEGER ,INTENT(IN) :: NSENSOR
85 CHARACTER(LEN=NCHARLINE100):: KEYWORD
86 my_real
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(*)
96 my_real
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
103C-----------------------------------------------
104C L o c a l V a r i a b l e s
105C-----------------------------------------------
106 INTEGER NL, N1, ISK, N2, N3, N4, N5,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
114 my_real
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
119 my_real
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
130C=======================================================================
131 IF (keyword /= 'VECT/PEXT') RETURN
132 nskin = numskin - numskinp
133 is_written_skin(nskin+1:numskin) = 0
134 skin_vector(1:3,nskin+1:numskin)=zero
135 np0= 0
136 nskin0 = nskin
137C---- fill SKIN_SCALAR(*) w/ IS_WRITTEN_SKIN(*)=1
138C-----Force (pressure) first
139 n_old = 0
140 x_old = zero
141 DO nl=1,nconld-nploadpinch
142 n1 = ib(1,nl)
143 n2 = ib(2,nl)
144 n3 = ib(3,nl)
145 n4 = ib(4,nl)
146 n5 = ib(5,nl)
147 idel = ib(8,nl)
148 functype= ib(9,nl)
149 fcy = forc(1,nl)
150 fcx = forc(2,nl)
151 IF (n1==0.OR.n2==0.OR.n3==0.OR.n4==-1) cycle
152C--------default zero
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 ! SEGMENT DELETED
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 ! IF (ISMOOTH == 0)
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 ! IF (ISMOOTH == 0)
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
215C----------load_pressure
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)
221 ninterp = iloadp(5,np)
222 isens = iloadp(7,np)
223 fcy = fac(1,np)
224 fcx = fac(2,np)
225C--------default zero
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
237C
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
255C----------------
256C Check if segment is in contact
257C----------------
258 segcont = 0
259
260 tagn1 = 0
261 tagn2 = 0
262 tagn3 = 0
263 tagn4 = 0
264 fp = one
265 IF(ninterp > 0 ) THEN
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 . tagncont(npl,n3)==1) THEN
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 !N=1, NPRES/4
295 END DO !NP=1,NLOADP_HYD
296C---------pfluid
297 DO nl=1,nloadp_f
298C--------default zero
299 isiz_seg = iloadp(1,nl)/4
300 iad = iloadp(4,nl)
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
312 fun_hsp=iloadp(7,nl)
313 dir_hsp=iloadp(8,nl)
314 ifra1=iloadp(9,nl)
315 fcy = fac(1,nl)
316 fcx = fac(2,nl)
317 fun_cx=iloadp(10,nl)
318 fcy1 = fac(3,nl)
319 fcx1 = fac(4,nl)
320 fun_vel=iloadp(11,nl)
321 fcy2 = fac(5,nl)
322 fcx2 = fac(6,nl)
323 ! To avoid a check bound issue when the velocity options are not set in the input,
324 ! the DIR_VEL variable is bounded to a minimal value of 1
325 dir_vel=max(iloadp(12,nl),1)
326 ifra2=iloadp(13,nl)
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)
346C
347 aa = zero
348 vel = zero
349 pvel=zero
350C------ ----------
351C
352 IF(n4/=0 .AND. n1/=n2 .AND. n1/=n3 .AND. n1/=n4 .AND.
353 . n2/=n3 .AND. n2/=n4 .AND. n3/=n4 )THEN
354C
355 k1=3*dir_hsp-2
356 k2=3*dir_hsp-1
357 k3=3*dir_hsp
358 ! hydrostatic pressure
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
370C vel pressure
371 k1=3*dir_vel-2
372 k2=3*dir_vel-1
373 k3=3*dir_vel
374c
375 nsign = (nx * xframe(k1,ifra2) +
376 . ny * xframe(k2,ifra2) +
377 . nz * xframe(k3,ifra2))
378 IF(nsign/=zero) nsign = sign(one,nsign)
379C
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
397C
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
419C true triangles.
420 IF(fun_hsp /=0)THEN
421 k1=3*dir_hsp-2
422 k2=3*dir_hsp-1
423 k3=3*dir_hsp
424 ! hydrostatic pressure
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
435C vel pressure
436 k1=3*dir_vel-2
437 k2=3*dir_vel-1
438 k3=3*dir_vel
439c
440 nsign = (nx * xframe(k1,ifra2) +
441 . ny * xframe(k2,ifra2) +
442 . nz * xframe(k3,ifra2))
443 IF(nsign/=zero) nsign = sign(one,nsign)
444C
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
467C---------pblast
468 DO nl=1+nloadp_f,nloadp_f+pblast%NLOADP_B
469C--------default zero
470 isiz_seg = iloadp(1,nl)/4
471 iad = iloadp(4,nl)
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
483 il = nl-nloadp_f
484 tdet = fac(01,nl)
485 id = iloadp(08,nl) !user_id
486 !---------------------------------------------
487 ! LOOP ON SEGMENTS (4N or 3N)
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
503C---------/LOAD/PCYL
504 DO nl=1,loads%NLOAD_CYL
505C--------default zero
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 ! SEGP beam axis
541 dir(1) = dirx / len
542 dir(2) = diry / len
543 dir(3) = dirz / len
544 p0(1) = x(1,m2)
545 p0(2) = x(2,m2)
546 p0(3) = x(3,m2)
547 !---------------------------------------------
548 ! LOOP ON SEGMENTS (4N or 3N)
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 ! 3 node segment
567 CALL press_seg3(a ,b ,c ,p0 ,dir ,
568 . ifun ,table ,xfacr ,xfact ,segp )
569 press = abs(segp) * yfac
570c
571 ELSE ! 4 node segment
572 d(1) = x(1,n4)
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) + x(2,n4)) * fourth
577 m(3) = (x(3,n1) + x(3,n2) + x(3,n3) + x(3,n4)) * fourth
578c 1st internal triangle
579 CALL press_seg3(a ,b ,m ,p0 ,dir ,
580 . ifun ,table ,xfacr ,xfact ,segp )
581 press = press + segp * fourth
582c 2nd internal triangle
583 CALL press_seg3(b ,c ,m ,p0 ,dir ,
584 . ifun ,table ,xfacr ,xfact ,segp )
585 press = press + segp * fourth
586c 3rd internal triangle
587 CALL press_seg3(c ,d ,m ,p0 ,dir ,
588 . ifun ,table ,xfacr ,xfact ,segp )
589 press = press + segp * fourth
590c 4th internal triangle
591 CALL press_seg3(d ,a ,m ,p0 ,dir ,
592 . ifun ,table ,xfacr ,xfact ,segp )
593 press = press + segp * fourth
594 press = abs(press) * yfac
595 END IF ! seg 4 node
596 nskin = nskin0+ imapskp(np0)
597 skin_vector(3,nskin)=skin_vector(3,nskin) +press
598 enddo!next N
599 END DO
600C
601 RETURN
602 END
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
subroutine h3d_skin_vector(skin_vector, nodal_ipart, nsensor, is_written_skin, h3d_part, info1, keyword, ib, iloadp, lloadp, fac, npc, tf, sensor_tab, tagncont, loadp_hyd_inter, forc, xframe, x, v, imapskp, loads, table, iframe, dis, pblast)
#define max(a, b)
Definition macros.h:21
initmumps id
integer numskinp
Definition h3d_inc_mod.F:44
integer, parameter ncharline100
integer nploadpinch
subroutine press_seg3(a, b, c, n1, dir, ifunc, table, xfacr, xfact, press)
Definition press_seg3.F:37