OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
nrf51ini.F File Reference
#include "implicit_f.inc"
#include "param_c.inc"
#include "vect01_c.inc"
#include "com04_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine nrf51ini (ipm, pm, x, nix, ix, ale_connectivity, bufmat, uparam, rho0, uvar, nuvar, nel, rho, numel)

Function/Subroutine Documentation

◆ nrf51ini()

subroutine nrf51ini ( integer, dimension(npropmi,nummat), intent(in) ipm,
pm,
x,
integer, intent(in) nix,
integer, dimension(nix,numel), intent(in) ix,
type(t_ale_connectivity), intent(inout) ale_connectivity,
dimension(*), target bufmat,
uparam,
rho0,
dimension(nel,nuvar), intent(inout) uvar,
integer, intent(in) nuvar,
integer, intent(in) nel,
dimension(nel), intent(inout) rho,
integer, intent(in) numel )

Definition at line 32 of file nrf51ini.F.

36C-----------------------------------------------
37C D e s c r i p t i o n
38C-----------------------------------------------
39C This subroutines is automatically initializing
40C law51-NRF parameter according state from adjacent
41C element. User is no longer asked to copy/paste
42C material parameters. Starter is doing it automatically
43C
44C NRF READS
45C AV, RHO, E, PMIN, P0, SSP
46C 0.0 means automatic initialisation
47C otherwise use read
48C
49C NUMEL : total number of solid elements (2d or 3d)
50C NEL : Number of element in current group (<=MVSIZ)
51C IX : element connectivity + mat_id
52C IPM : material properties (integer)
53C PM : material properties (real)
54C
55C-----------------------------------------------
56C M o d u l e s
57C-----------------------------------------------
58 USE message_mod
59 USE ale_connectivity_mod
60 USE multimat_param_mod , ONLY : m51_n0phas, m51_nvphas, m51_tcp_ref, m51_lset_iflg6, m51_lc0max, m51_ssp0max, m51_iloop_nrf
61C-----------------------------------------------
62C I m p l i c i t T y p e s
63C-----------------------------------------------
64#include "implicit_f.inc"
65C-----------------------------------------------
66C C o m m o n B l o c k s
67C-----------------------------------------------
68#include "param_c.inc"
69#include "vect01_c.inc"
70#include "com04_c.inc"
71C-----------------------------------------------
72C D u m m y A r g u m e n t s
73C-----------------------------------------------
74 INTEGER,INTENT(IN) :: NUMEL,NIX,IX(NIX,NUMEL),IPM(NPROPMI,NUMMAT),NUVAR,NEL
75 my_real :: pm(npropm,nummat), x(3,numnod), uparam(*), rho0
76 my_real,TARGET :: bufmat(*)
77 my_real,INTENT(INOUT) :: uvar(nel,nuvar), rho(nel)
78 TYPE(t_ale_connectivity), INTENT(INOUT) :: ALE_CONNECTIVITY
79C-----------------------------------------------
80C L o c a l V a r i a b l e s
81C-----------------------------------------------
82 INTEGER I, IVMAT,IMAT,IID,MLN,J,J2,IADBUF,IFLGv,IE,IEV,NUPARAM,K,IFORM,IAD1,LGTH
83 my_real
84 . c01,c02,c03,c04, c01v,c02v,c03v,c04v,
85 . c11,c12,c13,c14, c11v,c12v,c13v,c14v,
86 . rho01,rho02,rho03,rho04, rho01v,rho02v,rho03v,rho04v,
87 . c21,c22,c23,c24, c21v,c22v,c23v,c24v,
88 . c31,c32,c33,c34, c31v,c32v,c33v,c34v,
89 . c41,c42,c43,c44, c41v,c42v,c43v,c44v,
90 . c51,c52,c53,c54, c51v,c52v,c53v,c54v,
91 . e01,e02,e03,e04, e01v,e02v,e03v,e04v,
92 . p01,p02,p03,p04, p01v,p02v,p03v,p04v,
93 . pm1,pm2,pm3,pm4, pm1v,pm2v,pm3v,pm4v,
94 . av1,av2,av3,av4, av1v,av2v,av3v,av4v,
95 . ssp1,ssp2,ssp3,ssp4, ssp1v,ssp2v,ssp3v,ssp4v,
96 . gg1v,gg2v,gg3v,
97 . dpdmu1v, dpdmu2v,dpdmu3v,
98 . pext, pextv, p0_nrf, p0_nrfv, avsum,
99 . lc,xmin,ymin,zmin,xmax,ymax,zmax, ssp0, tmp,
100 . av(4)
101
102 my_real,POINTER,DIMENSION(:) :: uprm
103
104 INTEGER lFOUND_51, lFOUND_other, lSET
105C-----------------------------------------------
106
107C=======================================================================
108C GET ADJACENT VALUES AND SET INITIAL STATE
109C=======================================================================
110 av1 = uparam(04)
111 av2 = uparam(05)
112 av3 = uparam(06)
113 av4 = uparam(46)
114 rho01 = uparam(09)
115 rho02 = uparam(10)
116 rho03 = uparam(11)
117 rho04 = uparam(47)
118 e01 = uparam(32)
119 e02 = uparam(33)
120 e03 = uparam(34)
121 e04 = uparam(48)
122 c11 = uparam(12)
123 c12 = uparam(13)
124 c13 = uparam(14)
125 c14 = uparam(50)
126 c21 = uparam(15)
127 c22 = uparam(16)
128 c23 = uparam(17)
129 c24 = zero
130 c31 = uparam(18)
131 c32 = uparam(19)
132 c33 = uparam(20)
133 c34 = zero
134 c41 = uparam(22)
135 c42 = uparam(23)
136 c43 = uparam(24)
137 c44 = zero
138 c51 = uparam(25)
139 c52 = uparam(26)
140 c53 = uparam(27)
141 c54 = zero
142 pm1 = uparam(39)
143 pm2 = uparam(40)
144 pm3 = uparam(41)
145 pm4 = uparam(56)
146 pext = uparam(8)
147 c01 = uparam(35) !here C01v is EOS coefficient and C01 is initial pressure
148 c02 = uparam(36)
149 c03 = uparam(37)
150 c04 = uparam(49)
151 ssp1 = uparam(124)
152 ssp2 = uparam(174)
153 ssp3 = uparam(224)
154 ssp4 = uparam(273)
155 p0_nrf = av1*(c01+c41*e01) + av2*(c02+c42*e02) + av3*(c03+c43*e03) + av4*(c04+c44*e04)
156 ssp0 = zero
157 p01 = c01+c41*e01
158 p02 = c02+c42*e02
159 p03 = c03+c43*e03
160 p04 = c04
161
162 DO i=lft,llt
163 ie = i + nft
164 iid = ix(nix,ie)
165 imat = ix(1,ie)
166 lfound_51 = 0
167 lfound_other = 0
168 nuparam = ipm(9,imat)
169 uvar(i,4) = p0_nrf ! default initialization for corner elems
170 iad1 = ale_connectivity%ee_connect%iad_connect(ie)
171 lgth = ale_connectivity%ee_connect%iad_connect(ie+1) - iad1
172 ivmat = 0
173 DO j=1,lgth
174 iev = ale_connectivity%ee_connect%connected(iad1 + j - 1)
175 IF(iev == 0)cycle
176 ivmat = abs(ix(1,iev))
177 mln = nint(pm(19,ivmat))
178 iflgv = 0
179 IF(mln/=51)lfound_other = lfound_other + 1
180 iadbuf = ipm(7,ivmat)
181 iflgv = nint(bufmat(iadbuf-1+31))
182 IF(iflgv>=2 .AND. iflgv<=6 )cycle
183 lfound_51 = lfound_51 + 1
184 !-----------------------------------------!
185 ! GET ADJACENT LAW51 PARAMETER !
186 ! necessarily defined - check by lecm51 !
187 !-----------------------------------------!
188 av1v = bufmat(iadbuf-1+04) ;
189 av2v = bufmat(iadbuf-1+05) ;
190 av3v = bufmat(iadbuf-1+06) ;
191 av4v = bufmat(iadbuf-1+46) ;
192 rho01v = bufmat(iadbuf-1+09) ;
193 rho02v = bufmat(iadbuf-1+10) ;
194 rho03v = bufmat(iadbuf-1+11) ;
195 rho04v = bufmat(iadbuf-1+47) ;
196 e01v = bufmat(iadbuf-1+32) ;
197 e02v = bufmat(iadbuf-1+33) ;
198 e03v = bufmat(iadbuf-1+34) ;
199 e04v = bufmat(iadbuf-1+48) ;
200 c11v = bufmat(iadbuf-1+12) ;
201 c12v = bufmat(iadbuf-1+13) ;
202 c13v = bufmat(iadbuf-1+14) ;
203 c14v = bufmat(iadbuf-1+50) ;
204 c21v = bufmat(iadbuf-1+15) ;
205 c22v = bufmat(iadbuf-1+16) ;
206 c23v = bufmat(iadbuf-1+17) ;
207 c24v = zero ;
208 c31v = bufmat(iadbuf-1+18) ;
209 c32v = bufmat(iadbuf-1+19) ;
210 c33v = bufmat(iadbuf-1+20) ;
211 c34v = zero ;
212 c41v = bufmat(iadbuf-1+22) ;
213 c42v = bufmat(iadbuf-1+23) ;
214 c43v = bufmat(iadbuf-1+24) ;
215 c44v = zero ;
216 c51v = bufmat(iadbuf-1+25) ;
217 c52v = bufmat(iadbuf-1+26) ;
218 c53v = bufmat(iadbuf-1+27) ;
219 c54v = zero ;
220 pm1v = bufmat(iadbuf-1+39) ;
221 pm2v = bufmat(iadbuf-1+40) ;
222 pm3v = bufmat(iadbuf-1+41) ;
223 pm4v = bufmat(iadbuf-1+56) ;
224 pextv = bufmat(iadbuf-1+8) ;
225 c01v = bufmat(iadbuf-1+35) ;
226 c02v = bufmat(iadbuf-1+36) ;
227 c03v = bufmat(iadbuf-1+37) ;
228 c04v = bufmat(iadbuf-1+49) ;
229 p01v = c01v + c41v*e01v
230 p02v = c02v + c42v*e02v
231 p03v = c03v + c43v*e03v
232 p04v = c04v
233 ssp1v = zero
234 ssp2v = zero
235 ssp3v = zero
236 ssp4v = zero
237 dpdmu1v = (c11v+c51v*e01v) + c41v*(pextv+p01v)
238 dpdmu2v = (c12v+c52v*e02v) + c42v*(pextv+p02v)
239 dpdmu3v = (c13v+c53v*e03v) + c43v*(pextv+p03v)
240 gg1v = bufmat(iadbuf-1+28)
241 gg2v = bufmat(iadbuf-1+29)
242 gg3v = bufmat(iadbuf-1+30)
243 IF(rho01v /= zero) ssp1v = sqrt( (dpdmu1v + two_third*gg1v) / rho01v ) !warning GG1v = 2*G1v => TWO_THIRD*GG = FOUR_OVER_3*G
244 IF(rho02v /= zero) ssp2v = sqrt( (dpdmu2v + two_third*gg2v) / rho02v )
245 IF(rho03v /= zero) ssp3v = sqrt( (dpdmu3v + two_third*gg3v) / rho03v )
246 ssp4v = bufmat(iadbuf-1+42) !VDETv
247 !-----------------------------------------!
248 ! SET AUTOMATIC LAW51 PARAMETERS !
249 !-----------------------------------------!
250 uvar(i,4) = p0_nrf
251 p0_nrfv = av1v*p01v + av2v*p02v + av3v*p03v + av4v*p04v
252 IF(p0_nrf==zero)uvar(i,4) = p0_nrfv
253
254 avsum = av1+av2+av3+av4
255 IF(avsum==zero) THEN
256 av(1) = av1v
257 av(2) = av2v
258 av(3) = av3v
259 av(4) = av4v
260 !
261 k=m51_n0phas+(1-1)*m51_nvphas
262 uvar(i,k+01)=av1v
263 uvar(i,k+23)=av1v
264 k=m51_n0phas+(2-1)*m51_nvphas
265 uvar(i,k+01)=av2v
266 uvar(i,k+23)=av2v
267 k=m51_n0phas+(3-1)*m51_nvphas
268 uvar(i,k+01)=av3v
269 uvar(i,k+23)=av3v
270 k=m51_n0phas+(4-1)*m51_nvphas
271 uvar(i,k+01)=av4v
272 uvar(i,k+23)=av4v
273 ELSE
274 av(1) = av1
275 av(2) = av2
276 av(3) = av3
277 av(4) = av4
278 k=m51_n0phas+(1-1)*m51_nvphas
279 uvar(i,k+01)=av1
280 uvar(i,k+23)=av1
281 k=m51_n0phas+(2-1)*m51_nvphas
282 uvar(i,k+01)=av2
283 uvar(i,k+23)=av2
284 k=m51_n0phas+(3-1)*m51_nvphas
285 uvar(i,k+01)=av3
286 uvar(i,k+23)=av3
287 k=m51_n0phas+(4-1)*m51_nvphas
288 uvar(i,k+01)=av4
289 uvar(i,k+23)=av4
290 ENDIF
291
292 k=m51_n0phas+(1-1)*m51_nvphas
293 IF(rho01==zero)THEN
294 uvar(i,k+09)= rho01v
295 uvar(i,k+12)= rho01v
296 uvar(i,k+20)= rho01v
297 rho01 =rho01v
298 ELSE
299 uvar(i,k+09)= rho01
300 uvar(i,k+12)= rho01
301 uvar(i,k+20)= rho01
302 ENDIF
303 IF(e01==zero)THEN
304 uvar(i,k+08)= e01v
305 uvar(i,k+21)= e01v
306 ELSE
307 uvar(i,k+08)= e01
308 uvar(i,k+21)= e01
309 ENDIF
310 IF(c01==zero)THEN
311 uvar(i,k+18)= p01v
312 ELSE
313 uvar(i,k+18)= p01
314 ENDIF
315 IF(ssp1==zero)THEN
316 uvar(i,k+22)= ssp1v
317 ssp1 = ssp1v
318 ssp0 = max(ssp0,ssp1)
319 ELSE
320 uvar(i,k+22)= ssp1
321 ssp0 = max(ssp0,ssp1)
322 ENDIF
323
324 k=m51_n0phas+(2-1)*m51_nvphas
325 IF(rho02==zero)THEN
326 uvar(i,k+09)= rho02v
327 uvar(i,k+12)= rho02v
328 uvar(i,k+20)= rho02v
329 rho02 =rho02v
330 ELSE
331 uvar(i,k+09)= rho02
332 uvar(i,k+12)= rho02
333 uvar(i,k+20)= rho02
334 ENDIF
335 IF(e02==zero)THEN
336 uvar(i,k+08)= e02v
337 uvar(i,k+21)= e02v
338 ELSE
339 uvar(i,k+08)= e02
340 uvar(i,k+21)= e02
341 ENDIF
342 IF(c02==zero)THEN
343 uvar(i,k+18)= p02v
344 ELSE
345 uvar(i,k+18)= p02
346 ENDIF
347 IF(ssp2==zero)THEN
348 uvar(i,k+22)= ssp2v
349 ssp2 = ssp2v
350 ssp0 = max(ssp0,ssp2)
351 ELSE
352 uvar(i,k+22)= ssp2
353 ssp0 = max(ssp0,ssp2)
354 ENDIF
355
356 k=m51_n0phas+(3-1)*m51_nvphas
357 IF(rho03==zero)THEN
358 uvar(i,k+09)= rho03v
359 uvar(i,k+12)= rho03v
360 uvar(i,k+20)= rho03v
361 rho03 =rho03v
362 ELSE
363 uvar(i,k+09)= rho03
364 uvar(i,k+12)= rho03
365 uvar(i,k+20)= rho03
366 ENDIF
367 IF(e03==zero)THEN
368 uvar(i,k+08)= e03v
369 uvar(i,k+21)= e03v
370 ELSE
371 uvar(i,k+08)= e03
372 uvar(i,k+21)= e03
373 ENDIF
374 IF(c03==zero)THEN
375 uvar(i,k+18)= p03v
376 ELSE
377 uvar(i,k+18)= p03
378 ENDIF
379 IF(ssp3==zero)THEN
380 uvar(i,k+22)= ssp3v
381 ssp3 = ssp3v
382 ssp0 = max(ssp0,ssp3)
383 ELSE
384 uvar(i,k+22)= ssp3
385 ssp0 = max(ssp0,ssp3)
386 ENDIF
387
388 k=m51_n0phas+(4-1)*m51_nvphas
389 IF(rho04==zero)THEN
390 uvar(i,k+09)= rho04v
391 uvar(i,k+12)= rho04v
392 uvar(i,k+20)= rho04v
393 rho04 =rho04v
394 ELSE
395 uvar(i,k+09)= rho04
396 uvar(i,k+12)= rho04
397 uvar(i,k+20)= rho04
398 ENDIF
399 IF(e04==zero)THEN
400 uvar(i,k+08)= e04v
401 uvar(i,k+21)= e04v
402 ELSE
403 uvar(i,k+08)= e04
404 uvar(i,k+21)= e04
405 ENDIF
406 IF(c04==zero)THEN
407 uvar(i,k+18)= p04v
408 ELSE
409 uvar(i,k+18)= p04
410 ENDIF
411 IF(ssp4==zero)THEN
412 uvar(i,k+22)= ssp4v
413 ssp4 = ssp4v
414 ssp0 = max(ssp0,ssp4)
415 ELSE
416 uvar(i,k+22)= ssp4
417 ssp0 = max(ssp0,ssp4)
418 ENDIF
419
420 rho0 = av(1)*rho01+av(2)*rho02+av(3)*rho03+av(4)*rho04
421 pm( 1,imat) = rho0 !NINT(PM(1,IVMAT))
422 pm(89,imat) = rho0 !NINT(PM(19,IVMAT))
423 rho(i) = rho0
424
425 EXIT
426 enddo! next J
427 !checking input file
428 !2 adjacent elem in computation domain: this is not expected. 1 boundary elem => 1 adjacent elem
429 DO j2=j+1,lgth
430 iev = ale_connectivity%ee_connect%connected(iad1 + j2 - 1)
431 IF(iev /= 0)THEN
432 mln = nint(pm(19,ivmat))
433 iflgv = 0
434 ivmat = abs(ix(1,iev))
435 IF(mln/=51 )lfound_other = lfound_other + 1
436 iadbuf = ipm(7,ivmat)
437 iflgv = nint(bufmat(iadbuf-1+31))
438 IF(iflgv>=2 .AND. iflgv<=6 )cycle
439 lfound_51 = lfound_51 + 1
440 ENDIF
441 ENDDO
442 IF(lfound_other /= 0)THEN
443 CALL ancmsg(msgid = 122 ,
444 . msgtype = msgerror,
445 . anmode = aninfo ,
446 . i1 = iid ,
447 . c1 = "MUST BE ADJACENT TO MM-ALE LAW51 PART" )
448 ENDIF
449 IF(lfound_51 == 0)THEN
450 CALL ancmsg(msgid = 123 ,
451 . msgtype = msgwarning,
452 . anmode = aninfo ,
453 . i1 = iid ,
454 . c1 = "HAS NO ADJACENT FACE IN COMPUTATION DOMAIN" )
455 ENDIF
456 IF(lfound_51 >= 2)THEN
457 CALL ancmsg(msgid = 122 ,
458 . msgtype = msgerror,
459 . anmode = aninfo ,
460 . i1 = iid ,
461 . c1 = "MUST HAVE ONLY ONE FACE ADJACENT TO COMPUTATION DOMAIN" )
462 ENDIF
463 ENDDO !next I
464
465 m51_ssp0max=max(m51_ssp0max,ssp0)
466
467C=======================================================================
468C AUTOMATIC CHARACTERISTIC LENGTH COMPUTATION
469C=======================================================================
470
471 !-----------------------------!
472 ! SSP0MAX (sound speed) !
473 !-----------------------------!
474 IF(m51_iloop_nrf==0)THEN
475 DO i=1,numels+numelq
476 imat = ix(1,i)
477 mln = ipm(2,imat)
478 iadbuf = ipm(7,imat)
479 IF(mln/=51)cycle
480 uprm => bufmat(iadbuf:iadbuf+275)
481 tmp = uprm(31)
482 iform = nint(tmp)
483 IF(iform>1)cycle
484 ssp0 = uprm(124)
485 ssp0 = max(ssp0,uprm(174))
486 ssp0 = max(ssp0,uprm(224))
487 ssp0 = max(ssp0,uprm(273))
488 ENDDO
489 m51_ssp0max = ssp0
490 m51_iloop_nrf = 1
491 ENDIF
492 !-----------------------------!
493 ! Characteristic Length LC !
494 !-----------------------------!
495 IF(m51_lc0max==zero)THEN
496 xmin = x(1,1)
497 ymin = x(2,1)
498 zmin = x(3,1)
499 xmax = x(1,1)
500 ymax = x(2,1)
501 zmax = x(3,1)
502 DO i=1,numnod
503 xmin = min(xmin,x(1,i))
504 ymin = min(ymin,x(2,i))
505 zmin = min(zmin,x(3,i))
506 xmax = max(xmax,x(1,i))
507 ymax = max(ymax,x(2,i))
508 zmax = max(zmax,x(3,i))
509 ENDDO
510 lc = xmax-xmin
511 lc = max(lc,ymax-ymin)
512 lc = max(lc,zmax-zmin)
513 m51_lc0max = lc
514 ENDIF
515 !-----------------------------!
516 ! Tcp : characteric time !
517 !-----------------------------!
518 m51_tcp_ref = m51_lc0max/two/m51_ssp0max/log(two)
519 !Si Tcp=0 alors calcul automatique
520 IF(uparam(70)==zero)THEN
521 uparam(70) = m51_tcp_ref
522 lset = 1
523 ELSE
524 lset = 0
525 ENDIF
526
527 if(lset==1) m51_lset_iflg6 = 1
528C-----------------------------------------------
529 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
ymax
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
Definition law100_upd.F:272
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
ale_connectivity_mod
Definition ale_connectivity_mod.F:223
message_mod
Definition message_mod.F:1249
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889