OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
law70_upd.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!|| law70_upd ../starter/source/materials/mat/mat070/law70_upd.F
25!||--- called by ------------------------------------------------------
26!|| updmat ../starter/source/materials/updmat.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../starter/source/output/message/message.F
29!|| law70_table ../starter/source/materials/mat/mat070/law70_table.F
30!|| table_slope ../starter/source/materials/tools/table_slope.F
31!|| vw_smooth ../starter/source/tools/curve/vw_smooth.F
32!||--- uses -----------------------------------------------------
33!|| message_mod ../starter/share/message_module/message_mod.f
34!||====================================================================
35 SUBROUTINE law70_upd(MAT_PARAM,TITR ,MAT_ID ,NUPARAM ,UPARAM ,
36 . NFUNC ,IFUNC ,NPC ,PLD ,IOUT ,
37 . NFUNCT ,FUNC_ID ,NPROPM ,PM )
38C-----------------------------------------------
39C M o d u l e s
40C-----------------------------------------------
41 USE message_mod
42 USE matparam_def_mod
43 USE names_and_titles_mod , ONLY : nchartitle
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48C-----------------------------------------------
49C D u m m y A r g u m e n t s
50C-----------------------------------------------
51 INTEGER ,INTENT(IN) :: MAT_ID
52 INTEGER ,INTENT(IN) :: NFUNC
53 INTEGER ,INTENT(IN) :: NFUNCT
54 INTEGER ,INTENT(IN) :: NUPARAM
55 INTEGER ,INTENT(IN) :: IOUT
56 INTEGER ,INTENT(IN) :: NPROPM
57 INTEGER ,DIMENSION(NFUNC) ,INTENT(IN) :: IFUNC
58 INTEGER ,DIMENSION(NFUNCT) ,INTENT(IN) :: FUNC_ID
59 INTEGER ,INTENT(IN) :: NPC(*)
60 my_real ,INTENT(IN) :: pld(*)
61 CHARACTER(LEN=NCHARTITLE) ,INTENT(IN):: TITR
62 my_real ,DIMENSION(NPROPM) ,INTENT(INOUT) :: pm
63 my_real ,DIMENSION(NUPARAM) ,INTENT(INOUT) :: uparam
64 TYPE (MATPARAM_STRUCT_) ,INTENT(INOUT) :: MAT_PARAM
65C-----------------------------------------------
66C L o c a l V a r i a b l e s
67C-----------------------------------------------
68 INTEGER :: I,J,K,NDIM,NLOAD,NULOAD,NCURV,IPT,NPT,LMAX,FUNC_N,
69 . IC1,IC2,IZERO,IERROR,NTABLE,IS_ENCRYPTED,IFLAG0,IFLAG,ICHK
70 INTEGER ,PARAMETER :: LOAD = 1
71 integer ,PARAMETER :: unload = 2
72 INTEGER ,PARAMETER :: NPTMAX = 100 ! max number of function input points
73 my_real :: e0,emax,epsmax,stiffmin,stiffmax,stiffini,c1,g,nu,xmax,
74 . s1,s2,t1,t2,x1,x2,y1,y2,deri,yy,eps0,epst
75 INTEGER ,DIMENSION(NFUNC) :: FUNC,PERM,LEN
76 my_real ,DIMENSION(NFUNC*2) :: RATE,YFAC
77 my_real ,DIMENSION(:) ,ALLOCATABLE :: xf
78 my_real ,DIMENSION(:,:) ,ALLOCATABLE :: xi,yi,yf
79C=======================================================================
80 s2 = -huge(s2)
81 t2 = -huge(t2)
82 y1 = -huge(y1)
83 y2 = -huge(y2)
84 ncurv = int(uparam(1))
85 nload = int(uparam(7))
86 nuload = int(uparam(8))
87 e0 = uparam(2)
88 emax = uparam(2*nfunc + 12)
89 epsmax = uparam(4)
90 nu = uparam(6)
91 is_encrypted = nint(uparam(2*nfunc + 16))
92 IF (nuload > 0) THEN
93 ntable = 2
94 ELSE
95 ntable = 1
96 END IF
97 mat_param%NTABLE = ntable
98 ALLOCATE (mat_param%TABLE(ntable))
99 mat_param%TABLE(load)%NOTABLE = load
100 IF (ntable == 2) mat_param%TABLE(unload)%NOTABLE = unload
101c
102c------------------------------------------------------
103c Loading functions
104c------------------------------------------------------
105 lmax = 0
106 DO i = 1,nload
107 func_n = ifunc(i)
108 rate(i) = uparam(i + 8)
109 yfac(i) = uparam(i + 8 + nfunc)
110 len(i) = (npc(func_n+1) - npc(func_n)) / 2
111 lmax = max(lmax,len(i))
112 END DO
113 ALLOCATE (xi(lmax,nload))
114 ALLOCATE (yi(lmax,nload))
115c build X,Y vectors from loading functions
116c skip negative function values (only positive part is taken into account and symmetrized)
117c check and enforce passing through (0,0)
118 DO i = 1,nload
119 func_n = ifunc(i)
120 ic1 = npc(func_n)
121 ic2 = npc(func_n+1) - 2
122 ipt = 0
123 izero = 0
124 s1 = pld(ic1)
125 DO j = ic1,ic2-2,2
126 s1 = pld(j)
127 s2 = pld(j+2)
128 t1 = pld(j+1) * yfac(i)
129 t2 = pld(j+3) * yfac(i)
130 IF (s1 < zero .and. s2 < zero) cycle
131 IF (j == ic1 .and. s1 > zero) THEN
132 s1 = zero
133 t1 = zero
134 izero = 1
135 ELSE IF (s1 <= zero .and. s2 > zero) THEN
136 IF (t1 /= zero ) THEN
137 s1 = zero
138 t1 = zero
139 izero = 1
140 END IF
141 END IF
142 ipt = ipt + 1
143 xi(ipt,i) = s1
144 yi(ipt,i) = t1
145 END DO
146 ipt = ipt + 1
147 xi(ipt,i) = s2
148 yi(ipt,i) = t2
149 len(i) = ipt
150 END DO
151c------------------------------------------------------
152c interpolate if not passing through (0,0)
153c------------------------------------------------------
154 DO i = 1,nload
155 DO j = 1,len(i)
156 IF (xi(j,i) == zero .and. yi(j,i) /= zero) THEN
157 yi(j,i) = zero
158 ELSE IF (xi(j,i) == zero .and. xi(j,i) /= zero) THEN
159 xi(j,i) = zero
160 END IF
161 END DO
162 END DO
163c reduce number of points if necessary
164 DO i = 1,nload
165 IF (len(i) > nptmax) THEN
166 CALL vw_smooth(len(i),nptmax,xi(1:len(i),i),yi(1:len(i),i))
167 len(i) = nptmax
168 END IF
169 END DO
170c
171c------------------------------------------------------
172c
173 CALL law70_table(
174 . mat_param%TABLE(load) ,nload ,len ,lmax ,rate ,
175 . xi ,yi )
176
177 DEALLOCATE (yi)
178 DEALLOCATE (xi)
179c------------------------------------------------------
180c Unloading functions -> tables
181c------------------------------------------------------
182 IF (nuload > 0) THEN
183c
184 lmax = 0
185 DO i = 1,nuload
186 k = nload + i
187 func(i) = ifunc(k)
188 rate(i+nload) = uparam(k + 8)
189 yfac(i) = uparam(k + 8 + nfunc)
190 len(i) = (npc(func(i)+1) - npc(func(i))) / 2
191 lmax = max(lmax,len(i))
192 END DO
193 ALLOCATE (xi(lmax,nuload))
194 ALLOCATE (yi(lmax,nuload))
195c
196 DO i = 1,nuload
197 func_n = func(i)
198 ic1 = npc(func_n)
199 ic2 = npc(func_n+1) - 2
200 ipt = 0
201 izero = 0
202 s2 = zero
203 t2 = zero
204 DO j = ic1,ic2-2,2
205 s1 = pld(j)
206 s2 = pld(j+2)
207 t1 = pld(j+1) * yfac(i)
208 t2 = pld(j+3) * yfac(i)
209 IF (s1 < zero .and. s2 < zero) cycle
210 IF (j == ic1 .and. s1 > zero) THEN
211 s1 = zero
212 t1 = zero
213 izero = 1
214 ELSE IF (s1 <= zero .and. s2 > zero) THEN
215 IF (t1 /= zero ) THEN
216 s1 = zero
217 t1 = zero
218 izero = 1
219 END IF
220 END IF
221 ipt = ipt + 1
222 xi(ipt,i) = s1
223 yi(ipt,i) = t1
224 END DO
225 ipt = ipt + 1
226 xi(ipt,i) = s2
227 yi(ipt,i) = t2
228 len(i) = ipt
229 END DO
230c------------------------------------------------------
231c interpolate if not passing through (0,0)
232c------------------------------------------------------
233 DO i = 1,nuload
234 DO j = 1,len(i)
235 IF (xi(j,i) == zero .and. yi(j,i) /= zero) THEN
236 yi(j,i) = zero
237 ELSE IF (xi(j,i) == zero .and. xi(j,i) /= zero) THEN
238 xi(j,i) = zero
239 END IF
240 END DO
241 END DO
242c reduce number of points if necessary
243 DO i = 1,nuload
244 IF (len(i) > nptmax) THEN
245 CALL vw_smooth(len(i),nptmax,xi(1:len(i),i),yi(1:len(i),i))
246 len(i) = nptmax
247 END IF
248 END DO
249c--------------------------------------------------------
250 CALL law70_table(
251 . mat_param%TABLE(unload) ,nuload ,len ,lmax ,rate(nload+1) ,
252 . xi , yi )
253
254 DEALLOCATE (yi)
255 DEALLOCATE (xi)
256
257 END IF
258c--------------------------------------------------------
259c Automatic calculation of max function slope when Emax=0
260c--------------------------------------------------------
261 CALL table_slope(mat_param%TABLE(load),stiffini,stiffmin,stiffmax,xmax)
262c
263 IF (emax == zero) THEN
264 emax = stiffmax
265 uparam(3) = (emax - e0) / epsmax
266 uparam(2*nfunc + 12) = emax
267c
268 CALL ancmsg(msgid=1219, msgtype=msginfo, anmode=aninfo_blind_1,
269 . i1 = mat_id,
270 . c1 = titr ,
271 . r1 = emax )
272 END IF
273c
274 iflag = 0
275 iflag0 = 0
276 IF (e0 < stiffini) THEN
277 e0 = stiffini
278 IF (emax < e0) emax = e0
279 iflag0 = 1
280 END IF
281c--------------------------------------------------------
282c Automatic modification of EPST and E0
283c--------------------------------------------------------
284 eps0 = one
285 epst = one
286 x1 = mat_param%TABLE(load)%X(1)%VALUES(1)
287 ndim = mat_param%TABLE(load)%NDIM
288 npt = SIZE(mat_param%TABLE(load)%X(1)%VALUES)
289 DO k=1,nload
290 DO i = 1,npt-1
291 j = i+1
292 x2 = mat_param%TABLE(load)%X(1)%VALUES(j)
293 IF (ndim == 1) THEN
294 y1 = mat_param%TABLE(load)%Y1D(i)
295 y2 = mat_param%TABLE(load)%Y1D(j)
296 ELSE IF (ndim == 2) THEN
297 y1 = mat_param%TABLE(load)%Y2D(i,k)
298 y2 = mat_param%TABLE(load)%Y2D(j,k)
299 END IF
300 deri = (y2 - y1) / (x2 - x1)
301 IF (deri >= emax .and. x1 > zero) THEN
302 eps0 = min(eps0, x1 )
303 iflag = 1
304 IF (x1 == eps0) THEN
305 epst = min(epst,abs(eps0 - y1/emax))
306 ENDIF
307 ENDIF
308 x1 = x2
309 ENDDO
310 ENDDO ! NLOAD
311c
312 IF (iflag == 1) THEN
313 e0 = min(e0, emax)
314 uparam(3) = (emax - e0) / epst
315 uparam(4) = eps0
316 CALL ancmsg(msgid=864, msgtype=msginfo, anmode=aninfo_blind_1,
317 . i1=mat_id,
318 . c1=titr,
319 . r1=eps0)
320 ENDIF
321 IF (iflag0 == 1) THEN
322 e0 = min(e0, emax)
323 uparam(2) = e0
324 uparam(3) = (emax - e0)/epst
325 CALL ancmsg(msgid=865, msgtype=msgwarning, anmode=aninfo_blind_1,
326 . i1=mat_id,
327 . c1=titr,
328 . r1=e0)
329 ENDIF
330c------------------------------------
331c static function value at EPSMAX for engine
332c------------------------------------
333 k = 1
334 DO WHILE (mat_param%TABLE(load)%X(1)%VALUES(k) < epsmax .and. k < len(1)-1)
335 k = k + 1
336 IF(k >= len(1) -1 ) EXIT
337 END DO
338 x1 = mat_param%TABLE(load)%X(1)%VALUES(k-1)
339 x2 = mat_param%TABLE(load)%X(1)%VALUES(k)
340 y1 = -huge(y1)
341 y2 = -huge(y2)
342 IF (mat_param%TABLE(load)%NDIM == 1) THEN
343 y1 = mat_param%TABLE(load)%Y1D(k-1)
344 y2 = mat_param%TABLE(load)%Y1D(k)
345 ELSE IF (mat_param%TABLE(load)%NDIM == 2) THEN
346 y1 = mat_param%TABLE(load)%Y2D(k-1,1)
347 y2 = mat_param%TABLE(load)%Y2D(k,1)
348 END IF
349 deri = (y2 - y1) / (x2 - x1)
350 uparam(2*nfunc + 15) = y1 + deri * (epsmax - x1)
351c------------------------------------
352c update stiffness values in PM
353c------------------------------------
354 g = half *e0 / (one + nu)
355 c1 = third*e0 / (one - two*nu)
356 uparam(5) = g
357 pm(20) = e0
358 pm(22) = g
359 pm(24) = emax ! E0/(ONE - NU**2)
360 pm(32) = c1
361c------------------------------------
362c Output new tables definition
363c------------------------------------
364 IF (is_encrypted == 0)THEN
365 WRITE(iout,1000)
366 WRITE(iout,1001)
367 ndim = mat_param%TABLE(load)%NDIM
368 IF (ndim == 1) THEN
369 WRITE(iout,1101) func_id(ifunc(1))
370 DO j=1,SIZE(mat_param%TABLE(load)%X(1)%VALUES)
371 WRITE(iout,2000) mat_param%TABLE(load)%X(1)%VALUES(j),
372 . mat_param%TABLE(load)%Y1D(j)
373 END DO
374 ELSE
375 DO i=1,nload
376 WRITE(iout,1102) func_id(ifunc(i)),rate(i)
377 DO j=1,SIZE(mat_param%TABLE(load)%X(1)%VALUES)
378 WRITE(iout,2000) mat_param%TABLE(load)%X(1)%VALUES(j),
379 . mat_param%TABLE(load)%Y2D(j,i)
380 END DO
381 END DO
382 END IF
383c
384 IF (nuload == 1) THEN
385 WRITE(iout,1002)
386 WRITE(iout,1101) func_id(ifunc(1+nload))
387 DO j=1,SIZE(mat_param%TABLE(unload)%X(1)%VALUES)
388 WRITE(iout,2000) mat_param%TABLE(unload)%X(1)%VALUES(j),
389 . mat_param%TABLE(unload)%Y1D(j)
390 END DO
391 ELSE IF (nuload > 1) THEN
392 WRITE(iout,1002)
393 ndim = mat_param%TABLE(unload)%NDIM
394 IF (ndim == 1) THEN
395 WRITE(iout,1101) func_id(ifunc(1+nload))
396 DO j=1,SIZE(mat_param%TABLE(unload)%X(1)%VALUES)
397 WRITE(iout,2000) mat_param%TABLE(unload)%X(1)%VALUES(j),
398 . mat_param%TABLE(unload)%Y1D(j)
399 END DO
400 ELSE
401 DO i=1,nuload
402 WRITE(iout,1102) func_id(ifunc(i+nload)),rate(i+nload)
403 DO j=1,SIZE(mat_param%TABLE(unload)%X(1)%VALUES)
404 WRITE(iout,2000) mat_param%TABLE(unload)%X(1)%VALUES(j),
405 . mat_param%TABLE(unload)%Y2D(j,i)
406 END DO
407 END DO
408 END IF
409 END IF
410c
411 END IF
412c-----------
413 RETURN
414c-----------
415 1000 FORMAT(/,'------------------------------------------',/,
416 . 'MATERIAL LAW70 : UPDATE OF INPUT FUNCTIONS',/,
417 . '------------------------------------------',/)
418 1001 FORMAT(5x,'LOADING :')
419 1002 FORMAT(5x,'UNLOADING :')
420 1101 FORMAT(5x,/,'YIELD STRESS FUNCTION=',i10,
421 . 5x,/' X, Y')
422 1102 FORMAT(5x,/,'YIELD STRESS FUNCTION=',i10,
423 . 5x,'STRAIN RATE = ',1pg20.13,/,
424 . 5x,' X Y')
425 2000 FORMAT(2g20.13)
426c------------------------------------
427 END
my_real
#define my_real
Definition cppsort.cpp:32
law70_table
subroutine law70_table(table, nfunc, length, lmax, rate, xi, yi)
Definition law70_table.F:35
law70_upd
subroutine law70_upd(mat_param, titr, mat_id, nuparam, uparam, nfunc, ifunc, npc, pld, iout, nfunct, func_id, npropm, pm)
Definition law70_upd.F:38
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
message_mod
Definition message_mod.F:1249
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1003
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
starter
program starter
Definition starter.F:39
table_slope
subroutine table_slope(table, stiffini, stiffmin, stiffmax, xmax)
Definition table_slope.F:32
vw_smooth
subroutine vw_smooth(npt, ntarget, x, y)
Definition vw_smooth.F:30