OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat112_xia_newton.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!|| mat112_xia_newton ../engine/source/materials/mat/mat112/mat112_xia_newton.F
25!||--- called by ------------------------------------------------------
26!|| sigeps112 ../engine/source/materials/mat/mat112/sigeps112.F
27!||--- calls -----------------------------------------------------
28!|| table2d_vinterp_log ../engine/source/tools/curve/table2d_vinterp_log.f
29!||--- uses -----------------------------------------------------
30!|| interface_table_mod ../engine/share/modules/table_mod.F
31!|| table_mod ../engine/share/modules/table_mod.F
32!||====================================================================
33 SUBROUTINE mat112_xia_newton(
34 1 NEL ,NGL ,NUPARAM ,NUVAR ,TIME ,TIMESTEP,
35 2 UPARAM ,UVAR ,JTHE ,OFF ,RHO0 ,RHO ,
36 3 PLA ,DPLA ,EPSD ,SOUNDSP ,EPSZZ ,
37 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 5 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
39 6 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
40 7 SIGY ,ET ,
41 8 NVARTMP ,NUMTABL ,VARTMP ,ITABLE ,TABLE )
42 !=======================================================================
43 ! Modules
44 !=======================================================================
45 USE table_mod
46 USE interface_table_mod
47 !=======================================================================
48 ! Implicit types
49 !=======================================================================
50#include "implicit_f.inc"
51 !=======================================================================
52 ! Common
53 !=======================================================================
54#include "com04_c.inc"
55 !=======================================================================
56 ! Dummy arguments
57 !=======================================================================
58 INTEGER NEL,NUPARAM,NUVAR,JTHE,NUMTABL,ITABLE(NUMTABL),NVARTMP
59 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
60 my_real
61 . TIME,TIMESTEP
62 INTEGER :: VARTMP(NEL,NVARTMP)
63 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
64 . UPARAM
65 my_real,DIMENSION(NEL), INTENT(IN) ::
66 . rho,rho0,epszz,
67 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
68 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
69c
70 my_real ,DIMENSION(NEL), INTENT(OUT) ::
71 . soundsp,sigy,et,
72 . signxx,signyy,signzz,signxy,signyz,signzx
73c
74 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
75 . dpla,off
76 my_real ,DIMENSION(NEL,4),INTENT(INOUT) ::
77 . pla,epsd
78 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
79 . uvar
80c
81 TYPE(ttable), DIMENSION(NTABLE) :: TABLE
82 !=======================================================================
83 ! Local Variables
84 !=======================================================================
85 INTEGER I,K,II,NINDX,INDEX(NEL),NINDX2,INDEX2(NEL),NINDX3,INDEX3(NEL),
86 . niter,iter,itab,ismooth,ipos(nel,2)
87c
88 my_real
89 . young1,young2,young3,nu12,nu21,
90 . a11,a12,a21,a22,
91 . g12,g23,g31,deuxk,e3c,cc,c1,
92 . nu1p,nu2p,nu4p,nu5p,s01,a01,b01,c01,
93 . s02,a02,b02,c02,s03,a03,b03,c03,
94 . s04,a04,b04,c04,s05,a05,b05,c05,
95 . asig,bsig,csig,tau0,atau,btau,n3,n6,
96 . bulk,nu13,nu31,nu23,nu32
97 my_real
98 . normsig,h,dpdt,dlam,ddep,depxx,depyy
99 my_real
100 . normxx,normyy,normxy,normpxx,normpyy,normpxy,
101 . dphi_dlam,dphi,dfdsig2,dphi_dpla,dpxx,dpyy,dpxy,
102 . dphi_dsigy1,dphi_dsigy2,dphi_dsigy3,dphi_dsigy4,
103 . dphi_dsigy5
104 my_real
105 . normzz,dtheta_dlam,dtheta,dtheta_dpla,
106 . dtheta_dsigyo
107 my_real
108 . normyz,normzx,dpsi_dlam,dpsi,dpsi_dpla,
109 . dpyz,dpzx,dpsi_dsigys
110 my_real
111 . xscale1,yscale1,xscale2,yscale2,xscale3,yscale3,xscale4,yscale4,
112 . xscale5,yscale5,xscalec,yscalec,xscales,yscales,xvec(nel,2),asrate
113c
114 my_real, DIMENSION(2) ::
115 . n1,n2,n4,n5
116c
117 my_real, DIMENSION(NEL) ::
118 . khi1,khi2,khi3,khi4,khi5,khi6,sigy1,sigy2,sigy3,sigy4,sigy5,dpla2,dpla3,
119 . dpla4,sigys,sigyo,dsigxx,dsigyy,dsigxy,dsigyz,dsigzx,hardp,phi,phi0,psi,psi0,
120 . theta,theta0,dsigzz,eezz,gam1,gam2,gam3,gam4,gam5,gam6,dsigy1_dp,dsigy2_dp,
121 . dsigy3_dp,dsigy4_dp,dsigy5_dp,dsigyo_dp,dsigys_dp,hardr
122c
123 !=======================================================================
124 ! - INITIALISATION OF COMPUTATION ON TIME STEP
125 !=======================================================================
126 ! Recovering model parameters
127 ! Elastic parameters
128 young1 = uparam(1) ! Young modulus in direction 1 (MD)
129 young2 = uparam(2) ! Young modulus in direction 2 (CD)
130 young3 = uparam(3) ! Young modulus in direction 3 (Out of plane)
131 nu12 = uparam(4) ! Poisson's ratio in 12
132 nu21 = uparam(5) ! Poisson's ratio in 21
133 a11 = uparam(6) ! Component 11 of orthotropic shell elasticity matrix
134 a12 = uparam(7) ! Component 12 of orthotropic shell elasticity matrix
135 a21 = uparam(8) ! Component 21 of orthotropic shell elasticity matrix
136 a22 = uparam(9) ! Component 22 of orthotropic shell elasticity matrix
137 g12 = uparam(10) ! Shear modulus in 12
138 g23 = uparam(11) ! Shear modulus in 23
139 g31 = uparam(12) ! Shear modulus in 31
140 itab = int(uparam(14)) ! Tabulated yield stress flag
141 deuxk = uparam(15) ! Yield criterion exponent
142 e3c = uparam(16) ! Compressive elasticity parameter
143 cc = uparam(17) ! Exponent compressive elasticity parameter
144 nu1p = uparam(18) ! Tensile plastic Poisson ratio in direction 1 (MD)
145 nu2p = uparam(19) ! Tensile plastic Poisson ratio in direction 2 (CD)
146 nu4p = uparam(20) ! Compressive plastic Poisson ratio in direction 1 (MD)
147 nu5p = uparam(21) ! Compressive plastic Poisson ratio in direction 2 (CD)
148 IF (itab == 0) THEN
149 s01 = uparam(22) ! 1st Tensile plasticity parameter direction 1 (MD)
150 a01 = uparam(23) ! 2nd Tensile plasticity parameter direction 1 (MD)
151 b01 = uparam(24) ! 3rd Tensile plasticity parameter direction 1 (MD)
152 c01 = uparam(25) ! 4th Tensile plasticity parameter direction 1 (MD)
153 s02 = uparam(26) ! 1st Tensile plasticity parameter direction 2 (CD)
154 a02 = uparam(27) ! 2nd Tensile plasticity parameter direction 2 (CD)
155 b02 = uparam(28) ! 3rd Tensile plasticity parameter direction 2 (CD)
156 c02 = uparam(29) ! 4th Tensile plasticity parameter direction 2 (CD)
157 s03 = uparam(30) ! 1st Positive shear plasticity parameter
158 a03 = uparam(31) ! 2nd Positive shear plasticity parameter
159 b03 = uparam(32) ! 3rd Positive shear plasticity parameter
160 c03 = uparam(33) ! 4th Positive shear plasticity parameter
161 s04 = uparam(34) ! 1st Compressive plasticity parameter direction 1 (MD)
162 a04 = uparam(35) ! 2nd Compressive plasticity parameter direction 1 (MD)
163 b04 = uparam(36) ! 3rd Compressive plasticity parameter direction 1 (MD)
164 c04 = uparam(37) ! 4th Compressive plasticity parameter direction 1 (MD)
165 s05 = uparam(38) ! 1st Compressive plasticity parameter direction 2 (CD)
166 a05 = uparam(39) ! 2nd Compressive plasticity parameter direction 2 (CD)
167 b05 = uparam(40) ! 3rd Compressive plasticity parameter direction 2 (CD)
168 c05 = uparam(41) ! 4th Compressive plasticity parameter direction 2 (CD)
169 asig = uparam(42) ! 1st Out-of-plane plasticity parameter
170 bsig = uparam(43) ! 2nd Out-of-plane plasticity parameter
171 csig = uparam(44) ! 3rd Out-of-plane plasticity parameter
172 tau0 = uparam(45) ! 1st Transverse shear plasticity parameter
173 atau = uparam(46) ! 2nd Transverse shear plasticity parameter
174 btau = uparam(47) ! 3rd Transverse shear plasticity parameter
175 ELSE
176 xscale1 = uparam(22)
177 yscale1 = uparam(23)
178 xscale2 = uparam(24)
179 yscale2 = uparam(25)
180 xscale3 = uparam(26)
181 yscale3 = uparam(27)
182 xscale4 = uparam(28)
183 yscale4 = uparam(29)
184 xscale5 = uparam(30)
185 yscale5 = uparam(31)
186 xscalec = uparam(32)
187 yscalec = uparam(33)
188 xscales = uparam(34)
189 yscales = uparam(35)
190 asrate = uparam(36)
191 asrate = (two*pi*asrate*timestep)/(two*pi*asrate*timestep + one)
192 ismooth = int(uparam(37))
193 ENDIF
194c
195 ! Yield planes normal
196 n1(1) = one/(sqrt(one+nu1p**2))
197 n1(2) = -nu1p/(sqrt(one+nu1p**2))
198 n2(1) = -nu2p/(sqrt(one+nu2p**2))
199 n2(2) = one/(sqrt(one+nu2p**2))
200 n3 = one
201 n4(1) = -one/(sqrt(one+nu4p**2))
202 n4(2) = nu4p/(sqrt(one+nu4p**2))
203 n5(1) = nu5p/(sqrt(one+nu5p**2))
204 n5(2) = -one/(sqrt(one+nu5p**2))
205 n6 = -one
206c
207 ! Recovering internal variables
208 DO i=1,nel
209 ! OFF parameter for element deletion
210 IF (off(i) < em01) off(i) = zero
211 IF (off(i) < one) off(i) = off(i)*four_over_5
212 ! User inputs
213 eezz(i) = epszz(i) + pla(i,3) ! out-of-plane elastic strain
214 ! Standard inputs
215 dpla(i) = zero ! Initialization of the "global" plastic strain increment
216 dpla2(i) = zero ! Initialization of the in-plane plastic strain increment
217 dpla3(i) = zero ! Initialization of the out-of-plane plastic strain increment
218 dpla4(i) = zero ! Initialization of the transverse shear plastic strain increment
219 et(i) = one ! Initialization of hourglass coefficient
220 hardp(i) = zero ! Initialization of hourglass coefficient
221 ENDDO
222c
223 ! Computation of the initial yield stress
224 IF (itab > 0) THEN
225 ! In-plane tabulation with strain-rate
226 xvec(1:nel,1) = pla(1:nel,2)
227 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
228 ! -> Tensile yield stress in direction 1 (MD)
229 ipos(1:nel,1) = vartmp(1:nel,1)
230 ipos(1:nel,2) = 1
231 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,sigy1,dsigy1_dp,hardr)
232 sigy1(1:nel) = sigy1(1:nel) * yscale1
233 dsigy1_dp(1:nel)= dsigy1_dp(1:nel) * yscale1
234 vartmp(1:nel,1) = ipos(1:nel,1)
235 ! -> Tensile yield stress in direction 2 (CD)
236 xvec(1:nel,2) = epsd(1:nel,2) * xscale2
237 ipos(1:nel,1) = vartmp(1:nel,2)
238 ipos(1:nel,2) = 1
239 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,sigy2,dsigy2_dp,hardr)
240 sigy2(1:nel) = sigy2(1:nel) * yscale2
241 dsigy2_dp(1:nel)= dsigy2_dp(1:nel) * yscale2
242 vartmp(1:nel,2) = ipos(1:nel,1)
243 ! -> Positive shear yield stress
244 xvec(1:nel,2) = epsd(1:nel,2) * xscale3
245 ipos(1:nel,1) = vartmp(1:nel,3)
246 ipos(1:nel,2) = 1
247 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,sigy3,dsigy3_dp,hardr)
248 sigy3(1:nel) = sigy3(1:nel) * yscale3
249 dsigy3_dp(1:nel)= dsigy3_dp(1:nel) * yscale3
250 vartmp(1:nel,3) = ipos(1:nel,1)
251 ! -> Compressive yield stress in direction 1 (MD)
252 xvec(1:nel,2) = epsd(1:nel,2) * xscale4
253 ipos(1:nel,1) = vartmp(1:nel,4)
254 ipos(1:nel,2) = 1
255 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,sigy4,dsigy4_dp,hardr)
256 sigy4(1:nel) = sigy4(1:nel) * yscale4
257 dsigy4_dp(1:nel)= dsigy4_dp(1:nel) * yscale4
258 vartmp(1:nel,4) = ipos(1:nel,1)
259 ! -> Compressive yield stress in direction 2 (CD)
260 xvec(1:nel,2) = epsd(1:nel,2) * xscale5
261 ipos(1:nel,1) = vartmp(1:nel,5)
262 ipos(1:nel,2) = 1
263 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,sigy5,dsigy5_dp,hardr)
264 sigy5(1:nel) = sigy5(1:nel) * yscale5
265 dsigy5_dp(1:nel)= dsigy5_dp(1:nel) * yscale5
266 vartmp(1:nel,5) = ipos(1:nel,1)
267 ! -> Out-of-plane tabulation with strain-rate
268 xvec(1:nel,1) = pla(1:nel,3)
269 xvec(1:nel,2) = epsd(1:nel,3) * xscalec
270 ! -> Transverse shear yield stress
271 ipos(1:nel,1) = vartmp(1:nel,6)
272 ipos(1:nel,2) = 1
273 CALL table2d_vinterp_log(table(itable(6)),ismooth,nel,nel,ipos,xvec,sigyo,dsigyo_dp,hardr)
274 sigyo(1:nel) = sigyo(1:nel) * yscalec
275 dsigyo_dp(1:nel)= dsigyo_dp(1:nel) * yscalec
276 vartmp(1:nel,6) = ipos(1:nel,1)
277 ! -> Transverse shear tabulation with strain-rate
278 xvec(1:nel,1) = pla(1:nel,4)
279 xvec(1:nel,2) = epsd(1:nel,4) * xscales
280 ! -> Transverse shear yield stress
281 ipos(1:nel,1) = vartmp(1:nel,7)
282 ipos(1:nel,2) = 1
283 CALL table2d_vinterp_log(table(itable(7)),ismooth,nel,nel,ipos,xvec,sigys,dsigys_dp,hardr)
284 sigys(1:nel) = sigys(1:nel) * yscales
285 dsigys_dp(1:nel)= dsigys_dp(1:nel) * yscales
286 vartmp(1:nel,7) = ipos(1:nel,1)
287 ELSE
288 DO i = 1,nel
289 ! Tensile yield stress in direction 1 (MD)
290 sigy1(i) = s01 + a01*tanh(b01*pla(i,2)) + c01*pla(i,2)
291 ! Tensile yield stress in direction 2 (CD)
292 sigy2(i) = s02 + a02*tanh(b02*pla(i,2)) + c02*pla(i,2)
293 ! Positive shear yield stress
294 sigy3(i) = s03 + a03*tanh(b03*pla(i,2)) + c03*pla(i,2)
295 ! Compressive yield stress in direction 1 (MD)
296 sigy4(i) = s04 + a04*tanh(b04*pla(i,2)) + c04*pla(i,2)
297 ! Compressive yield stress in direction 2 (CD)
298 sigy5(i) = s05 + a05*tanh(b05*pla(i,2)) + c05*pla(i,2)
299 ! Out-of-plane yield stress
300 sigyo(i) = asig + bsig*exp(csig*pla(i,3))
301 ! Transverse shear yield stress
302 sigys(i) = tau0 + (atau - min(zero,sigozz(i))*btau)*pla(i,4)
303 ENDDO
304 ENDIF
305c
306 !========================================================================
307 ! - COMPUTATION OF TRIAL VALUES
308 !========================================================================
309 DO i=1,nel
310c
311 ! Computation of the trial stress tensor
312 ! - in-plane stresses
313 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
314 signyy(i) = sigoyy(i) + a21*depsxx(i) + a22*depsyy(i)
315 signxy(i) = sigoxy(i) + depsxy(i)*g12
316 ! - out-of-plane stress
317 IF (eezz(i) >= zero) THEN
318 signzz(i) = sigozz(i) + young3*depszz(i)
319 ELSE
320 signzz(i) = sigozz(i) + e3c*cc*exp(-cc*eezz(i))*depszz(i)
321 ENDIF
322 ! - transverse shear stresses
323 signyz(i) = sigoyz(i) + depsyz(i)*g23
324 signzx(i) = sigozx(i) + depszx(i)*g31
325c
326 ! Computation of trial khi coefficients
327 khi1(i) = zero
328 khi2(i) = zero
329 khi3(i) = zero
330 khi4(i) = zero
331 khi5(i) = zero
332 khi6(i) = zero
333 IF ((n1(1)*signxx(i)+n1(2)*signyy(i)) > zero) khi1(i) = one
334 IF ((n2(1)*signxx(i)+n2(2)*signyy(i)) > zero) khi2(i) = one
335 IF (n3*signxy(i) > zero) khi3(i) = one
336 IF ((n4(1)*signxx(i)+n4(2)*signyy(i)) > zero) khi4(i) = one
337 IF ((n5(1)*signxx(i)+n5(2)*signyy(i)) > zero) khi5(i) = one
338 IF (n6*signxy(i) > zero) khi6(i) = one
339c
340 ! Computation of the in-plane yield function
341 gam1(i) = (n1(1)*signxx(i)+n1(2)*signyy(i))/sigy1(i)
342 gam2(i) = (n2(1)*signxx(i)+n2(2)*signyy(i))/sigy2(i)
343 gam3(i) = n3*signxy(i)/sigy3(i)
344 gam4(i) = (n4(1)*signxx(i)+n4(2)*signyy(i))/sigy4(i)
345 gam5(i) = (n5(1)*signxx(i)+n5(2)*signyy(i))/sigy5(i)
346 gam6(i) = n6*signxy(i)/sigy3(i)
347 phi(i) = khi1(i)*gam1(i)**deuxk + khi2(i)*gam2(i)**deuxk +
348 . khi3(i)*gam3(i)**deuxk + khi4(i)*gam4(i)**deuxk +
349 . khi5(i)*gam5(i)**deuxk + khi6(i)*gam6(i)**deuxk - one
350c
351 ! Computation of the out-of-plane shear yield function
352 theta(i) = - signzz(i) - sigyo(i)
353c
354 ! Computation of the transverse shear yield function
355 psi(i) = (sqrt(signyz(i)**2 + signzx(i)**2)/(sigys(i))) - one
356c
357 ENDDO
358c
359 !========================================================================
360 ! - CHECKING THE PLASTIC BEHAVIOR
361 !========================================================================
362c
363 ! Checking plastic behavior
364 nindx = 0
365 nindx2 = 0
366 nindx3 = 0
367 DO i=1,nel
368 ! In plane
369 IF (phi(i) > zero .AND. off(i) == one) THEN
370 nindx=nindx+1
371 index(nindx)=i
372 ENDIF
373 ! Out-of-plane
374 IF (theta(i) > zero .AND. off(i) == one) THEN
375 nindx2=nindx2+1
376 index2(nindx2)=i
377 ENDIF
378 ! Transverse shear
379 IF (psi(i) > zero .AND. off(i) == one) THEN
380 nindx3=nindx3+1
381 index3(nindx3)=i
382 ENDIF
383 ENDDO
384c
385 !====================================================================
386 ! - I IN-PLANE PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON RESOLUTION)
387 !====================================================================
388c
389 ! Number of Newton iterations
390 niter = 3
391c
392 ! Loop over the iterations
393 DO iter = 1, niter
394#include "vectorize.inc"
395 ! Loop over yielding elements
396 DO ii=1,nindx
397 i = index(ii)
398c
399 ! Note : in this part, the purpose is to compute for each iteration
400 ! a plastic multiplier allowing to update internal variables to satisfy
401 ! the consistency condition using the cutting plane method
402 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
403 ! -> PHI : current value of yield function (known)
404 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
405 ! into account of internal variables kinetic :
406 ! hardening parameters
407c
408 ! 1 - Computation of DPHI_DSIG the normal to the yield criterion
409 !-------------------------------------------------------------
410 ! Computation of derivatives to the yield criterion
411 ! - in-plane normal
412 normxx = deuxk*khi1(i)*(gam1(i)**(deuxk-1))*(n1(1)/sigy1(i)) +
413 . deuxk*khi2(i)*(gam2(i)**(deuxk-1))*(n2(1)/sigy2(i)) +
414 . deuxk*khi4(i)*(gam4(i)**(deuxk-1))*(n4(1)/sigy4(i)) +
415 . deuxk*khi5(i)*(gam5(i)**(deuxk-1))*(n5(1)/sigy5(i))
416 normyy = deuxk*khi1(i)*(gam1(i)**(deuxk-1))*(n1(2)/sigy1(i)) +
417 . deuxk*khi2(i)*(gam2(i)**(deuxk-1))*(n2(2)/sigy2(i)) +
418 . deuxk*khi4(i)*(gam4(i)**(deuxk-1))*(n4(2)/sigy4(i)) +
419 . deuxk*khi5(i)*(gam5(i)**(deuxk-1))*(n5(2)/sigy5(i))
420 normxy = deuxk/two*khi3(i)*(gam3(i)**(deuxk-1))*(n3/sigy3(i)) +
421 . deuxk/two*khi6(i)*(gam6(i)**(deuxk-1))*(n6/sigy3(i))
422c
423 ! Normalization of the derivatives
424 ! - in-plane normal
425 normsig = sqrt(normxx**2 + normyy**2 + two*normxy**2)
426 normsig = max(normsig,em20)
427 normpxx = normxx/normsig
428 normpyy = normyy/normsig
429 normpxy = normxy/normsig
430c
431 ! 3 - Computation of DPHI_DLAMBDA
432 !---------------------------------------------------------
433c
434 ! a) Derivative with respect stress increments tensor DSIG
435 ! --------------------------------------------------------
436 dfdsig2 = normxx * (a11*normpxx + a12*normpyy)
437 . + normyy * (a21*normpxx + a22*normpyy)
438 . + two*normxy * two*normpxy * g12
439c
440 ! b) Derivatives with respect to hardening parameters
441 ! ---------------------------------------------------
442c
443 ! i) Derivatives of PHI with respect to the yield stresses
444 dphi_dsigy1 = deuxk*khi1(i)*(gam1(i)**(deuxk-1))*(-(n1(1)*signxx(i)+n1(2)*signyy(i))/(sigy1(i)**2))
445 dphi_dsigy2 = deuxk*khi2(i)*(gam2(i)**(deuxk-1))*(-(n2(1)*signxx(i)+n2(2)*signyy(i))/(sigy2(i)**2))
446 dphi_dsigy3 = deuxk*khi3(i)*(gam3(i)**(deuxk-1))*(-n3*signxy(i)/(sigy3(i)**2)) +
447 . deuxk*khi6(i)*(gam6(i)**(deuxk-1))*(-n6*signxy(i)/(sigy3(i)**2))
448 dphi_dsigy4 = deuxk*khi4(i)*(gam4(i)**(deuxk-1))*(-(n4(1)*signxx(i)+n4(2)*signyy(i))/(sigy4(i)**2))
449 dphi_dsigy5 = deuxk*khi5(i)*(gam5(i)**(deuxk-1))*(-(n5(1)*signxx(i)+n5(2)*signyy(i))/(sigy5(i)**2))
450 ! ii) Derivatives of yield stresses with respect to hardening parameters
451 IF (itab == 0) THEN
452 dsigy1_dp(i) = a01*b01*(one-(tanh(b01*pla(i,2)))**2) + c01
453 dsigy2_dp(i) = a02*b02*(one-(tanh(b02*pla(i,2)))**2) + c02
454 dsigy3_dp(i) = a03*b03*(one-(tanh(b03*pla(i,2)))**2) + c03
455 dsigy4_dp(i) = a04*b04*(one-(tanh(b04*pla(i,2)))**2) + c04
456 dsigy5_dp(i) = a05*b05*(one-(tanh(b05*pla(i,2)))**2) + c05
457 ENDIF
458 ! iii) Assembling derivatives of PHI with respect to hardening parameter
459 hardp(i) = sqrt(khi1(i)*dsigy1_dp(i)**2 + khi2(i)*dsigy2_dp(i)**2 +
460 . khi3(i)*two*dsigy3_dp(i)**2 + khi4(i)*dsigy4_dp(i)**2 +
461 . khi5(i)*dsigy5_dp(i)**2)
462 dphi_dpla = dphi_dsigy1*dsigy1_dp(i) + dphi_dsigy2*dsigy2_dp(i) +
463 . dphi_dsigy3*dsigy3_dp(i) + dphi_dsigy4*dsigy4_dp(i) +
464 . dphi_dsigy5*dsigy5_dp(i)
465c
466 ! iv) Derivative of PHI with respect to DLAM ( = -DENOM)
467 dphi_dlam = -dfdsig2 + dphi_dpla
468 dphi_dlam = sign(max(abs(dphi_dlam),em20) ,dphi_dlam)
469c
470 ! 4 - Computation of plastic multiplier and variables update
471 !----------------------------------------------------------
472c
473 ! a) Computation of the plastic multiplier
474 dlam = -phi(i) / dphi_dlam
475c
476 ! b) Plastic strains tensor update
477 dpxx = dlam * normpxx
478 dpyy = dlam * normpyy
479 dpxy = two * dlam * normpxy
480c
481 ! c) Elasto-plastic stresses update
482 signxx(i) = signxx(i) - (a11*dpxx + a12*dpyy)
483 signyy(i) = signyy(i) - (a21*dpxx + a22*dpyy)
484 signxy(i) = signxy(i) - g12*dpxy
485c
486 ! d) Cumulated plastic strain and hardening parameter update
487 ddep = dlam
488 dpla2(i) = max(zero, dpla2(i) + ddep)
489 pla(i,2) = pla(i,2) + ddep
490c
491 ! e) Yield stresses update
492 IF (itab == 0) THEN
493 ! Tensile yield stress in direction 1 (MD)
494 sigy1(i) = s01 + a01*tanh(b01*pla(i,2)) + c01*pla(i,2)
495 ! Tensile yield stress in direction 2 (CD)
496 sigy2(i) = s02 + a02*tanh(b02*pla(i,2)) + c02*pla(i,2)
497 ! Positive shear yield stress
498 sigy3(i) = s03 + a03*tanh(b03*pla(i,2)) + c03*pla(i,2)
499 ! compressive yield stress in direction 1 (md)
500 sigy4(i) = s04 + a04*tanh(b04*pla(i,2)) + c04*pla(i,2)
501 ! Compressive yield stress in direction 2 (CD)
502 sigy5(i) = s05 + a05*tanh(b05*pla(i,2)) + c05*pla(i,2)
503 ENDIF
504c
505 ! i) Update of trial alpha coefficients
506 khi1(i) = zero
507 khi2(i) = zero
508 khi3(i) = zero
509 khi4(i) = zero
510 khi5(i) = zero
511 khi6(i) = zero
512 IF ((n1(1)*signxx(i)+n1(2)*signyy(i)) > zero) khi1(i) = one
513 IF ((n2(1)*signxx(i)+n2(2)*signyy(i)) > zero) khi2(i) = one
514 IF (n3*signxy(i) > zero) khi3(i) = one
515 IF ((n4(1)*signxx(i)+n4(2)*signyy(i)) > zero) khi4(i) = one
516 IF ((n5(1)*signxx(i)+n5(2)*signyy(i)) > zero) khi5(i) = one
517 IF (n6*signxy(i) > zero) khi6(i) = one
518c
519 ! j) Yield function value update
520 IF (itab == 0) THEN
521 gam1(i) = (n1(1)*signxx(i)+n1(2)*signyy(i))/sigy1(i)
522 gam2(i) = (n2(1)*signxx(i)+n2(2)*signyy(i))/sigy2(i)
523 gam3(i) = n3*signxy(i)/sigy3(i)
524 gam4(i) = (n4(1)*signxx(i)+n4(2)*signyy(i))/sigy4(i)
525 gam5(i) = (n5(1)*signxx(i)+n5(2)*signyy(i))/sigy5(i)
526 gam6(i) = n6*signxy(i)/sigy3(i)
527 phi(i) = khi1(i)*gam1(i)**deuxk + khi2(i)*gam2(i)**deuxk +
528 . khi3(i)*gam3(i)**deuxk + khi4(i)*gam4(i)**deuxk +
529 . khi5(i)*gam5(i)**deuxk + khi6(i)*gam6(i)**deuxk - one
530 ENDIF
531c
532 ENDDO
533 ! End of the loop over yielding elements
534c
535 ! If tabulated yield function, update of the yield stress for all element
536 IF (itab > 0) THEN
537 IF (nindx > 0) THEN
538 ! In-plane tabulation with strain-rate
539 xvec(1:nel,1) = pla(1:nel,2)
540 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
541 ! -> Tensile yield stress in direction 1 (MD)
542 ipos(1:nel,1) = vartmp(1:nel,1)
543 ipos(1:nel,2) = 1
544 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,sigy1,dsigy1_dp,hardr)
545 sigy1(1:nel) = sigy1(1:nel) * yscale1
546 dsigy1_dp(1:nel)= dsigy1_dp(1:nel) * yscale1
547 vartmp(1:nel,1) = ipos(1:nel,1)
548 ! -> Tensile yield stress in direction 2 (CD)
549 xvec(1:nel,2) = epsd(1:nel,2) * xscale2
550 ipos(1:nel,1) = vartmp(1:nel,2)
551 ipos(1:nel,2) = 1
552 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,sigy2,dsigy2_dp,hardr)
553 sigy2(1:nel) = sigy2(1:nel) * yscale2
554 dsigy2_dp(1:nel)= dsigy2_dp(1:nel) * yscale2
555 vartmp(1:nel,2) = ipos(1:nel,1)
556 ! -> Positive shear yield stress
557 xvec(1:nel,2) = epsd(1:nel,2) * xscale3
558 ipos(1:nel,1) = vartmp(1:nel,3)
559 ipos(1:nel,2) = 1
560 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,sigy3,dsigy3_dp,hardr)
561 sigy3(1:nel) = sigy3(1:nel) * yscale3
562 dsigy3_dp(1:nel)= dsigy3_dp(1:nel) * yscale3
563 vartmp(1:nel,3) = ipos(1:nel,1)
564 ! -> Compressive yield stress in direction 1 (MD)
565 xvec(1:nel,2) = epsd(1:nel,2) * xscale4
566 ipos(1:nel,1) = vartmp(1:nel,4)
567 ipos(1:nel,2) = 1
568 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,sigy4,dsigy4_dp,hardr)
569 sigy4(1:nel) = sigy4(1:nel) * yscale4
570 dsigy4_dp(1:nel)= dsigy4_dp(1:nel) * yscale4
571 vartmp(1:nel,4) = ipos(1:nel,1)
572 ! -> Compressive yield stress in direction 2 (CD)
573 xvec(1:nel,2) = epsd(1:nel,2) * xscale5
574 ipos(1:nel,1) = vartmp(1:nel,5)
575 ipos(1:nel,2) = 1
576 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,sigy5,dsigy5_dp,hardr)
577 sigy5(1:nel) = sigy5(1:nel) * yscale5
578 dsigy5_dp(1:nel)= dsigy5_dp(1:nel) * yscale5
579 vartmp(1:nel,5) = ipos(1:nel,1)
580 ! Yield function value update
581#include "vectorize.inc"
582 DO ii=1,nindx
583 i = index(ii)
584 gam1(i) = (n1(1)*signxx(i)+n1(2)*signyy(i))/sigy1(i)
585 gam2(i) = (n2(1)*signxx(i)+n2(2)*signyy(i))/sigy2(i)
586 gam3(i) = n3*signxy(i)/sigy3(i)
587 gam4(i) = (n4(1)*signxx(i)+n4(2)*signyy(i))/sigy4(i)
588 gam5(i) = (n5(1)*signxx(i)+n5(2)*signyy(i))/sigy5(i)
589 gam6(i) = n6*signxy(i)/sigy3(i)
590 phi(i) = khi1(i)*gam1(i)**deuxk + khi2(i)*gam2(i)**deuxk +
591 . khi3(i)*gam3(i)**deuxk + khi4(i)*gam4(i)**deuxk +
592 . khi5(i)*gam5(i)**deuxk + khi6(i)*gam6(i)**deuxk - one
593 ENDDO
594 ENDIF
595 ENDIF
596 ENDDO
597 ! End of the loop over the iterations
598 !============================================================================
599 ! - END OF IN-PLANE PLASTIC CORRECTION WITH NEWTON IMPLICIT ITERATIVE METHOD
600 !============================================================================
601c
602 !============================================================================
603 ! - II OUT-OF-PLANE PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON RESOLUTION)
604 !============================================================================
605c
606 ! Loop over the iterations
607 DO iter = 1, niter
608#include "vectorize.inc"
609 ! Loop over yielding elements
610 DO ii=1,nindx2
611 i = index2(ii)
612c
613 ! Note : in this part, the purpose is to compute for each iteration
614 ! a plastic multiplier allowing to update internal variables to satisfy
615 ! the consistency condition using the cutting plane method
616 ! Its expression at each iteration is : DLAMBDA = - THETA/DTHETA_DLAMBDA
617 ! -> THETA : current value of yield function (known)
618 ! -> DTHETA_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
619 ! into account of internal variables kinetic :
620 ! hardening parameters
621c
622 ! 1 - Computation of DTHETA_DSIG the normal to the yield criterion
623 !-------------------------------------------------------------
624 ! Computation of derivatives to the yield criterion
625 ! - in-plane normal
626 normzz = -one
627c
628 ! 3 - Computation of DTHETA_DLAMBDA
629 !---------------------------------------------------------
630c
631 ! a) Derivative with respect stress increments tensor DSIG
632 ! --------------------------------------------------------
633 IF (eezz(i) >= zero) THEN
634 dfdsig2 = normzz * young3 * normzz
635 ELSE
636 dfdsig2 = normzz * e3c*cc*exp(-cc*eezz(i)) * normzz
637 ENDIF
638c
639 ! b) Derivatives with respect to hardening parameters
640 ! ---------------------------------------------------
641c
642 ! i) Derivatives of THETA with respect to the yield stresses
643 dtheta_dsigyo = -one
644 ! ii) Derivatives of yield stresses with respect to hardening parameters
645 IF (itab == 0) dsigyo_dp(i) = csig*bsig*exp(csig*pla(i,3))
646 ! iii) Assembling derivatives of THETA with respect to hardening parameter
647 dtheta_dpla = dtheta_dsigyo*dsigyo_dp(i)
648c
649 ! iv) Derivative of THETA with respect to DLAM ( = -DENOM)
650 dtheta_dlam = -dfdsig2 + dtheta_dpla
651 dtheta_dlam = sign(max(abs(dtheta_dlam),em20) ,dtheta_dlam)
652c
653 ! 4 - Computation of plastic multiplier and variables update
654 !----------------------------------------------------------
655c
656 ! a) Computation of the plastic multiplier
657 dlam = -theta(i) / dtheta_dlam
658c
659 ! c) Elasto-plastic stresses update
660 IF (eezz(i)>=zero) THEN
661 signzz(i) = signzz(i) - young3*dlam*normzz
662 ELSE
663 signzz(i) = signzz(i) - e3c*cc*exp(-cc*eezz(i))*dlam*normzz
664 ENDIF
665c
666 ! d) Cumulated plastic strain and hardening parameter update
667 ddep = dlam
668 dpla3(i) = max(zero, dpla3(i) + ddep)
669 pla(i,3) = pla(i,3) + ddep
670 eezz(i) = epszz(i) + pla(i,3)
671c
672 ! e) Yield stresses update
673 IF (itab == 0) THEN
674 ! Out-of-plane yield stress
675 sigyo(i) = asig + bsig*exp(csig*pla(i,3))
676 ! Yield function value update
677 theta(i) = - signzz(i) - sigyo(i)
678 ENDIF
679c
680 ENDDO
681 ! End of the loop over yielding elements
682c
683 ! If tabulated yield function, update of the yield stress for all element
684 IF (itab > 0) THEN
685 IF (nindx2 > 0) THEN
686 ! -> Transverse shear tabulation with strain-rate
687 xvec(1:nel,1) = pla(1:nel,3)
688 xvec(1:nel,2) = epsd(1:nel,3) * xscalec
689 ! -> Transverse shear yield stress
690 ipos(1:nel,1) = vartmp(1:nel,6)
691 ipos(1:nel,2) = 1
692 CALL table2d_vinterp_log(table(itable(6)),ismooth,nel,nel,ipos,xvec,sigyo,dsigyo_dp,hardr)
693 sigyo(1:nel) = sigyo(1:nel) * yscalec
694 dsigyo_dp(1:nel)= dsigyo_dp(1:nel) * yscalec
695 vartmp(1:nel,6) = ipos(1:nel,1)
696 ! Yield function value update
697#include "vectorize.inc"
698 DO ii=1,nindx2
699 i = index2(ii)
700 theta(i) = - signzz(i) - sigyo(i)
701 ENDDO
702 ENDIF
703 ENDIF
704 ENDDO
705 ! End of the loop over the iterations
706 !===============================================================================
707 ! - END OF OUT-OF-PLANE PLASTIC CORRECTION WITH NEWTON IMPLICIT ITERATIVE METHOD
708 !===============================================================================
709c
710 !================================================================================
711 ! - III TRANSVERSE SHEAR PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON RESOLUTION)
712 !================================================================================
713c
714 ! Number of Newton iterations
715 niter = 3
716c
717 ! Loop over the iterations
718 DO iter = 1, niter
719#include "vectorize.inc"
720 ! Loop over yielding elements
721 DO ii=1,nindx3
722 i = index3(ii)
723c
724 ! Note : in this part, the purpose is to compute for each iteration
725 ! a plastic multiplier allowing to update internal variables to satisfy
726 ! the consistency condition using the cutting plane method
727 ! Its expression at each iteration is : DLAMBDA = - PSI/DPSI_DLAMBDA
728 ! -> PSI : current value of yield function (known)
729 ! -> DPSI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
730 ! into account of internal variables kinetic :
731 ! hardening parameters
732c
733 ! 1 - Computation of DPSI_DSIG the normal to the yield criterion
734 !-------------------------------------------------------------
735 ! Computation of derivatives to the yield criterion
736 ! - transverse shear normal
737 normyz = signyz(i)/max((sqrt(signyz(i)**2 + signzx(i)**2)*sigys(i)),em20)
738 normzx = signzx(i)/max((sqrt(signyz(i)**2 + signzx(i)**2)*sigys(i)),em20)
739c
740 ! 2 - Computation of DPSI_DLAMBDA
741 !---------------------------------------------------------
742c
743 ! a) Derivative with respect stress increments tensor DSIG
744 ! --------------------------------------------------------
745 dfdsig2 = two*normyz * g23 * two*normyz +
746 . two*normzx * g31 * two*normzx
747c
748 ! b) Derivatives with respect to hardening parameters
749 ! ---------------------------------------------------
750c
751 ! i) Derivatives of PSI with respect to the yield stresses
752 dpsi_dsigys = -sqrt(signyz(i)**2 + signzx(i)**2)/(sigys(i)**2)
753 ! ii) Derivatives of yield stresses with respect to hardening parameters
754 IF (itab == 0) dsigys_dp(i) = atau - min(zero,sigozz(i))*btau
755 ! iii) Assembling derivatives of PSI with respect to hardening parameter
756 dpsi_dpla = dpsi_dsigys*dsigys_dp(i)
757c
758 ! iv) Derivative of PSI with respect to DLAM ( = -DENOM)
759 dpsi_dlam = -dfdsig2 + dpsi_dpla
760 dpsi_dlam = sign(max(abs(dpsi_dlam),em20),dpsi_dlam)
761c
762 ! 4 - Computation of plastic multiplier and variables update
763 !----------------------------------------------------------
764c
765 ! a) Computation of the plastic multiplier
766 dlam = -psi(i) / dpsi_dlam
767c
768 ! b) Plastic strains tensor update
769 dpyz = two*dlam * normyz
770 dpzx = two*dlam * normzx
771c
772 ! c) Elasto-plastic stresses update
773 signyz(i) = signyz(i) - g23*dpyz
774 signzx(i) = signzx(i) - g31*dpzx
775c
776 ! d) Cumulated plastic strain and hardening parameter update
777 ddep = dlam
778 dpla4(i) = max(zero, dpla4(i) + ddep)
779 pla(i,4) = pla(i,4) + ddep
780c
781 ! e) Yield stresses update
782 IF (itab == 0) THEN
783 ! Transverse shear yield stress
784 sigys(i) = tau0 + (atau - min(zero,sigozz(i))*btau)*pla(i,4)
785 ! f) Yield function value update
786 psi(i) = (sqrt(signyz(i)**2 + signzx(i)**2)/sigys(i)) - one
787 ENDIF
788c
789 ENDDO
790 ! End of the loop over yielding elements
791c
792 ! If tabulated yield function, update of the yield stress for all element
793 IF (itab > 0) THEN
794 IF (nindx3 > 0) THEN
795 ! Transverse shear tabulation with strain-rate
796 xvec(1:nel,1) = pla(1:nel,4)
797 xvec(1:nel,2) = epsd(1:nel,4) * xscales
798 ! transverse shear yield stress
799 ipos(1:nel,1) = vartmp(1:nel,7)
800 ipos(1:nel,2) = 1
801 CALL table2d_vinterp_log(table(itable(7)),ismooth,nel,nel,ipos,xvec,sigys,dsigys_dp,hardr)
802 sigys(1:nel) = sigys(1:nel) * yscales
803 dsigys_dp(1:nel)= dsigys_dp(1:nel) * yscales
804 vartmp(1:nel,7) = ipos(1:nel,1)
805 ! Yield function value update
806#include "vectorize.inc"
807 DO ii=1,nindx3
808 i = index3(ii)
809 psi(i) = (sqrt(signyz(i)**2 + signzx(i)**2)/sigys(i)) - one
810 ENDDO
811 ENDIF
812 ENDIF
813 ENDDO
814 ! End of the loop over the iterations
815 !===================================================================================
816 ! - END OF TRANSVERSE SHEAR PLASTIC CORRECTION WITH NEWTON IMPLICIT ITERATIVE METHOD
817 !===================================================================================
818c
819 ! Storing new values
820 DO i=1,nel
821 ! Plastic strain
822 pla(i,1) = sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2)
823 dpla(i) = (pla(i,2)/(max(sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2),em20)))*dpla2(i) +
824 . (pla(i,3)/(max(sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2),em20)))*dpla3(i) +
825 . (pla(i,4)/(max(sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2),em20)))*dpla4(i)
826 ! Plastic strain-rate
827 IF (itab == 0) THEN
828 epsd(i,1) = dpla(i) / max(em20,timestep)
829 epsd(i,2) = dpla2(i) / max(em20,timestep)
830 epsd(i,3) = dpla3(i) / max(em20,timestep)
831 epsd(i,4) = dpla4(i) / max(em20,timestep)
832 ELSE
833 dpdt = dpla(i)/max(em20,timestep)
834 epsd(i,1) = asrate * dpdt + (one - asrate) * epsd(i,1)
835 dpdt = dpla2(i)/max(em20,timestep)
836 epsd(i,2) = asrate * dpdt + (one - asrate) * epsd(i,2)
837 dpdt = dpla3(i)/max(em20,timestep)
838 epsd(i,3) = asrate * dpdt + (one - asrate) * epsd(i,3)
839 dpdt = dpla4(i)/max(em20,timestep)
840 epsd(i,4) = asrate * dpdt + (one - asrate) * epsd(i,4)
841 ENDIF
842 ! Coefficient for hourglass and soundspeed
843 IF (eezz(i) >= zero) THEN
844 IF (dpla(i) > zero) THEN
845 et(i) = hardp(i) / (hardp(i) + max(young1,young2,young3))
846 ELSE
847 et(i) = one
848 ENDIF
849 soundsp(i) = sqrt(max(a11,a12,a21,a22,young3,g12,g23,g31)/ rho(i))
850 ELSE
851 IF (dpla(i) > zero) THEN
852 et(i) = hardp(i) / (hardp(i) + max(young1,young2,e3c*cc*exp(-cc*eezz(i))))
853 ELSE
854 et(i) = one
855 ENDIF
856 soundsp(i) = sqrt(max(a11,a12,a21,a22,e3c*cc*exp(-cc*eezz(i)),g12,g23,g31)/ rho(i))
857 ENDIF
858 ! Storing the yield stress
859 sigy(i) = sqrt(khi1(i)*sigy1(i)**2 + khi2(i)*sigy2(i)**2 + khi4(i)*sigy4(i)**2 +
860 . khi5(i)*sigy5(i)**2 + two*khi3(i)*sigy3(i)**2 + sigyo(i)**2
861 . + two*sigys(i)**2)
862 ENDDO
863c
864 END
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat112_xia_newton
subroutine mat112_xia_newton(nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, jthe, off, rho0, rho, pla, dpla, epsd, soundsp, epszz, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, nvartmp, numtabl, vartmp, itable, table)
Definition mat112_xia_newton.F:42
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112
table_mod::ttable
Definition table_mod.F:125
table2d_vinterp_log
subroutine table2d_vinterp_log(table, ismooth, dimx, nel, ipos, xx, yy, dydx1, dydx2)
Definition table2d_vinterp_log.F:52