OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
hm_read_mat110.F File Reference
#include "implicit_f.inc"
#include "units_c.inc"
#include "param_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine hm_read_mat110 (uparam, maxuparam, nuparam, nuvar, ntabl, mtag, parmat, unitab, pm, lsubmodel, israte, mat_id, titr, itable, maxtabl, nvartmp, matparam)

Function/Subroutine Documentation

◆ hm_read_mat110()

subroutine hm_read_mat110 ( intent(inout) uparam,
integer, intent(in) maxuparam,
integer, intent(inout) nuparam,
integer, intent(inout) nuvar,
integer, intent(inout) ntabl,
type(mlaw_tag_), intent(inout) mtag,
intent(inout) parmat,
type (unit_type_), intent(in) unitab,
intent(inout) pm,
type(submodel_data), dimension(*), intent(in) lsubmodel,
integer, intent(inout) israte,
integer, intent(in) mat_id,
character(len=nchartitle), intent(in) titr,
integer, dimension(maxtabl), intent(inout) itable,
integer, intent(in) maxtabl,
integer, intent(inout) nvartmp,
type(matparam_struct_), intent(inout) matparam )

Definition at line 40 of file hm_read_mat110.F.

45C-----------------------------------------------
46C M o d u l e s
47C-----------------------------------------------
48 USE unitab_mod
49 USE message_mod
50 USE submodel_mod
51 USE elbuftag_mod
52 USE matparam_def_mod
53 USE names_and_titles_mod , ONLY : nchartitle
54C-----------------------------------------------
55C I m p l i c i t T y p e sXM
56C-----------------------------------------------
57#include "implicit_f.inc"
58C-----------------------------------------------
59C C o m m o n B l o c k s
60C-----------------------------------------------
61#include "units_c.inc"
62#include "param_c.inc"
63C-----------------------------------------------
64C D u m m y A r g u m e n t s
65C-----------------------------------------------
66 TYPE (UNIT_TYPE_),INTENT(IN) ::UNITAB
67 INTEGER, INTENT(IN) :: MAT_ID,MAXUPARAM,MAXTABL
68 my_real, DIMENSION(NPROPM) ,INTENT(INOUT) :: pm
69 CHARACTER(LEN=NCHARTITLE) ,INTENT(IN) :: TITR
70 INTEGER, INTENT(INOUT) :: ISRATE,ITABLE(MAXTABL)
71 INTEGER, INTENT(INOUT) :: NUPARAM,NUVAR,NTABL,NVARTMP
72 my_real, DIMENSION(MAXUPARAM) ,INTENT(INOUT) :: uparam
73 my_real, DIMENSION(100),INTENT(INOUT) :: parmat
74 TYPE(SUBMODEL_DATA), DIMENSION(*),INTENT(IN) :: LSUBMODEL
75 TYPE(MLAW_TAG_), INTENT(INOUT) :: MTAG
76 TYPE(MATPARAM_STRUCT_),INTENT(INOUT) :: MATPARAM
77C
78C-----------------------------------------------
79C L o c a l V a r i a b l e s
80C-----------------------------------------------
81 INTEGER I, J, K, ILAW, Ivflag, NANGLE, INFO, Icrit, TAB_YLD, Ismooth,
82 . TAB_TEMP,ITER,Ires
83C REAL
84 my_real
85 . rho0,young,nu,yld0,dsigm,beta,omega,nexp,eps0,sigstar,dg0,
86 . deps0,mexp,fbi,rhobi,rhor,bulk,a11,g,n1,n2,m1,m2,a1,a2,c1,c2,
87 . asrate,kboltz,xscale,yscale,fisokin,tini,rhobi_theta,b1,b2,mu,
88 . xscale_unit,yscale_unit
89 my_real
90 . f1,f2,df1_dmu,df2_dmu,d2f1_d2mu,d2f2_d2mu,
91 . dphi_dsig1,dphi_dsig2,dmu_dsig1,dmu_dsig2,mineig,thet,da_dcos2(2),
92 . db_dcos2(2),dc_dcos2(2),d2a_d2cos2(2),d2b_d2cos2(2),d2c_d2cos2(2),
93 . df1_dcos2,df2_dcos2,dphi_dcos2,
94 . d2f1_d2cos2,d2f2_d2cos2,d2f1_d2mucos2,d2f2_d2mucos2,mineigmu,
95 . mineigthet,denom,
96 . dmu_dcos2,t(3,3),d(3,3),weight,dweight_dcos2,d2weight_d2cos2,u,uprim,upprim,
97 . du_dmu,duprim_dmu,v,vprim,vpprim,dv_dmu,dvprim_dmu,w0,fun_av,rlank_av,
98 . fbi_min,fbi_max,fbi_trans,adirac,fdirac,fps_min,fps_max,fps_trans,
99 . fsh_min,fsh_max,fsh_trans,fh(2),fps1,rhoone_theta,phi,dfps1_dweight
100 real*8 :: cos2(10,10),d2phi_d2(3,3),work(102),wr(3),wi(3),vl(3,3),vr(3,3)
101 real*8 , DIMENSION(:,:), ALLOCATABLE :: amat,bvec
102 my_real, DIMENSION(:,:), ALLOCATABLE :: fun,fsh,fps,hips,hiun,hish,
103 . q_fsh,q_fun,q_fps,q_hiun,q_hips,q_hish
104 my_real, DIMENSION(:) , ALLOCATABLE :: rlank,theta,theta_rad,
105 . wps,wsh,q_wsh,q_wps,alps,rm,ag,
106 . epsag,sigag,cag,epseq
107 INTEGER, DIMENSION(:), ALLOCATABLE :: IPIV
108 my_real, PARAMETER :: tol = 1.0d-6
109C
110 LOGICAL :: IS_AVAILABLE,IS_ENCRYPTED
111C
112 DATA cos2/
113 1 1. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
114 2 0. ,1. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
115 3 -1. ,0. ,2. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,
116 4 0. ,-3. ,0. ,4. ,0. ,0. ,0. ,0. ,0. ,0. ,
117 5 1. ,0. ,-8. ,0. ,8. ,0. ,0. ,0. ,0. ,0. ,
118 6 0. ,5. ,0. ,-20. ,0. ,16. ,0. ,0. ,0. ,0. ,
119 7 -1. ,0. ,18. ,0. ,-48. ,0. ,32. ,0. ,0. ,0. ,
120 8 0. ,-7. ,0. ,56. ,0. ,-112.,0. ,64. ,0. ,0. ,
121 9 1. ,0. ,-32. ,0. ,160. ,0. ,-256.,0. ,128. ,0. ,
122 a 0. ,9. ,0. ,-120.,0. ,432. ,0. ,-576 ,0. ,256. /
123C=======================================================================
124 is_encrypted = .false.
125 is_available = .false.
126 ilaw = 110
127c------------------------------------------
128 CALL hm_option_is_encrypted(is_encrypted)
129c------------------------------------------
130c
131card1 - Density
132 CALL hm_get_floatv('MAT_RHO' ,rho0 ,is_available, lsubmodel, unitab)
133 CALL hm_get_floatv('Refer_Rho' ,rhor ,is_available, lsubmodel, unitab)
134card2 - Elasticity
135 CALL hm_get_floatv('MAT_E' ,young ,is_available, lsubmodel, unitab)
136 CALL hm_get_floatv('MAT_NU' ,nu ,is_available, lsubmodel, unitab)
137 CALL hm_get_intv ('MAT_Ires' ,ires ,is_available, lsubmodel)
138card3
139 CALL hm_get_intv ('MAT_Icrit' ,icrit ,is_available, lsubmodel)
140 ! Default : classical Vegter yield locus
141 IF (icrit == 0) icrit = 1
142 IF (icrit > 4) THEN
143 CALL ancmsg(msgid=1776,msgtype=msgerror,
144 . anmode=aninfo_blind_1,i1=mat_id,c1=titr,
145 . i2=icrit)
146 ENDIF
147 CALL hm_get_intv ('MAT_TAB_YLD',tab_yld ,is_available, lsubmodel)
148 CALL hm_get_floatv('MAT_Xscale',xscale ,is_available, lsubmodel, unitab)
149 CALL hm_get_floatv('MAT_Yscale',yscale ,is_available, lsubmodel, unitab)
150 CALL hm_get_floatv('MAT_fBI' ,fbi ,is_available, lsubmodel, unitab)
151 CALL hm_get_floatv('MAT_rhoBI' ,rhobi ,is_available, lsubmodel, unitab)
152c
153card4 - Hardening
154 CALL hm_get_floatv('SIGMA_r' ,yld0 ,is_available, lsubmodel, unitab)
155 IF ((tab_yld == 0).AND.(yld0==zero)) THEN
156 CALL ancmsg(msgid=1777,msgtype=msgerror,
157 . anmode=aninfo_blind_1,i1=mat_id,c1=titr)
158 ENDIF
159 CALL hm_get_floatv('MAT_DSIGM' ,dsigm ,is_available, lsubmodel, unitab)
160 CALL hm_get_floatv('MAT_BETA' ,beta ,is_available, lsubmodel, unitab)
161 CALL hm_get_floatv('Omega' ,omega ,is_available, lsubmodel, unitab)
162 CALL hm_get_floatv('mat_hard' ,NEXP ,IS_AVAILABLE, LSUBMODEL, UNITAB)
163c
164card5 - Strain-rate and temperature dependency
165 CALL HM_GET_FLOATV('epsilon_0' ,EPS0 ,IS_AVAILABLE, LSUBMODEL, UNITAB)
166 CALL HM_GET_FLOATV('mat_sigs' ,SIGSTAR ,IS_AVAILABLE, LSUBMODEL, UNITAB)
167 CALL HM_GET_FLOATV('mat_dg0' ,DG0 ,IS_AVAILABLE, LSUBMODEL, UNITAB)
168 CALL HM_GET_FLOATV('mat_deps0' ,Deps0 ,IS_AVAILABLE, LSUBMODEL, UNITAB)
169 CALL HM_GET_FLOATV('mat_strainrate_m',MEXP ,IS_AVAILABLE, LSUBMODEL, UNITAB)
170c
171card6 - Strain-rate computation
172 CALL HM_GET_FLOATV('t_initial' ,TINI ,IS_AVAILABLE, LSUBMODEL, UNITAB)
173 IF (TINI==ZERO) THEN
174 CALL ANCMSG(MSGID=1778,MSGTYPE=MSGWARNING,
175 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,C1=TITR)
176 ENDIF
177 CALL HM_GET_FLOATV('mat_chard' ,FISOKIN ,IS_AVAILABLE, LSUBMODEL, UNITAB)
178 CALL HM_GET_FLOATV('fcut' ,ASRATE ,IS_AVAILABLE, LSUBMODEL, UNITAB)
179 CALL HM_GET_INTV ('vflag' ,Ivflag ,IS_AVAILABLE, LSUBMODEL)
180 CALL HM_GET_INTV ('mat_ismooth',Ismooth ,IS_AVAILABLE, LSUBMODEL)
181 CALL HM_GET_INTV ('mat_tab_temp' ,TAB_TEMP ,IS_AVAILABLE, LSUBMODEL)
182card7 - Yield criterion parameters
183 ! Classical Vegter yield locus
184 NANGLE = 0
185 IF (Icrit == 1) THEN
186 CALL HM_GET_INTV('mat_refanglemax',NANGLE,IS_AVAILABLE,LSUBMODEL)
187 IF (NANGLE > 10) THEN
188 CALL ANCMSG(MSGID=1779,MSGTYPE=MSGERROR,
189 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,C1=TITR)
190 ENDIF
191.NOT. IF (ALLOCATED(fUN)) ALLOCATE(fUN(NANGLE,2))
192.NOT. IF (ALLOCATED(fSH)) ALLOCATE(fSH(NANGLE,2))
193.NOT. IF (ALLOCATED(fPS)) ALLOCATE(fPS(NANGLE,2))
194 fPS(1:NANGLE,1:2) = ZERO
195.NOT. IF (ALLOCATED(rLank)) ALLOCATE(rLank(NANGLE))
196 ! Loop over angles (must be equally distributed between 0 and pi/2)
197 DO J = 1, NANGLE
198 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_fun_theta' ,fUN(J,1),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
199 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_r_theta' ,rLank(J),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
200 IF (rLank(J) == 0) rLank(J) = ONE
201 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_fps1_theta',fPS(J,1),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
202 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_fps2_theta',fPS(J,2),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
203 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_fsh_theta' ,fSH(J,1),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
204 ENDDO
205 ! Standard Vegter yield locus
206 ELSEIF (Icrit == 2) THEN
207 CALL HM_GET_INTV('mat_refanglemax',NANGLE,IS_AVAILABLE,LSUBMODEL)
208 IF (NANGLE > 10) THEN
209 CALL ANCMSG(MSGID=1779,MSGTYPE=MSGERROR,
210 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,C1=TITR)
211 ENDIF
212.NOT. IF (ALLOCATED(fUN)) ALLOCATE(fUN(NANGLE,2))
213.NOT. IF (ALLOCATED(fSH)) ALLOCATE(fSH(NANGLE,2))
214.NOT. IF (ALLOCATED(fPS)) ALLOCATE(fPS(NANGLE,2))
215 fPS(1:NANGLE,1:2) = ZERO
216.NOT. IF (ALLOCATED(ALPS)) ALLOCATE(ALPS(NANGLE))
217.NOT. IF (ALLOCATED(rLank)) ALLOCATE(rLank(NANGLE))
218 ! Loop over angles (must be equally distributed between 0 and pi/2)
219 DO J = 1, NANGLE
220 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_fun_theta' ,fUN(J,1),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
221 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_r_theta' ,rLank(J),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
222 IF (rLank(J) == 0) rLank(J) = ONE
223 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_fps1_theta',fPS(J,1),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
224 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_alps_theta',ALPS(J) ,J,IS_AVAILABLE, LSUBMODEL, UNITAB)
225 IF (ALPS(J) == ZERO) ALPS(J) = HALF
226 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_fsh_theta' ,fSH(J,1),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
227 ENDDO
228 ! Vegter - 2017 yield locus
229 ELSEIF (Icrit == 3) THEN
230 ! Number of angle is fixed to three in this case
231 NANGLE = 3
232.NOT. IF (ALLOCATED(RM)) ALLOCATE(RM(NANGLE))
233.NOT. IF (ALLOCATED(AG)) ALLOCATE(AG(NANGLE))
234.NOT. IF (ALLOCATED(EPSAG)) ALLOCATE(EPSAG(NANGLE))
235.NOT. IF (ALLOCATED(SIGAG)) ALLOCATE(SIGAG(NANGLE))
236.NOT. IF (ALLOCATED(CAG)) ALLOCATE(CAG(NANGLE))
237.NOT. IF (ALLOCATED(EPSEQ)) ALLOCATE(EPSEQ(NANGLE))
238.NOT. IF (ALLOCATED(fUN)) ALLOCATE(fUN(NANGLE,2))
239.NOT. IF (ALLOCATED(fSH)) ALLOCATE(fSH(NANGLE,2))
240.NOT. IF (ALLOCATED(fPS)) ALLOCATE(fPS(NANGLE,2))
241.NOT. IF (ALLOCATED(rLank)) ALLOCATE(rLank(NANGLE))
242.NOT. IF (ALLOCATED(THETA)) ALLOCATE(THETA(NANGLE))
243.NOT. IF (ALLOCATED(THETA_RAD)) ALLOCATE(THETA_RAD(NANGLE))
244 CALL HM_GET_FLOATV('mat_rm_0' ,RM(1) ,IS_AVAILABLE, LSUBMODEL, UNITAB)
245 CALL HM_GET_FLOATV('mat_rm_45' ,RM(2) ,IS_AVAILABLE, LSUBMODEL, UNITAB)
246 CALL HM_GET_FLOATV('mat_rm_90' ,RM(3) ,IS_AVAILABLE, LSUBMODEL, UNITAB)
247 CALL HM_GET_FLOATV('mat_ag_0' ,AG(1) ,IS_AVAILABLE, LSUBMODEL, UNITAB)
248 CALL HM_GET_FLOATV('mat_ag_45' ,AG(2) ,IS_AVAILABLE, LSUBMODEL, UNITAB)
249 CALL HM_GET_FLOATV('mat_ag_90' ,AG(3) ,IS_AVAILABLE, LSUBMODEL, UNITAB)
250 CALL HM_GET_FLOATV('mat_r_0' ,rLank(1) ,IS_AVAILABLE, LSUBMODEL, UNITAB)
251 IF (rLank(1) == ZERO) rLank(1) = ONE
252 CALL HM_GET_FLOATV('mat_r_45' ,rLank(2) ,IS_AVAILABLE, LSUBMODEL, UNITAB)
253 IF (rLank(2) == ZERO) rLank(2) = ONE
254 CALL HM_GET_FLOATV('mat_r_90' ,rLank(3) ,IS_AVAILABLE, LSUBMODEL, UNITAB)
255 IF (rLank(3) == ZERO) rLank(3) = ONE
256c
257 ! Computation of uniaxial reference point
258 ! Uniaxial strain at AG
259 EPSAG(1) = LOG(ONE + AG(1)/HUNDRED)
260 EPSAG(2) = LOG(ONE + AG(2)/HUNDRED)
261 EPSAG(3) = LOG(ONE + AG(3)/HUNDRED)
262 ! Uniaxial stress at AG
263 SIGAG(1) = RM(1)*(ONE + AG(1)/HUNDRED)
264 SIGAG(2) = RM(2)*(ONE + AG(2)/HUNDRED)
265 SIGAG(3) = RM(3)*(ONE + AG(3)/HUNDRED)
266 ! Energy at 0 deg direction
267 CAG(1) = SIGAG(1)/(EPSAG(1)**EPSAG(1))
268 CAG(2) = SIGAG(2)/(EPSAG(2)**EPSAG(2))
269 CAG(3) = SIGAG(3)/(EPSAG(3)**EPSAG(3))
270 W0 = (CAG(1)/(EPSAG(1) + ONE))*((EPSAG(1))**(EPSAG(1)+ONE))
271 ! Equivalent strain
272 EPSEQ(1) = ((EPSAG(1)+ONE)*(W0/CAG(1)))**(ONE/(EPSAG(1)+ONE))
273 EPSEQ(2) = ((EPSAG(2)+ONE)*(W0/CAG(2)))**(ONE/(EPSAG(2)+ONE))
274 EPSEQ(3) = ((EPSAG(3)+ONE)*(W0/CAG(3)))**(ONE/(EPSAG(3)+ONE))
275 ! Reference point
276 fUN(1,1) = EPSAG(1)/EPSEQ(1)
277 fUN(2,1) = EPSAG(1)/EPSEQ(2)
278 fUN(3,1) = EPSAG(1)/EPSEQ(3)
279c
280 ! Computation of biaxial reference point
281 fUN_AV = (fUN(1,1) + TWO*fUN(2,1) + fUN(3,1))/FOUR
282 rLank_AV = (rLank(1) + TWO*rLank(2) + rLank(3))/FOUR
283 fBI_MIN = 0.97D0
284 fBI_MAX = 1.14D0
285 fBI_TRANS = 1.22D0
286 ADIRAC = 3.4D0
287 FDIRAC = EXP(ADIRAC*(rLank_AV-fBI_TRANS))
288 fBI = fUN_AV*((fBI_MIN/(ONE+FDIRAC)) + (fBI_MAX/(ONE+(ONE/FDIRAC))))
289 IF (rhoBI == ZERO) THEN
290 rhoBI = rLank(1)/rLank(3)
291 ENDIF
292c
293 ! Computation of plane strain reference point first coordinate
294 fPS_MIN = 0.827D0
295 fPS_MAX = 1.315D0
296 fPS_TRANS = 0.5D0
297 ADIRAC = 1.2D0
298 DO J = 1,NANGLE
299 FDIRAC = EXP(ADIRAC*(rLank(J)-fPS_TRANS))
300 fPS(J,1) = fUN(J,1)*((fPS_MIN/(ONE+FDIRAC)) + (fPS_MAX/(ONE+(ONE/FDIRAC))))
301 ENDDO
302c
303 ! Computation of shear reference point first coordinate
304 fSH_MIN = 0.757D0
305 fSH_MAX = 0.525D0
306 fSH_TRANS = ZERO
307 ADIRAC = 1.6D0
308 fUN_AV = (fUN(1,1)+fUN(3,1))/TWO
309 rLank_AV = (rLank(1)+rLank(3))/TWO
310 FDIRAC = EXP(ADIRAC*(rLank_AV - fSH_TRANS))
311 fSH(1,1) = fUN_AV*((fSH_MIN/(ONE+FDIRAC)) + (fSH_MAX/(ONE+(ONE/FDIRAC))))
312 fSH(3,1) = fUN_AV*((fSH_MIN/(ONE+FDIRAC)) + (fSH_MAX/(ONE+(ONE/FDIRAC))))
313 FDIRAC = EXP(ADIRAC*(rLank(2)-fSH_TRANS))
314 fSH(2,1) = fUN(2,1)*((fSH_MIN/(ONE+FDIRAC)) + (fSH_MAX/(ONE+(ONE/FDIRAC))))
315c
316 ! Computation of plane strain reference point second coordinate
317 DO J = 1, NANGLE
318 ! Angle theta with respect to rolling direction
319 THETA(J) = (J-1)*(90.0D0/(NANGLE-1))
320 THETA_RAD(J) = THETA(J)*(PI/180.0D0)
321 ! Strain-rate ratio
322 rhoBI_theta = (rhoBI + ONE) + (rhoBI - ONE)*COS(TWO*THETA_RAD(J))
323 rhoBI_theta = rhoBI_theta / ((rhoBI + ONE) - (rhoBI - ONE)*COS(TWO*THETA_RAD(J)))
324 rhoONE_theta = -rLank(J)/(rLank(J)+ONE)
325 ! Virtual hinge point
326 fH(2) = (fBI + rhoBI_theta*fBI - fUN(J,1))/(rhoBI_theta - rhoONE_theta)
327 fH(1) = fBI - rhoBI_theta*(fH(2)-fBI)
328 ! Initialization of the computation of Weight
329 MU = HALF
330 WEIGHT = ZERO
331 fPS1 = (fUN(J,1)*((ONE - MU)**2) + TWO*fH(1)*MU*WEIGHT*(ONE-MU) + fBI*(MU**2))/
332 . ((ONE-MU)**2 + TWO*WEIGHT*MU*(ONE-MU) + MU**2)
333 PHI = fPS1-fPS(J,1)
334 ! Newton iteration to compute WEIGHT
335 DO ITER = 1,10
336 ! Derivative of fPS1 with respect to WEIGHT
337 U = fUN(J,1)*((ONE - MU)**2) + TWO*fH(1)*MU*WEIGHT*(ONE-MU) + fBI*(MU**2)
338 UPRIM = TWO*MU*(ONE-MU)*fH(1)
339 V = (ONE-MU)**2 + TWO*WEIGHT*MU*(ONE-MU) + MU**2
340 VPRIM = TWO*MU*(ONE-MU)
341 DfPS1_DWEIGHT = (UPRIM*V - U*VPRIM)/(MAX(V**2,EM20))
342 ! Updating the Weight
343 WEIGHT = WEIGHT - (PHI / DfPS1_DWEIGHT)
344 ! Updating the fPS1 value
345 fPS1 = (fUN(J,1)*((ONE - MU)**2) + TWO*fH(1)*MU*WEIGHT*(ONE-MU) + fBI*(MU**2))/
346 . ((ONE-MU)**2 + TWO*WEIGHT*MU*(ONE-MU) + MU**2)
347 ! Updating the PHI value
348 PHI = fPS1-fPS(J,1)
349 ENDDO
350 ! Filling the second coordinate of the plane-strain point
351 fPS(J,2) = (TWO*fH(2)*MU*WEIGHT*(ONE-MU) + fBI*(MU**2))/
352 . ((ONE-MU)**2 + TWO*WEIGHT*MU*(ONE-MU) + MU**2)
353 ENDDO
354c
355 ! Simplified Vegter lite yield locus
356 ELSE
357 CALL HM_GET_INTV('mat_refanglemax',NANGLE,IS_AVAILABLE,LSUBMODEL)
358 IF (NANGLE > 10) THEN
359 CALL ANCMSG(MSGID=1779,MSGTYPE=MSGERROR,
360 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,C1=TITR)
361 ENDIF
362.NOT. IF (ALLOCATED(fUN)) ALLOCATE(fUN(NANGLE,2))
363.NOT. IF (ALLOCATED(rLank)) ALLOCATE(rLank(NANGLE))
364.NOT. IF (ALLOCATED(WPS)) ALLOCATE(WPS(NANGLE))
365.NOT. IF (ALLOCATED(WSH)) ALLOCATE(WSH(NANGLE))
366 ! Loop over angles (must be equally distributed between 0 and pi/2)
367 DO J = 1, NANGLE
368 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_fun_theta' ,fUN(J,1),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
369 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_r_theta' ,rLank(J),J,IS_AVAILABLE, LSUBMODEL, UNITAB)
370 IF (rLank(J) == 0) rLank(J) = ONE
371 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_w_ps' ,WPS(J) ,J,IS_AVAILABLE, LSUBMODEL, UNITAB)
372 CALL HM_GET_FLOAT_ARRAY_INDEX('mat_w_sh' ,WSH(J) ,J,IS_AVAILABLE, LSUBMODEL, UNITAB)
373 ENDDO
374 DO J = 1, NANGLE
375 IF (WSH(J) /= WSH(NANGLE-(J-1))) THEN
376 CALL ANCMSG(MSGID=1800,MSGTYPE=MSGWARNING,
377 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,C1=TITR)
378 WSH(NANGLE-(J-1)) = WSH(J)
379 ENDIF
380 ENDDO
381 ENDIF
382c
383c---------------------
384c Default values
385c---------------------
386 ! Density
387 IF (RHOR == ZERO) RHOR = RHO0
388 IF (Ires == 0) Ires = 2
389 ! Poisson's ratio
390 IF (nu < zero .OR. nu >= half) THEN
391 CALL ancmsg(msgid=49,
392 . msgtype=msgerror,
393 . anmode=aninfo_blind_2,
394 . r1=nu,
395 . i1=mat_id,
396 . c1=titr)
397 ENDIF
398 ! Elasticity parameter
399 bulk = young / three / (one - two*nu)
400 a11 = young / (one-nu*nu)
401 g = young / two / (one + nu)
402 ! Default strain-rate computation
403 IF (ivflag == 0) ivflag = 2
404 IF (ivflag > 3) THEN
405 CALL ancmsg(msgid=1802,msgtype=msgerror,
406 . anmode=aninfo_blind_1,i1=mat_id,c1=titr,
407 . i2=ivflag)
408 ENDIF
409 ! No Strain-rate dependence
410 IF ((tab_yld == 0).AND.((deps0 == zero).OR.(dg0 == zero).OR.(sigstar == zero))) THEN
411 sigstar = zero
412 deps0 = one
413 dg0 = one
414 israte = 0
415 asrate = zero
416 ! Info message
417 CALL ancmsg(msgid=1801,msgtype=msginfo,
418 . anmode=aninfo_blind_1,i1=mat_id,c1=titr)
419 ELSE
420 ! Set strain-rate filtering
421 israte = 1
422 ! Default cutting frequency
423 IF ((asrate == zero).OR.(ivflag == 1)) asrate = 10000.0d0*unitab%FAC_T_WORK
424 ENDIF
425 ! Boltzmann constant
426 kboltz = 8.6173303*em5 ! eV/K (~J/K)
427 kboltz = kboltz/(unitab%FAC_M_WORK*unitab%FAC_L_WORK*unitab%FAC_L_WORK/(unitab%FAC_T_WORK*unitab%FAC_T_WORK))
428 ! Kinematic hardening factor
429 IF (fisokin > one) fisokin = one
430 IF (fisokin < zero) fisokin = zero
431 ! Default interpolation computation
432 IF (ismooth == 0) ismooth = 1
433 ! Tabulated function scale factor
434 IF (tab_yld == 0) THEN
435 xscale = zero
436 yscale = zero
437 tab_temp = 0
438 ENDIF
439 IF ((tab_temp > 0).AND.(tini == zero)) tini = 293.0d0
440 IF ((yscale == zero).AND.(tab_yld>0)) THEN
441 CALL hm_get_floatv_dim('MAT_Yscale' ,yscale_unit ,is_available, lsubmodel, unitab)
442 yscale = one * yscale_unit
443 ENDIF
444 IF ((xscale == zero).AND.(tab_yld>0)) THEN
445 CALL hm_get_floatv_dim('MAT_Xscale' ,xscale_unit ,is_available, lsubmodel, unitab)
446 xscale = one * xscale_unit
447 ENDIF
448c----------------------------------------------
449c Computation of yield criterion parameter
450c----------------------------------------------
451 IF (.NOT.ALLOCATED(hips)) ALLOCATE(hips(nangle,2))
452 IF (.NOT.ALLOCATED(hiun)) ALLOCATE(hiun(nangle,2))
453 IF (.NOT.ALLOCATED(hish)) ALLOCATE(hish(nangle,2))
454 IF (.NOT.ALLOCATED(theta)) ALLOCATE(theta(nangle))
455 IF (.NOT.ALLOCATED(theta_rad)) ALLOCATE(theta_rad(nangle))
456c
457 ! Computation of angles
458 DO j = 1, nangle
459 IF (nangle > 1) THEN
460 theta(j) = (j-1)*(90.0d0/(nangle-1))
461 theta_rad(j) = theta(j)*(pi/180.0d0)
462 ELSE
463 theta(j) = zero
464 theta_rad(j) = zero
465 ENDIF
466 ENDDO
467c
468 ! Loop over angles
469 DO j = 1, nangle
470c
471 ! hinge points for the classical corus-vegter yield locus(3 hinge points)
472 IF ((icrit == 1).OR.(icrit == 2).OR.(icrit == 3)) THEN
473 ! Computation of Hinge point between Bi-axial and Plane-strain point
474 ! Plane strain point and normal
475 a1 = fps(j,1)
476 a2 = fps(j,2)
477 n1 = one
478 n2 = zero
479 ! Biaxial point and normal
480 rhobi_theta = (rhobi + one) + (rhobi - one)*cos(two*theta_rad(j))
481 rhobi_theta = rhobi_theta / ((rhobi + one) - (rhobi - one)*cos(two*theta_rad(j)))
482 c1 = fbi
483 c2 = fbi
484 m1 = one
485 m2 = rhobi_theta
486 ! Hinge point
487 hips(j,1) = (m2*(n1*a1+n2*a2) - n2*(m1*c1+m2*c2))/(n1*m2 - m1*n2)
488 hips(j,2) = (n1*(m1*c1+m2*c2) - m1*(n1*a1+n2*a2))/(n1*m2 - m1*n2)
489c
490 ! Computation of Hinge point between Plane-strain and Uniaxial point
491 ! Uniaxial point and normal
492 a1 = fun(j,1)
493 fun(j,2) = zero
494 a2 = fun(j,2)
495 n1 = one
496 n2 = -rlank(j)/(rlank(j)+one)
497 ! Plane strain point and normal
498 c1 = fps(j,1)
499 c2 = fps(j,2)
500 m1 = one
501 m2 = zero
502 ! Hinge point
503 hiun(j,1) = (m2*(n1*a1+n2*a2) - n2*(m1*c1+m2*c2))/(n1*m2 - m1*n2)
504 hiun(j,2) = (n1*(m1*c1+m2*c2) - m1*(n1*a1+n2*a2))/(n1*m2 - m1*n2)
505c
506 ! If fPS2 is not measured or if Icrit = 2
507 IF ((icrit == 1).AND.(fps(j,2) == zero)) THEN
508 fps(j,2) = hiun(j,2) + half*(hips(j,2) - hiun(j,2))
509 ELSEIF (icrit == 2) THEN
510 fps(j,2) = hiun(j,2) + alps(j)*(hips(j,2) - hiun(j,2))
511 ENDIF
512c
513 ! Computation of Hinge point between Uniaxial point and Shear point
514 ! Shear point and normal
515 a1 = fsh(j,1)
516 fsh(j,2) = -fsh(nangle-(j-1),1)
517 a2 = fsh(j,2)
518 n1 = one
519 n2 = -one
520 ! Plane strain point and normal
521 c1 = fun(j,1)
522 c2 = zero
523 m1 = one
524 m2 = -rlank(j)/(rlank(j)+one)
525 ! Hinge point
526 hish(j,1) = (m2*(n1*a1+n2*a2) - n2*(m1*c1+m2*c2))/(n1*m2 - m1*n2)
527 hish(j,2) = (n1*(m1*c1+m2*c2) - m1*(n1*a1+n2*a2))/(n1*m2 - m1*n2)
528c
529 ! Hinge points for the simplified Corus-Vegter lite yield locus (2 hinge points)
530 ELSE
531 ! Computation of Hinge point between Bi-axial and Uniaxial point
532 ! Uniaxial point and normal
533 a1 = fun(j,1)
534 fun(j,2) = zero
535 a2 = fun(j,2)
536 n1 = one
537 n2 = -rlank(j)/(rlank(j)+one)
538 ! Biaxial point and normal
539 rhobi_theta = (rhobi + one) + (rhobi - one)*cos(two*theta_rad(j))
540 rhobi_theta = rhobi_theta / ((rhobi + one) - (rhobi - one)*cos(two*theta_rad(j)))
541 c1 = fbi
542 c2 = fbi
543 m1 = one
544 m2 = rhobi_theta
545 ! Hinge point
546 hips(j,1) = (m2*(n1*a1+n2*a2) - n2*(m1*c1+m2*c2))/(n1*m2 - m1*n2)
547 hips(j,2) = (n1*(m1*c1+m2*c2) - m1*(n1*a1+n2*a2))/(n1*m2 - m1*n2)
548c
549 ! Computation of Hinge point between Uniaxial point in tension and Uniaxial point in compression
550 ! Uniaxial point in compression
551 a1 = fun(j,2)
552 a2 = -fun(j,1)
553 n1 = rlank(j)/(rlank(j)+one)
554 n2 = -one
555 ! Uniaxial point in tension
556 c1 = fun(j,1)
557 c2 = fun(j,2)
558 m1 = one
559 m2 = -rlank(j)/(rlank(j)+one)
560 ! Hinge point
561 hish(j,1) = (m2*(n1*a1+n2*a2) - n2*(m1*c1+m2*c2))/(n1*m2 - m1*n2)
562 hish(j,2) = (n1*(m1*c1+m2*c2) - m1*(n1*a1+n2*a2))/(n1*m2 - m1*n2)
563 ENDIF
564 ENDDO
565c
566c----------------------------------------------
567c Computation of cosine interpolation factor
568c----------------------------------------------
569c
570 ! Allocation of the A-MATRIX and the Pivot vector
571 ALLOCATE (amat(nangle,nangle),ipiv(nangle))
572c
573 ! Filling the A-MATRIX
574 DO j = 1,nangle
575 DO i = 1,nangle
576 amat(j,i) = zero
577 DO k = 1,i
578 amat(j,i) = amat(j,i) + cos2(k,i)*(cos(two*theta_rad(j)))**(k-1)
579 ENDDO
580 ENDDO
581 ENDDO
582c
583 ! Classical Vegter yield locus and Standard Vegter yield locus
584 IF ((icrit == 1).OR.(icrit == 2).OR.(icrit == 3)) THEN
585c
586 ! Allocation of factors
587 ALLOCATE(q_fsh(nangle,2),q_fun(nangle,2),q_fps(nangle,2),
588 . q_hish(nangle,2),q_hiun(nangle,2),q_hips(nangle,2))
589c
590 ! Initialization of tables
591 q_fsh(1:nangle,1:2) = zero
592 q_fun(1:nangle,1:2) = zero
593 q_fps(1:nangle,1:2) = zero
594 q_hish(1:nangle,1:2) = zero
595 q_hiun(1:nangle,1:2) = zero
596 q_hips(1:nangle,1:2) = zero
597c
598 ! Filling the B vector with experimental points
599 ALLOCATE (bvec(nangle,12))
600 bvec(1:nangle,1:2) = fsh(1:nangle,1:2)
601 bvec(1:nangle,3:4) = fun(1:nangle,1:2)
602 bvec(1:nangle,5:6) = fps(1:nangle,1:2)
603 bvec(1:nangle,7:8) = hish(1:nangle,1:2)
604 bvec(1:nangle,9:10) = hiun(1:nangle,1:2)
605 bvec(1:nangle,11:12) = hips(1:nangle,1:2)
606c
607 ! Initialization of the Pivot vector
608 ipiv(1:nangle) = 0
609c
610 ! Solving the A-MATRIX * x = B vector system
611#ifndef WITHOUT_LINALG
612 CALL dgesv(nangle, 12, amat, nangle, ipiv, bvec, nangle, info)
613#else
614 WRITE(6,*) "Error: Blas/Lapack required for MAT110"
615#endif
616c
617 ! Recovering the solution
618 q_fsh(1:nangle,1:2) = bvec(1:nangle,1:2)
619 q_fun(1:nangle,1:2) = bvec(1:nangle,3:4)
620 q_fps(1:nangle,1:2) = bvec(1:nangle,5:6)
621 q_hish(1:nangle,1:2) = bvec(1:nangle,7:8)
622 q_hiun(1:nangle,1:2) = bvec(1:nangle,9:10)
623 q_hips(1:nangle,1:2) = bvec(1:nangle,11:12)
624c
625 ! Simplified Vegter Lite yield locus
626 ELSE
627c
628 ! Allocation of factors
629 ALLOCATE(q_fun(nangle,2),q_hish(nangle,2),q_hips(nangle,2),
630 . q_wsh(nangle),q_wps(nangle))
631c
632 ! Initialization of factors
633 q_fun(1:nangle,1:2) = zero
634 q_hish(1:nangle,1:2) = zero
635 q_hips(1:nangle,1:2) = zero
636 q_wsh(1:nangle) = zero
637 q_wps(1:nangle) = zero
638c
639 ! Filling the B vector with experimental points
640 ALLOCATE (bvec(nangle,8))
641 bvec(1:nangle,1:2) = fun(1:nangle,1:2)
642 bvec(1:nangle,3:4) = hish(1:nangle,1:2)
643 bvec(1:nangle,5:6) = hips(1:nangle,1:2)
644 bvec(1:nangle,7) = wsh(1:nangle)
645 bvec(1:nangle,8) = wps(1:nangle)
646c
647 ! Initialization of the Pivot vector
648 ipiv(1:nangle) = 0
649c
650 ! Solving the A-MATRIX * x = B vector system
651#ifndef WITHOUT_LINALG
652 CALL dgesv(nangle, 8, amat, nangle, ipiv, bvec, nangle, info)
653#else
654 WRITE(6,*) "Error: Blas/Lapack required for MAT110"
655#endif
656
657c
658 ! Recovering the solution
659 q_fun(1:nangle,1:2) = bvec(1:nangle,1:2)
660 q_hish(1:nangle,1:2) = bvec(1:nangle,3:4)
661 q_hips(1:nangle,1:2) = bvec(1:nangle,5:6)
662 q_wsh(1:nangle) = bvec(1:nangle,7)
663 q_wps(1:nangle) = bvec(1:nangle,8)
664c
665 ENDIF
666c-----------------------------------
667c Checking yield locus convexity
668c-----------------------------------
669c
670 ! Classical Vegter yield locus and Standard Yield locus
671 IF ((icrit == 1).OR.(icrit == 2).OR.(icrit == 3)) THEN
672
673 ! ---------------------------------------------------------
674 ! 1. Checking convexity between shear and uniaxial point
675 ! ---------------------------------------------------------
676c
677 ! Minimum eigen value of the Hessian matrix of the yield locus
678 mineig = infinity
679 mineigmu = zero
680 mineigthet = zero
681c
682 ! Initialization of the THETA angle
683 thet = zero
684c
685 ! Loop over the angles
686 DO j = 1,61
687c
688 ! Initialization of the MU parameter
689 mu = zero
690 ! Computing the reference and hinge point for a giving angle
691 ! Vector point
692 a1 = zero
693 a2 = zero
694 b1 = zero
695 b2 = zero
696 c1 = zero
697 c2 = zero
698 ! Initialization of the first derivative with respect to COS2THET
699 da_dcos2(1:2) = zero
700 db_dcos2(1:2) = zero
701 dc_dcos2(1:2) = zero
702 ! Initialization of the second derivative with respect to COS2THET
703 d2a_d2cos2(1:2) = zero
704 d2b_d2cos2(1:2) = zero
705 d2c_d2cos2(1:2) = zero
706 ! Loop over angles
707 DO i = 1,nangle
708 ! Computation of the reference points
709 DO k = 1,i
710 a1 = a1 + q_fsh(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
711 a2 = a2 + q_fsh(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
712 b1 = b1 + q_hish(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
713 b2 = b2 + q_hish(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
714 c1 = c1 + q_fun(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
715 c2 = c2 + q_fun(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
716 ENDDO
717 ENDDO
718 ! If anisotropic, compute derivatives with respect to COS2THET
719 IF (nangle > 1) THEN
720 ! Computation of their first derivative with respect to COS2THET
721 DO i = 2, nangle
722 DO k = 2,i
723 da_dcos2(1:2) = da_dcos2(1:2) + q_fsh(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
724 db_dcos2(1:2) = db_dcos2(1:2) + q_hish(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
725 dc_dcos2(1:2) = dc_dcos2(1:2) + q_fun(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
726 ENDDO
727 ENDDO
728 ! Computation of their second derivative with respect to COS2THET
729 IF (nangle > 2) THEN
730 DO i = 3, nangle
731 DO k = 3,i
732 d2a_d2cos2(1:2) = d2a_d2cos2(1:2) + q_fsh(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
733 d2b_d2cos2(1:2) = d2b_d2cos2(1:2) + q_hish(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
734 d2c_d2cos2(1:2) = d2c_d2cos2(1:2) + q_fun(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
735 ENDDO
736 ENDDO
737 ENDIF
738 ENDIF
739c
740 !----------------------------------------------
741 ! Computing the Hessian matrix and eigen values
742 !----------------------------------------------
743 DO i = 1,21
744c
745 ! Computation of F1 and F2
746 f1 = ((one-mu)**2)*a1 + two*mu*(one-mu)*b1 + (mu**2)*c1
747 f2 = ((one-mu)**2)*a2 + two*mu*(one-mu)*b2 + (mu**2)*c2
748 ! Derivatives of Fi with respect to MU
749 df1_dmu = -two*(one-mu)*a1 + two*(one-two*mu)*b1 + two*mu*c1
750 df2_dmu = -two*(one-mu)*a2 + two*(one-two*mu)*b2 + two*mu*c2
751 ! Second derivatives of Fi with respect to MU
752 d2f1_d2mu = two*(a1+c1-two*b1)
753 d2f2_d2mu = two*(a2+c2-two*b2)
754 ! Denominator and its derivative with respect to MU
755 denom = f1*df2_dmu - f2*df1_dmu
756 denom = sign(max(abs(denom),em20),denom)
757 ! Derivatives of Phi with respect to SIG1 and SIG2
758 dphi_dsig1 = df2_dmu/denom
759 dphi_dsig2 = -df1_dmu/denom
760 ! Derivatives of MU with respect to SIG1 and SIG2
761 dmu_dsig1 = -f2/denom
762 dmu_dsig2 = f1/denom
763c
764 ! Derivatives of Fi with respect to COS2
765 df1_dcos2 = ((one-mu)**2)*da_dcos2(1) + two*mu*(one-mu)*db_dcos2(1) + (mu**2)*dc_dcos2(1)
766 df2_dcos2 = ((one-mu)**2)*da_dcos2(2) + two*mu*(one-mu)*db_dcos2(2) + (mu**2)*dc_dcos2(2)
767 dphi_dcos2 = (df2_dcos2*df1_dmu - df1_dcos2*df2_dmu)/denom
768 dmu_dcos2 = (f2*df1_dcos2 - f1*df2_dcos2)/denom
769c
770 ! Second derivatives of Fi with respect to COS2
771 d2f1_d2mucos2 = -two*(one-mu)*da_dcos2(1) + two*(one-two*mu)*db_dcos2(1) + two*mu*dc_dcos2(1)
772 d2f2_d2mucos2 = -two*(one-mu)*da_dcos2(2) + two*(one-two*mu)*db_dcos2(2) + two*mu*dc_dcos2(2)
773 d2f1_d2cos2 = ((one-mu)**2)*d2a_d2cos2(1) + two*mu*(one-mu)*d2b_d2cos2(1) + (mu**2)*d2c_d2cos2(1)
774 d2f2_d2cos2 = ((one-mu)**2)*d2a_d2cos2(2) + two*mu*(one-mu)*d2b_d2cos2(2) + (mu**2)*d2c_d2cos2(2)
775c
776 ! Second derivatives of PHI with respect to SIG1 and SIG2 (Hessian matrix)
777 !-------------------------------------------------------------------------
778 ! FIRST LINE OF THE MATRIX
779 ! D2PHI_D2SIG1
780 d2phi_d2(1,1) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig1*dmu_dsig1
781 ! D2PHI_DSIG1SIG2
782 d2phi_d2(1,2) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig1*dmu_dsig2
783 ! D2PHI_DSIG1COS2
784 d2phi_d2(1,3) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
785 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig1
786 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig1)
787c
788 ! SECOND LINE OF THE MATRIX
789 ! D2PHI_DSIG2SIG1
790 d2phi_d2(2,1) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig2*dmu_dsig1
791 ! D2PHI_D2SIG2
792 d2phi_d2(2,2) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig2*dmu_dsig2
793 ! D2PHI_DSIG2COS2
794 d2phi_d2(2,3) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
795 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig2
796 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig2)
797c
798 ! THIRD LINE OF THE MATRIX
799 ! D2PHI_D2COS2SIG1
800 d2phi_d2(3,1) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
801 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig1
802 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig1)
803 ! D2PHI_D2COS2SIG2
804 d2phi_d2(3,2) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
805 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig2
806 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig2)
807 ! D2PHI_D2COS2
808 d2phi_d2(3,3) = (-one/denom)*(df2_dmu*(d2f1_d2cos2 + two*d2f1_d2mucos2*dmu_dcos2 + d2f1_d2mu*(dmu_dcos2**2))
809 . - df1_dmu*(d2f2_d2cos2 + two*d2f2_d2mucos2*dmu_dcos2 + d2f2_d2mu*(dmu_dcos2**2))
810 . + two*(df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dcos2)
811c
812 ! Filling the T matrix
813 t(1,1) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig1-dphi_dsig2) +
814 . dphi_dcos2*(-three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
815 t(1,2) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig2 - dphi_dsig1) +
816 . dphi_dcos2*(three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
817 t(1,3) = (one/(f1-f2))*(sin(two*thet)*cos(two*thet))*(dphi_dsig2 - dphi_dsig1) +
818 . dphi_dcos2*(one/(f1-f2)**2)*(four*sin(two*thet)*(cos(two*thet))**2 -
819 . two*(sin(two*thet)**3))
820 t(2,1) = t(1,2)
821 t(2,2) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig1-dphi_dsig2) +
822 . dphi_dcos2*(-three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
823 t(2,3) = (one/(f1-f2))*(sin(two*thet)*cos(two*thet))*(dphi_dsig1-dphi_dsig2) +
824 . dphi_dcos2*(one/(f1-f2)**2)*(-four*sin(two*thet)*(cos(two*thet))**2 +
825 . two*(sin(two*thet)**3))
826 t(3,1) = t(1,3)
827 t(3,2) = t(2,3)
828 t(3,3) = (two/(f1-f2))*((cos(two*thet))**2)*(dphi_dsig1-dphi_dsig2)+
829 . dphi_dcos2*(one/(f1-f2)**2)*(eight*(sin(two*thet)**2)*cos(two*thet) -
830 . four*(cos(two*thet)**3))
831
832 ! Filling the rotation D matrix
833 d(1,1) = half*(one+cos(two*thet))
834 d(1,2) = half*(one-cos(two*thet))
835 d(1,3) = sin(two*thet)
836 d(2,1) = half*(one-cos(two*thet))
837 d(2,2) = half*(one+cos(two*thet))
838 d(2,3) = -sin(two*thet)
839 d(3,1) = (sin(two*thet)**2)/(f1-f2)
840 d(3,2) = -(sin(two*thet)**2)/(f1-f2)
841 d(3,3) = -(two*sin(two*thet)*cos(two*thet))/(f1-f2)
842
843 ! Assembling the Hessian Matrix Dt*DPHI_DSIG*D + T
844 d2phi_d2 = matmul(d2phi_d2,d)
845 d2phi_d2 = matmul(transpose(d),d2phi_d2)
846 d2phi_d2(1:3,1:3) = d2phi_d2(1:3,1:3) + t(1:3,1:3)
847c
848 ! Computation of eigen values
849 wr = zero
850 wi = zero
851 vl = zero
852 vr = zero
853 work = zero
854#ifndef WITHOUT_LINALG
855 CALL dgeev('N','N',3,d2phi_d2,3,wr,wi,vl,3,vr,3,work,102,info)
856#else
857 WRITE(6,*) "Error: Blas/Lapack required for MAT110"
858#endif
859c
860 ! Storing the minimal value
861 IF (minval(wr)<mineig) THEN
862 mineig = minval(wr)
863 mineigmu = mu
864 mineigthet = thet*(180.0d0/pi)
865 ENDIF
866 ! Updating the MU parameter
867 mu = mu + one/twenty
868 ENDDO
869 ! Updating the theta angle
870 thet = thet + (pi/two)/60.0d0
871 ENDDO
872c
873 ! Convexity check and error message
874 IF (mineig<-tol) THEN
875 CALL ancmsg(msgid=1803,msgtype=msgerror,
876 . anmode=aninfo_blind_1,i1=mat_id,c1=titr,
877 . r1=mineig,r2=mineigmu,r3=mineigthet)
878 ENDIF
879c
880 ! -------------------------------------------------------------
881 ! 2. Checking convexity between uniaxial and plane strain point
882 ! -------------------------------------------------------------
883c
884 ! Minimum eigen value of the Hessian matrix of the yield locus
885 mineig = infinity
886 mineigmu = zero
887 mineigthet = zero
888c
889 ! Initialization of the THETA angle
890 thet = zero
891c
892 ! Loop over the angles
893 DO j = 1,61
894c
895 ! Initialization of the MU parameter
896 mu = zero
897 ! Computing the reference and hinge point for a giving angle
898 ! Vector point
899 a1 = zero
900 a2 = zero
901 b1 = zero
902 b2 = zero
903 c1 = zero
904 c2 = zero
905 ! Initialization of the first derivative with respect to COS2THET
906 da_dcos2(1:2) = zero
907 db_dcos2(1:2) = zero
908 dc_dcos2(1:2) = zero
909 ! Initialization of the second derivative with respect to COS2THET
910 d2a_d2cos2(1:2) = zero
911 d2b_d2cos2(1:2) = zero
912 d2c_d2cos2(1:2) = zero
913 ! Loop over angles
914 DO i = 1,nangle
915 ! Computation of the reference points
916 DO k = 1,i
917 a1 = a1 + q_fun(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
918 a2 = a2 + q_fun(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
919 b1 = b1 + q_hiun(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
920 b2 = b2 + q_hiun(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
921 c1 = c1 + q_fps(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
922 c2 = c2 + q_fps(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
923 ENDDO
924 ENDDO
925 ! If anisotropic, compute derivatives with respect to COS2THET
926 IF (nangle > 1) THEN
927 ! Computation of their first derivative with respect to COS2THET
928 DO i = 2, nangle
929 DO k = 2,i
930 da_dcos2(1:2) = da_dcos2(1:2) + q_fun(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
931 db_dcos2(1:2) = db_dcos2(1:2) + q_hiun(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
932 dc_dcos2(1:2) = dc_dcos2(1:2) + q_fps(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
933 ENDDO
934 ENDDO
935 ! Computation of their second derivative with respect to COS2THET
936 IF (nangle > 2) THEN
937 DO i = 3, nangle
938 DO k = 3,i
939 d2a_d2cos2(1:2) = d2a_d2cos2(1:2) + q_fun(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
940 d2b_d2cos2(1:2) = d2b_d2cos2(1:2) + q_hiun(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
941 d2c_d2cos2(1:2) = d2c_d2cos2(1:2) + q_fps(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
942 ENDDO
943 ENDDO
944 ENDIF
945 ENDIF
946c
947 !----------------------------------------------
948 ! Computing the Hessian matrix and eigen values
949 !----------------------------------------------
950 DO i = 1,21
951c
952 ! Computation of F1 and F2
953 f1 = ((one-mu)**2)*a1 + two*mu*(one-mu)*b1 + (mu**2)*c1
954 f2 = ((one-mu)**2)*a2 + two*mu*(one-mu)*b2 + (mu**2)*c2
955 ! Derivatives of Fi with respect to MU
956 df1_dmu = -two*(one-mu)*a1 + two*(one-two*mu)*b1 + two*mu*c1
957 df2_dmu = -two*(one-mu)*a2 + two*(one-two*mu)*b2 + two*mu*c2
958 ! Second derivatives of Fi with respect to MU
959 d2f1_d2mu = two*(a1+c1-two*b1)
960 d2f2_d2mu = two*(a2+c2-two*b2)
961 ! Denominator and its derivative with respect to MU
962 denom = f1*df2_dmu - f2*df1_dmu
963 denom = sign(max(abs(denom),em20),denom)
964 ! Derivatives of Phi with respect to SIG1 and SIG2
965 dphi_dsig1 = df2_dmu/denom
966 dphi_dsig2 = -df1_dmu/denom
967 ! Derivatives of MU with respect to SIG1 and SIG2
968 dmu_dsig1 = -f2/denom
969 dmu_dsig2 = f1/denom
970c
971 ! Derivatives of Fi with respect to COS2
972 df1_dcos2 = ((one-mu)**2)*da_dcos2(1) + two*mu*(one-mu)*db_dcos2(1) + (mu**2)*dc_dcos2(1)
973 df2_dcos2 = ((one-mu)**2)*da_dcos2(2) + two*mu*(one-mu)*db_dcos2(2) + (mu**2)*dc_dcos2(2)
974 dphi_dcos2 = (df2_dcos2*df1_dmu - df1_dcos2*df2_dmu)/denom
975 dmu_dcos2 = (f2*df1_dcos2 - f1*df2_dcos2)/denom
976c
977 ! Second derivatives of Fi with respect to COS2
978 d2f1_d2mucos2 = -two*(one-mu)*da_dcos2(1) + two*(one-two*mu)*db_dcos2(1) + two*mu*dc_dcos2(1)
979 d2f2_d2mucos2 = -two*(one-mu)*da_dcos2(2) + two*(one-two*mu)*db_dcos2(2) + two*mu*dc_dcos2(2)
980 d2f1_d2cos2 = ((one-mu)**2)*d2a_d2cos2(1) + two*mu*(one-mu)*d2b_d2cos2(1) + (mu**2)*d2c_d2cos2(1)
981 d2f2_d2cos2 = ((one-mu)**2)*d2a_d2cos2(2) + two*mu*(one-mu)*d2b_d2cos2(2) + (mu**2)*d2c_d2cos2(2)
982c
983 ! Second derivatives of PHI with respect to SIG1 and SIG2 (Hessian matrix)
984 !-------------------------------------------------------------------------
985 ! FIRST LINE OF THE MATRIX
986 ! D2PHI_D2SIG1
987 d2phi_d2(1,1) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig1*dmu_dsig1
988 ! D2PHI_DSIG1SIG2
989 d2phi_d2(1,2) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig1*dmu_dsig2
990 ! D2PHI_DSIG1COS2
991 d2phi_d2(1,3) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
992 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig1
993 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig1)
994c
995 ! SECOND LINE OF THE MATRIX
996 ! D2PHI_DSIG2SIG1
997 d2phi_d2(2,1) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig2*dmu_dsig1
998 ! D2PHI_D2SIG2
999 d2phi_d2(2,2) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig2*dmu_dsig2
1000 ! D2PHI_DSIG2COS2
1001 d2phi_d2(2,3) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1002 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig2
1003 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig2)
1004c
1005 ! THIRD LINE OF THE MATRIX
1006 ! D2PHI_D2COS2SIG1
1007 d2phi_d2(3,1) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1008 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig1
1009 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig1)
1010 ! D2PHI_D2COS2SIG2
1011 d2phi_d2(3,2) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1012 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig2
1013 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig2)
1014 ! D2PHI_D2COS2
1015 d2phi_d2(3,3) = (-one/denom)*(df2_dmu*(d2f1_d2cos2 + two*d2f1_d2mucos2*dmu_dcos2 + d2f1_d2mu*(dmu_dcos2**2))
1016 . - df1_dmu*(d2f2_d2cos2 + two*d2f2_d2mucos2*dmu_dcos2 + d2f2_d2mu*(dmu_dcos2**2))
1017 . + two*(df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dcos2)
1018c
1019 ! Filling the T matrix
1020 t(1,1) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig1-dphi_dsig2) +
1021 . dphi_dcos2*(-three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1022 t(1,2) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig2 - dphi_dsig1) +
1023 . dphi_dcos2*(three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1024 t(1,3) = (one/(f1-f2))*(sin(two*thet)*cos(two*thet))*(dphi_dsig2 - dphi_dsig1) +
1025 . dphi_dcos2*(one/(f1-f2)**2)*(four*sin(two*thet)*(cos(two*thet))**2 -
1026 . two*(sin(two*thet)**3))
1027 t(2,1) = t(1,2)
1028 t(2,2) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig1-dphi_dsig2) +
1029 . dphi_dcos2*(-three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1030 t(2,3) = (one/(f1-f2))*(sin(two*thet)*cos(two*thet))*(dphi_dsig1-dphi_dsig2) +
1031 . dphi_dcos2*(one/(f1-f2)**2)*(-four*sin(two*thet)*(cos(two*thet))**2 +
1032 . two*(sin(two*thet)**3))
1033 t(3,1) = t(1,3)
1034 t(3,2) = t(2,3)
1035 t(3,3) = (two/(f1-f2))*((cos(two*thet))**2)*(dphi_dsig1-dphi_dsig2)+
1036 . dphi_dcos2*(one/(f1-f2)**2)*(eight*(sin(two*thet)**2)*cos(two*thet) -
1037 . four*(cos(two*thet)**3))
1038c
1039 ! Filling the rotation D matrix
1040 d(1,1) = half*(one+cos(two*thet))
1041 d(1,2) = half*(one-cos(two*thet))
1042 d(1,3) = sin(two*thet)
1043 d(2,1) = half*(one-cos(two*thet))
1044 d(2,2) = half*(one+cos(two*thet))
1045 d(2,3) = -sin(two*thet)
1046 d(3,1) = (sin(two*thet)**2)/(f1-f2)
1047 d(3,2) = -(sin(two*thet)**2)/(f1-f2)
1048 d(3,3) = -(two*sin(two*thet)*cos(two*thet))/(f1-f2)
1049c
1050 ! Assembling the Hessian Matrix Dt*DPHI_DSIG*D + T
1051 d2phi_d2 = matmul(d2phi_d2,d)
1052 d2phi_d2 = matmul(transpose(d),d2phi_d2)
1053 d2phi_d2(1:3,1:3) = d2phi_d2(1:3,1:3) + t(1:3,1:3)
1054c
1055 ! Computation of eigen values
1056 wr = zero
1057 wi = zero
1058 vl = zero
1059 vr = zero
1060 work = zero
1061
1062
1063
1064#ifndef WITHOUT_LINALG
1065 CALL dgeev('N','N',3,d2phi_d2,3,wr,wi,vl,3,vr,3,work,102,info)
1066#else
1067 WRITE(6,*) "Error: Blas/Lapack required for MAT110"
1068#endif
1069c
1070 ! Storing the minimal value
1071 IF (minval(wr)<mineig) THEN
1072 mineig = minval(wr)
1073 mineigmu = mu
1074 mineigthet = thet*(180.0d0/pi)
1075 ENDIF
1076 ! Updating the MU parameter
1077 mu = mu + one/twenty
1078 ENDDO
1079 ! Updating the theta angle
1080 thet = thet + (pi/two)/60.0d0
1081 ENDDO
1082c
1083 ! Convexity check and error message
1084 IF (mineig<-tol) THEN
1085 CALL ancmsg(msgid=1804,msgtype=msgerror,
1086 . anmode=aninfo_blind_1,i1=mat_id,c1=titr,
1087 . r1=mineig,r2=mineigmu,r3=mineigthet)
1088 ENDIF
1089c
1090 ! -------------------------------------------------------------
1091 ! 3. Checking convexity between plane-strain and biaxial point
1092 ! -------------------------------------------------------------
1093c
1094 ! Minimum eigen value of the Hessian matrix of the yield locus
1095 mineig = infinity
1096 mineigmu = zero
1097 mineigthet = zero
1098c
1099 ! Initialization of the THETA angle
1100 thet = zero
1101c
1102 ! Loop over the angles
1103 DO j = 1,61
1104c
1105 ! Initialization of the MU parameter
1106 mu = zero
1107 ! Computing the reference and hinge point for a giving angle
1108 ! Vector point
1109 a1 = zero
1110 a2 = zero
1111 b1 = zero
1112 b2 = zero
1113 c1 = fbi
1114 c2 = fbi
1115 ! Initialization of the first derivative with respect to COS2THET
1116 da_dcos2(1:2) = zero
1117 db_dcos2(1:2) = zero
1118 dc_dcos2(1:2) = zero
1119 ! Initialization of the second derivative with respect to COS2THET
1120 d2a_d2cos2(1:2) = zero
1121 d2b_d2cos2(1:2) = zero
1122 d2c_d2cos2(1:2) = zero
1123 ! Loop over angles
1124 DO i = 1,nangle
1125 ! Computation of the reference points
1126 DO k = 1,i
1127 a1 = a1 + q_fps(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
1128 a2 = a2 + q_fps(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
1129 b1 = b1 + q_hips(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
1130 b2 = b2 + q_hips(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
1131 ENDDO
1132 ENDDO
1133 ! If anisotropic, compute derivatives with respect to COS2THET
1134 IF (nangle > 1) THEN
1135 ! Computation of their first derivative with respect to COS2THET
1136 DO i = 2, nangle
1137 DO k = 2,i
1138 da_dcos2(1:2) = da_dcos2(1:2) + q_fps(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
1139 db_dcos2(1:2) = db_dcos2(1:2) + q_hips(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
1140 ENDDO
1141 ENDDO
1142 ! Computation of their second derivative with respect to COS2THET
1143 IF (nangle > 2) THEN
1144 DO i = 3, nangle
1145 DO k = 3,i
1146 d2a_d2cos2(1:2) = d2a_d2cos2(1:2) + q_fps(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
1147 d2b_d2cos2(1:2) = d2b_d2cos2(1:2) + q_hips(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
1148 ENDDO
1149 ENDDO
1150 ENDIF
1151 ENDIF
1152c
1153 !----------------------------------------------
1154 ! Computing the Hessian matrix and eigen values
1155 !----------------------------------------------
1156 DO i = 1,20
1157c
1158 ! Computation of F1 and F2
1159 f1 = ((one-mu)**2)*a1 + two*mu*(one-mu)*b1 + (mu**2)*c1
1160 f2 = ((one-mu)**2)*a2 + two*mu*(one-mu)*b2 + (mu**2)*c2
1161 ! Derivatives of Fi with respect to MU
1162 df1_dmu = -two*(one-mu)*a1 + two*(one-two*mu)*b1 + two*mu*c1
1163 df2_dmu = -two*(one-mu)*a2 + two*(one-two*mu)*b2 + two*mu*c2
1164 ! Second derivatives of Fi with respect to MU
1165 d2f1_d2mu = two*(a1+c1-two*b1)
1166 d2f2_d2mu = two*(a2+c2-two*b2)
1167 ! Denominator and its derivative with respect to MU
1168 denom = f1*df2_dmu - f2*df1_dmu
1169 denom = sign(max(abs(denom),em20),denom)
1170 ! Derivatives of Phi with respect to SIG1 and SIG2
1171 dphi_dsig1 = df2_dmu/denom
1172 dphi_dsig2 = -df1_dmu/denom
1173 ! Derivatives of MU with respect to SIG1 and SIG2
1174 dmu_dsig1 = -f2/denom
1175 dmu_dsig2 = f1/denom
1176c
1177 ! Derivatives of Fi with respect to COS2
1178 df1_dcos2 = ((one-mu)**2)*da_dcos2(1) + two*mu*(one-mu)*db_dcos2(1) + (mu**2)*dc_dcos2(1)
1179 df2_dcos2 = ((one-mu)**2)*da_dcos2(2) + two*mu*(one-mu)*db_dcos2(2) + (mu**2)*dc_dcos2(2)
1180 dphi_dcos2 = (df2_dcos2*df1_dmu - df1_dcos2*df2_dmu)/denom
1181 dmu_dcos2 = (f2*df1_dcos2 - f1*df2_dcos2)/denom
1182c
1183 ! Second derivatives of Fi with respect to COS2
1184 d2f1_d2mucos2 = -two*(one-mu)*da_dcos2(1) + two*(one-two*mu)*db_dcos2(1) + two*mu*dc_dcos2(1)
1185 d2f2_d2mucos2 = -two*(one-mu)*da_dcos2(2) + two*(one-two*mu)*db_dcos2(2) + two*mu*dc_dcos2(2)
1186 d2f1_d2cos2 = ((one-mu)**2)*d2a_d2cos2(1) + two*mu*(one-mu)*d2b_d2cos2(1) + (mu**2)*d2c_d2cos2(1)
1187 d2f2_d2cos2 = ((one-mu)**2)*d2a_d2cos2(2) + two*mu*(one-mu)*d2b_d2cos2(2) + (mu**2)*d2c_d2cos2(2)
1188c
1189 ! Second derivatives of PHI with respect to SIG1 and SIG2 (Hessian matrix)
1190 !-------------------------------------------------------------------------
1191 ! FIRST LINE OF THE MATRIX
1192 ! D2PHI_D2SIG1
1193 d2phi_d2(1,1) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig1*dmu_dsig1
1194 ! D2PHI_DSIG1SIG2
1195 d2phi_d2(1,2) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig1*dmu_dsig2
1196 ! D2PHI_DSIG1COS2
1197 d2phi_d2(1,3) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1198 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig1
1199 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig1)
1200c
1201 ! SECOND LINE OF THE MATRIX
1202 ! D2PHI_DSIG2SIG1
1203 d2phi_d2(2,1) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig2*dmu_dsig1
1204 ! D2PHI_D2SIG2
1205 d2phi_d2(2,2) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig2*dmu_dsig2
1206 ! D2PHI_DSIG2COS2
1207 d2phi_d2(2,3) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1208 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig2
1209 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig2)
1210c
1211 ! THIRD LINE OF THE MATRIX
1212 ! D2PHI_D2COS2SIG1
1213 d2phi_d2(3,1) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1214 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig1
1215 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig1)
1216 ! D2PHI_D2COS2SIG2
1217 d2phi_d2(3,2) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1218 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig2
1219 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig2)
1220 ! D2PHI_D2COS2
1221 d2phi_d2(3,3) = (-one/denom)*(df2_dmu*(d2f1_d2cos2 + two*d2f1_d2mucos2*dmu_dcos2 + d2f1_d2mu*(dmu_dcos2**2))
1222 . - df1_dmu*(d2f2_d2cos2 + two*d2f2_d2mucos2*dmu_dcos2 + d2f2_d2mu*(dmu_dcos2**2))
1223 . + two*(df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dcos2)
1224c
1225 ! Filling the T matrix
1226 t(1,1) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig1-dphi_dsig2) +
1227 . dphi_dcos2*(-three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1228 t(1,2) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig2 - dphi_dsig1) +
1229 . dphi_dcos2*(three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1230 t(1,3) = (one/(f1-f2))*(sin(two*thet)*cos(two*thet))*(dphi_dsig2 - dphi_dsig1) +
1231 . dphi_dcos2*(one/(f1-f2)**2)*(four*sin(two*thet)*(cos(two*thet))**2 -
1232 . two*(sin(two*thet)**3))
1233 t(2,1) = t(1,2)
1234 t(2,2) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig1-dphi_dsig2) +
1235 . dphi_dcos2*(-three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1236 t(2,3) = (one/(f1-f2))*(sin(two*thet)*cos(two*thet))*(dphi_dsig1-dphi_dsig2) +
1237 . dphi_dcos2*(one/(f1-f2)**2)*(-four*sin(two*thet)*(cos(two*thet))**2 +
1238 . two*(sin(two*thet)**3))
1239 t(3,1) = t(1,3)
1240 t(3,2) = t(2,3)
1241 t(3,3) = (two/(f1-f2))*((cos(two*thet))**2)*(dphi_dsig1-dphi_dsig2)+
1242 . dphi_dcos2*(one/(f1-f2)**2)*(eight*(sin(two*thet)**2)*cos(two*thet) -
1243 . four*(cos(two*thet)**3))
1244
1245 ! Filling the rotation D matrix
1246 d(1,1) = half*(one+cos(two*thet))
1247 d(1,2) = half*(one-cos(two*thet))
1248 d(1,3) = sin(two*thet)
1249 d(2,1) = half*(one-cos(two*thet))
1250 d(2,2) = half*(one+cos(two*thet))
1251 d(2,3) = -sin(two*thet)
1252 d(3,1) = (sin(two*thet)**2)/(f1-f2)
1253 d(3,2) = -(sin(two*thet)**2)/(f1-f2)
1254 d(3,3) = -(two*sin(two*thet)*cos(two*thet))/(f1-f2)
1255
1256 ! Assembling the Hessian Matrix Dt*DPHI_DSIG*D + T
1257 d2phi_d2 = matmul(d2phi_d2,d)
1258 d2phi_d2 = matmul(transpose(d),d2phi_d2)
1259 d2phi_d2(1:3,1:3) = d2phi_d2(1:3,1:3) + t(1:3,1:3)
1260c
1261 ! Computation of eigen values
1262 wr = zero
1263 wi = zero
1264 vl = zero
1265 vr = zero
1266 work = zero
1267#ifndef WITHOUT_LINALG
1268 CALL dgeev('N','N',3,d2phi_d2,3,wr,wi,vl,3,vr,3,work,102,info)
1269#else
1270 WRITE(6,*) "Error: Blas/Lapack required for MAT110"
1271#endif
1272c
1273 ! Storing the minimal value
1274 IF (minval(wr)<mineig) THEN
1275 mineig = minval(wr)
1276 mineigmu = mu
1277 mineigthet = thet*(180.0d0/pi)
1278 ENDIF
1279c
1280 ! Updating the MU parameter
1281 mu = mu + one/twenty
1282 ENDDO
1283 ! Updating the theta angle
1284 thet = thet + (pi/two)/60.0d0
1285 ENDDO
1286c
1287 ! Convexity check and error message
1288 IF (mineig<-tol) THEN
1289 CALL ancmsg(msgid=1805,msgtype=msgerror,
1290 . anmode=aninfo_blind_1,i1=mat_id,c1=titr,
1291 . r1=mineig,r2=mineigmu,r3=mineigthet)
1292 ENDIF
1293c
1294 ! Simplified Vegter Lite yield locus
1295 ELSE
1296c
1297 ! ------------------------------------------------------------------------
1298 ! 1. Checking convexity between uniaxial points in compression and tension
1299 ! ------------------------------------------------------------------------
1300c
1301 ! Minimum eigen value of the Hessian matrix of the yield locus
1302 mineig = infinity
1303 mineigmu = zero
1304 mineigthet = zero
1305c
1306 ! Initialization of the THETA angle
1307 thet = zero
1308c
1309 ! Loop over the angles
1310 DO j = 1,61
1311c
1312 ! Initialization of the MU parameter
1313 mu = zero
1314 ! Computing the reference and hinge point for a giving angle
1315 ! Vector point
1316 a1 = zero
1317 a2 = zero
1318 b1 = zero
1319 b2 = zero
1320 c1 = zero
1321 c2 = zero
1322 weight = zero
1323 ! Initialization of the first derivative with respect to COS2THET
1324 da_dcos2(1:2) = zero
1325 db_dcos2(1:2) = zero
1326 dc_dcos2(1:2) = zero
1327 dweight_dcos2 = zero
1328 ! Initialization of the second derivative with respect to COS2THET
1329 d2a_d2cos2(1:2) = zero
1330 d2b_d2cos2(1:2) = zero
1331 d2c_d2cos2(1:2) = zero
1332 d2weight_d2cos2 = zero
1333 ! Loop over angles
1334 DO i = 1,nangle
1335 ! Computation of the reference points
1336 DO k = 1,i
1337 b1 = b1 + q_hish(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
1338 b2 = b2 + q_hish(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
1339 c1 = c1 + q_fun(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
1340 c2 = c2 + q_fun(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
1341 weight = weight + q_wsh(i)*cos2(k,i)*(cos(two*thet))**(k-1)
1342 ENDDO
1343 ENDDO
1344 a1 = c2
1345 a2 = -c1
1346 ! If anisotropic, compute derivatives with respect to COS2THET
1347 IF (nangle > 1) THEN
1348 ! Computation of their first derivative with respect to COS2THET
1349 DO i = 2, nangle
1350 DO k = 2,i
1351 db_dcos2(1:2) = db_dcos2(1:2) + q_hish(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
1352 dc_dcos2(1:2) = dc_dcos2(1:2) + q_fun(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
1353 dweight_dcos2 = dweight_dcos2 + q_wsh(i)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
1354 ENDDO
1355 ENDDO
1356 da_dcos2(1) = dc_dcos2(2)
1357 da_dcos2(2) = -dc_dcos2(1)
1358 ! Computation of their second derivative with respect to COS2THET
1359 IF (nangle > 2) THEN
1360 DO i = 3, nangle
1361 DO k = 3,i
1362 d2b_d2cos2(1:2) = d2b_d2cos2(1:2) + q_hish(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
1363 d2c_d2cos2(1:2) = d2c_d2cos2(1:2) + q_fun(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
1364 d2weight_d2cos2 = d2weight_d2cos2 + q_wsh(i)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
1365 ENDDO
1366 ENDDO
1367 d2a_d2cos2(1) = d2c_d2cos2(2)
1368 d2a_d2cos2(2) = -d2c_d2cos2(1)
1369 ENDIF
1370 ENDIF
1371c
1372 !----------------------------------------------
1373 ! Computing the Hessian matrix and eigen values
1374 !----------------------------------------------
1375 DO i = 1,21
1376c
1377 ! Computation of F1 and F2
1378 f1 = ((one-mu)**2)*a1 + two*mu*weight*(one-mu)*b1 + (mu**2)*c1
1379 f1 = f1/(((one-mu)**2) + two*mu*weight*(one-mu) + (mu**2))
1380 f2 = ((one-mu)**2)*a2 + two*mu*weight*(one-mu)*b2 + (mu**2)*c2
1381 f2 = f2/(((one-mu)**2) + two*mu*weight*(one-mu) + (mu**2))
1382 ! Derivatives of Fi with respect to MU
1383 u = a1*((one - mu)**2) + two*b1*mu*weight*(one-mu) + c1*(mu**2)
1384 uprim = -two*(one-mu)*a1 + two*weight*b1*(one-two*mu) + two*mu*c1
1385 upprim = two*a1 - four*weight*b1 + two*c1
1386 v = (one-mu)**2 + two*weight*mu*(one-mu) + mu**2
1387 vprim = -two*(one-mu) + two*weight*(one-two*mu) + two*mu
1388 vpprim = two - four*weight + two
1389 df1_dmu = (uprim*v - u*vprim)/max((v**2),em20)
1390 d2f1_d2mu = ((upprim*v - u*vpprim)*(v**2) -
1391 . (uprim*v - u*vprim)*(two*v*vprim))/max((v**4),em20)
1392 u = a2*((one - mu)**2) + two*b2*mu*weight*(one-mu) + c2*(mu**2)
1393 uprim = -two*(one-mu)*a2 + two*weight*b2*(one-two*mu) + two*mu*c2
1394 upprim = two*a2 - four*weight*b2 + two*c2
1395 df2_dmu = (uprim*v - u*vprim)/max((v**2),em20)
1396 d2f2_d2mu = ((upprim*v - u*vpprim)*(v**2) -
1397 . (uprim*v - u*vprim)*(two*v*vprim))/max((v**4),em20)
1398 ! Denominator and its derivative with respect to MU
1399 denom = f1*df2_dmu - f2*df1_dmu
1400 denom = sign(max(abs(denom),em20),denom)
1401 ! Derivatives of Phi with respect to SIG1 and SIG2
1402 dphi_dsig1 = df2_dmu/denom
1403 dphi_dsig2 = -df1_dmu/denom
1404 ! Derivatives of MU with respect to SIG1 and SIG2
1405 dmu_dsig1 = -f2/denom
1406 dmu_dsig2 = f1/denom
1407c
1408 ! Derivatives of Fi with respect to COS2
1409 u = a1*((one - mu)**2) + two*b1*mu*weight*(one-mu) + c1*(mu**2)
1410 uprim = da_dcos2(1)*((one - mu)**2) + two*mu*(one-mu)*(dweight_dcos2*b1+
1411 . weight*db_dcos2(1)) + dc_dcos2(1)*(mu**2)
1412 upprim = d2a_d2cos2(1)*((one - mu)**2) + two*mu*(one-mu)*(d2weight_d2cos2*b1+
1413 . two*dweight_dcos2*db_dcos2(1) + weight*d2b_d2cos2(1)) +
1414 . d2c_d2cos2(1)*(mu**2)
1415 du_dmu = -two*(one-mu)*a1 + two*weight*b1*(one-two*mu) + two*mu*c1
1416 duprim_dmu = -two*(one-mu)*da_dcos2(1) + two*(one-two*mu)*(dweight_dcos2*b1+
1417 . weight*db_dcos2(1)) + two*mu*dc_dcos2(1)
1418 v = (one-mu)**2 + two*weight*mu*(one-mu) + mu**2
1419 dv_dmu = -two*(one-mu) + two*weight*(one-two*mu) + two*mu
1420 vprim = two*mu*(one-mu)*dweight_dcos2
1421 dvprim_dmu = two*(one-two*mu)*dweight_dcos2
1422 vpprim = two*mu*(one-mu)*d2weight_d2cos2
1423 df1_dcos2 = (uprim*v - u*vprim)/max((v**2),em20)
1424 d2f1_d2cos2 = ((upprim*v - u*vpprim)*(v**2) -
1425 . (uprim*v - u*vprim)*(two*v*vprim))/max((v**4),em20)
1426 d2f1_d2mucos2 = ((duprim_dmu*v + uprim*dv_dmu - du_dmu*vprim - u*dvprim_dmu)*(v**2) -
1427 . (uprim*v - u*vprim)*(two*v*dv_dmu))/max((v**4),em20)
1428c
1429 u = a2*((one - mu)**2) + two*b2*mu*weight*(one-mu) + c2*(mu**2)
1430 uprim = da_dcos2(2)*((one - mu)**2) + two*mu*(one-mu)*(dweight_dcos2*b2+
1431 . weight*db_dcos2(2)) + dc_dcos2(2)*(mu**2)
1432 upprim = d2a_d2cos2(2)*((one - mu)**2) + two*mu*(one-mu)*(d2weight_d2cos2*b2+
1433 . two*dweight_dcos2*db_dcos2(2) + weight*d2b_d2cos2(2)) +
1434 . d2c_d2cos2(2)*(mu**2)
1435 du_dmu = -two*(one-mu)*a2 + two*weight*b2*(one-two*mu) + two*mu*c2
1436 duprim_dmu = -two*(one-mu)*da_dcos2(2) + two*(one-two*mu)*(dweight_dcos2*b2+
1437 . weight*db_dcos2(2)) + two*mu*dc_dcos2(2)
1438 df2_dcos2 = (uprim*v - u*vprim)/max((v**2),em20)
1439 d2f2_d2cos2 = ((upprim*v - u*vpprim)*(v**2) -
1440 . (uprim*v - u*vprim)*(two*v*vprim))/max((v**4),em20)
1441 d2f2_d2mucos2 = ((duprim_dmu*v + uprim*dv_dmu - du_dmu*vprim - u*dvprim_dmu)*(v**2) -
1442 . (uprim*v - u*vprim)*(two*v*dv_dmu))/max((v**4),em20)
1443c
1444 dphi_dcos2 = (df2_dcos2*df1_dmu - df1_dcos2*df2_dmu)/denom
1445 dmu_dcos2 = (f2*df1_dcos2 - f1*df2_dcos2)/denom
1446c
1447 ! Second derivatives of PHI with respect to SIG1 and SIG2 (Hessian matrix)
1448 !-------------------------------------------------------------------------
1449 ! FIRST LINE OF THE MATRIX
1450 ! D2PHI_D2SIG1
1451 d2phi_d2(1,1) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig1*dmu_dsig1
1452 ! D2PHI_DSIG1SIG2
1453 d2phi_d2(1,2) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig1*dmu_dsig2
1454 ! D2PHI_DSIG1COS2
1455 d2phi_d2(1,3) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1456 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig1
1457 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig1)
1458c
1459 ! SECOND LINE OF THE MATRIX
1460 ! D2PHI_DSIG2SIG1
1461 d2phi_d2(2,1) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig2*dmu_dsig1
1462 ! D2PHI_D2SIG2
1463 d2phi_d2(2,2) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig2*dmu_dsig2
1464 ! D2PHI_DSIG2COS2
1465 d2phi_d2(2,3) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1466 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig2
1467 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig2)
1468c
1469 ! THIRD LINE OF THE MATRIX
1470 ! D2PHI_D2COS2SIG1
1471 d2phi_d2(3,1) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1472 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig1
1473 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig1)
1474 ! D2PHI_D2COS2SIG2
1475 d2phi_d2(3,2) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1476 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig2
1477 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig2)
1478 ! D2PHI_D2COS2
1479 d2phi_d2(3,3) = (-one/denom)*(df2_dmu*(d2f1_d2cos2 + two*d2f1_d2mucos2*dmu_dcos2 + d2f1_d2mu*(dmu_dcos2**2))
1480 . - df1_dmu*(d2f2_d2cos2 + two*d2f2_d2mucos2*dmu_dcos2 + d2f2_d2mu*(dmu_dcos2**2))
1481 . + two*(df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dcos2)
1482c
1483 ! Filling the T matrix
1484 t(1,1) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig1-dphi_dsig2) +
1485 . dphi_dcos2*(-three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1486 t(1,2) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig2 - dphi_dsig1) +
1487 . dphi_dcos2*(three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1488 t(1,3) = (one/(f1-f2))*(sin(two*thet)*cos(two*thet))*(dphi_dsig2 - dphi_dsig1) +
1489 . dphi_dcos2*(one/(f1-f2)**2)*(four*sin(two*thet)*(cos(two*thet))**2 -
1490 . two*(sin(two*thet)**3))
1491 t(2,1) = t(1,2)
1492 t(2,2) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig1-dphi_dsig2) +
1493 . dphi_dcos2*(-three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1494 t(2,3) = (one/(f1-f2))*(sin(two*thet)*cos(two*thet))*(dphi_dsig1-dphi_dsig2) +
1495 . dphi_dcos2*(one/(f1-f2)**2)*(-four*sin(two*thet)*(cos(two*thet))**2 +
1496 . two*(sin(two*thet)**3))
1497 t(3,1) = t(1,3)
1498 t(3,2) = t(2,3)
1499 t(3,3) = (two/(f1-f2))*((cos(two*thet))**2)*(dphi_dsig1-dphi_dsig2)+
1500 . dphi_dcos2*(one/(f1-f2)**2)*(eight*(sin(two*thet)**2)*cos(two*thet) -
1501 . four*(cos(two*thet)**3))
1502
1503 ! Filling the rotation D matrix
1504 d(1,1) = half*(one+cos(two*thet))
1505 d(1,2) = half*(one-cos(two*thet))
1506 d(1,3) = sin(two*thet)
1507 d(2,1) = half*(one-cos(two*thet))
1508 d(2,2) = half*(one+cos(two*thet))
1509 d(2,3) = -sin(two*thet)
1510 d(3,1) = (sin(two*thet)**2)/(f1-f2)
1511 d(3,2) = -(sin(two*thet)**2)/(f1-f2)
1512 d(3,3) = -(two*sin(two*thet)*cos(two*thet))/(f1-f2)
1513
1514 ! Assembling the Hessian Matrix Dt*DPHI_DSIG*D + T
1515 d2phi_d2 = matmul(d2phi_d2,d)
1516 d2phi_d2 = matmul(transpose(d),d2phi_d2)
1517 d2phi_d2(1:3,1:3) = d2phi_d2(1:3,1:3) + t(1:3,1:3)
1518c
1519 ! Computation of eigen values
1520 wr = zero
1521 wi = zero
1522 vl = zero
1523 vr = zero
1524 work = zero
1525#ifndef WITHOUT_LINALG
1526 CALL dgeev('N','N',3,d2phi_d2,3,wr,wi,vl,3,vr,3,work,102,info)
1527c
1528#else
1529 WRITE(6,*) "Error: Blas/Lapack required for MAT110"
1530#endif
1531
1532 ! Storing the minimal value
1533 IF (minval(wr)<mineig) THEN
1534 mineig = minval(wr)
1535 mineigmu = mu
1536 mineigthet = thet*(180.0d0/pi)
1537 ENDIF
1538 ! Updating the MU parameter
1539 mu = mu + one/twenty
1540 ENDDO
1541 ! Updating the theta angle
1542 thet = thet + (pi/two)/60.0d0
1543 ENDDO
1544c
1545 ! Convexity check and error message
1546 IF (mineig<-tol) THEN
1547 CALL ancmsg(msgid=1806,msgtype=msgerror,
1548 . anmode=aninfo_blind_1,i1=mat_id,c1=titr,
1549 . r1=mineig,r2=mineigmu,r3=mineigthet)
1550 ENDIF
1551c
1552 ! ------------------------------------------------------------------------
1553 ! 2. Checking convexity between uniaxial tension and equi-biaxial points
1554 ! ------------------------------------------------------------------------
1555c
1556 ! Minimum eigen value of the Hessian matrix of the yield locus
1557 mineig = infinity
1558 mineigmu = zero
1559 mineigthet = zero
1560c
1561 ! Initialization of the THETA angle
1562 thet = zero
1563c
1564 ! Loop over the angles
1565 DO j = 1,61
1566c
1567 ! Initialization of the MU parameter
1568 mu = zero
1569 ! Computing the reference and hinge point for a giving angle
1570 ! Vector point
1571 a1 = zero
1572 a2 = zero
1573 b1 = zero
1574 b2 = zero
1575 c1 = zero
1576 c2 = zero
1577 weight = zero
1578 ! Initialization of the first derivative with respect to COS2THET
1579 da_dcos2(1:2) = zero
1580 db_dcos2(1:2) = zero
1581 dc_dcos2(1:2) = zero
1582 dweight_dcos2 = zero
1583 ! Initialization of the second derivative with respect to COS2THET
1584 d2a_d2cos2(1:2) = zero
1585 d2b_d2cos2(1:2) = zero
1586 d2c_d2cos2(1:2) = zero
1587 d2weight_d2cos2 = zero
1588 ! Loop over angles
1589 DO i = 1,nangle
1590 ! Computation of the reference points
1591 DO k = 1,i
1592 a1 = a1 + q_fun(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
1593 a2 = a2 + q_fun(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
1594 b1 = b1 + q_hips(i,1)*cos2(k,i)*(cos(two*thet))**(k-1)
1595 b2 = b2 + q_hips(i,2)*cos2(k,i)*(cos(two*thet))**(k-1)
1596 weight = weight + q_wps(i)*cos2(k,i)*(cos(two*thet))**(k-1)
1597 ENDDO
1598 ENDDO
1599 c1 = fbi
1600 c2 = fbi
1601 ! If anisotropic, compute derivatives with respect to COS2THET
1602 IF (nangle > 1) THEN
1603 ! Computation of their first derivative with respect to COS2THET
1604 DO i = 2, nangle
1605 DO k = 2,i
1606 da_dcos2(1:2) = da_dcos2(1:2) + q_fun(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
1607 db_dcos2(1:2) = db_dcos2(1:2) + q_hips(i,1:2)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
1608 dweight_dcos2 = dweight_dcos2 + q_wps(i)*(k-1)*cos2(k,i)*(cos(two*thet))**(k-2)
1609 ENDDO
1610 ENDDO
1611 ! Computation of their second derivative with respect to COS2THET
1612 IF (nangle > 2) THEN
1613 DO i = 3, nangle
1614 DO k = 3,i
1615 d2a_d2cos2(1:2) = d2a_d2cos2(1:2) + q_fun(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
1616 d2b_d2cos2(1:2) = d2b_d2cos2(1:2) + q_hips(i,1:2)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
1617 d2weight_d2cos2 = d2weight_d2cos2 + q_wps(i)*(k-1)*(k-2)*cos2(k,i)*(cos(two*thet))**(k-3)
1618 ENDDO
1619 ENDDO
1620 ENDIF
1621 ENDIF
1622c
1623 !----------------------------------------------
1624 ! Computing the Hessian matrix and eigen values
1625 !----------------------------------------------
1626 DO i = 1,20
1627c
1628 ! Computation of F1 and F2
1629 f1 = ((one-mu)**2)*a1 + two*mu*weight*(one-mu)*b1 + (mu**2)*c1
1630 f1 = f1/(((one-mu)**2) + two*mu*weight*(one-mu) + (mu**2))
1631 f2 = ((one-mu)**2)*a2 + two*mu*weight*(one-mu)*b2 + (mu**2)*c2
1632 f2 = f2/(((one-mu)**2) + two*mu*weight*(one-mu) + (mu**2))
1633 ! Derivatives of Fi with respect to MU
1634 u = a1*((one - mu)**2) + two*b1*mu*weight*(one-mu) + c1*(mu**2)
1635 uprim = -two*(one-mu)*a1 + two*weight*b1*(one-two*mu) + two*mu*c1
1636 upprim = two*a1 - four*weight*b1 + two*c1
1637 v = (one-mu)**2 + two*weight*mu*(one-mu) + mu**2
1638 vprim = -two*(one-mu) + two*weight*(one-two*mu) + two*mu
1639 vpprim = two - four*weight + two
1640 df1_dmu = (uprim*v - u*vprim)/max((v**2),em20)
1641 d2f1_d2mu = ((upprim*v - u*vpprim)*(v**2) -
1642 . (uprim*v - u*vprim)*(two*v*vprim))/max((v**4),em20)
1643 u = a2*((one - mu)**2) + two*b2*mu*weight*(one-mu) + c2*(mu**2)
1644 uprim = -two*(one-mu)*a2 + two*weight*b2*(one-two*mu) + two*mu*c2
1645 upprim = two*a2 - four*weight*b2 + two*c2
1646 df2_dmu = (uprim*v - u*vprim)/max((v**2),em20)
1647 d2f2_d2mu = ((upprim*v - u*vpprim)*(v**2) -
1648 . (uprim*v - u*vprim)*(two*v*vprim))/max((v**4),em20)
1649 ! Denominator and its derivative with respect to MU
1650 denom = f1*df2_dmu - f2*df1_dmu
1651 denom = sign(max(abs(denom),em20),denom)
1652 ! Derivatives of Phi with respect to SIG1 and SIG2
1653 dphi_dsig1 = df2_dmu/denom
1654 dphi_dsig2 = -df1_dmu/denom
1655 ! Derivatives of MU with respect to SIG1 and SIG2
1656 dmu_dsig1 = -f2/denom
1657 dmu_dsig2 = f1/denom
1658c
1659 ! Derivatives of Fi with respect to COS2
1660 u = a1*((one - mu)**2) + two*b1*mu*weight*(one-mu) + c1*(mu**2)
1661 uprim = da_dcos2(1)*((one - mu)**2) + two*mu*(one-mu)*(dweight_dcos2*b1+
1662 . weight*db_dcos2(1)) + dc_dcos2(1)*(mu**2)
1663 upprim = d2a_d2cos2(1)*((one - mu)**2) + two*mu*(one-mu)*(d2weight_d2cos2*b1+
1664 . two*dweight_dcos2*db_dcos2(1) + weight*d2b_d2cos2(1)) +
1665 . d2c_d2cos2(1)*(mu**2)
1666 du_dmu = -two*(one-mu)*a1 + two*weight*b1*(one-two*mu) + two*mu*c1
1667 duprim_dmu = -two*(one-mu)*da_dcos2(1) + two*(one-two*mu)*(dweight_dcos2*b1+
1668 . weight*db_dcos2(1)) + two*mu*dc_dcos2(1)
1669 v = (one-mu)**2 + two*weight*mu*(one-mu) + mu**2
1670 dv_dmu = -two*(one-mu) + two*weight*(one-two*mu) + two*mu
1671 vprim = two*mu*(one-mu)*dweight_dcos2
1672 dvprim_dmu = two*(one-two*mu)*dweight_dcos2
1673 vpprim = two*mu*(one-mu)*d2weight_d2cos2
1674 df1_dcos2 = (uprim*v - u*vprim)/max((v**2),em20)
1675 d2f1_d2cos2 = ((upprim*v - u*vpprim)*(v**2) -
1676 . (uprim*v - u*vprim)*(two*v*vprim))/max((v**4),em20)
1677 d2f1_d2mucos2 = ((duprim_dmu*v + uprim*dv_dmu - du_dmu*vprim - u*dvprim_dmu)*(v**2) -
1678 . (uprim*v - u*vprim)*(two*v*dv_dmu))/max((v**4),em20)
1679c
1680 u = a2*((one - mu)**2) + two*b2*mu*weight*(one-mu) + c2*(mu**2)
1681 uprim = da_dcos2(2)*((one - mu)**2) + two*mu*(one-mu)*(dweight_dcos2*b2+
1682 . weight*db_dcos2(2)) + dc_dcos2(2)*(mu**2)
1683 upprim = d2a_d2cos2(2)*((one - mu)**2) + two*mu*(one-mu)*(d2weight_d2cos2*b2+
1684 . two*dweight_dcos2*db_dcos2(2) + weight*d2b_d2cos2(2)) +
1685 . d2c_d2cos2(2)*(mu**2)
1686 du_dmu = -two*(one-mu)*a2 + two*weight*b2*(one-two*mu) + two*mu*c2
1687 duprim_dmu = -two*(one-mu)*da_dcos2(2) + two*(one-two*mu)*(dweight_dcos2*b2+
1688 . weight*db_dcos2(2)) + two*mu*dc_dcos2(2)
1689 df2_dcos2 = (uprim*v - u*vprim)/max((v**2),em20)
1690 d2f2_d2cos2 = ((upprim*v - u*vpprim)*(v**2) -
1691 . (uprim*v - u*vprim)*(two*v*vprim))/max((v**4),em20)
1692 d2f2_d2mucos2 = ((duprim_dmu*v + uprim*dv_dmu - du_dmu*vprim - u*dvprim_dmu)*(v**2) -
1693 . (uprim*v - u*vprim)*(two*v*dv_dmu))/max((v**4),em20)
1694c
1695 dphi_dcos2 = (df2_dcos2*df1_dmu - df1_dcos2*df2_dmu)/denom
1696 dmu_dcos2 = (f2*df1_dcos2 - f1*df2_dcos2)/denom
1697c
1698 ! Second derivatives of PHI with respect to SIG1 and SIG2 (Hessian matrix)
1699 !-------------------------------------------------------------------------
1700 ! FIRST LINE OF THE MATRIX
1701 ! D2PHI_D2SIG1
1702 d2phi_d2(1,1) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig1*dmu_dsig1
1703 ! D2PHI_DSIG1SIG2
1704 d2phi_d2(1,2) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig1*dmu_dsig2
1705 ! D2PHI_DSIG1COS2
1706 d2phi_d2(1,3) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1707 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig1
1708 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig1)
1709c
1710 ! SECOND LINE OF THE MATRIX
1711 ! D2PHI_DSIG2SIG1
1712 d2phi_d2(2,1) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig2*dmu_dsig1
1713 ! D2PHI_D2SIG2
1714 d2phi_d2(2,2) = (-one/denom)*(df2_dmu*d2f1_d2mu-df1_dmu*d2f2_d2mu)*dmu_dsig2*dmu_dsig2
1715 ! D2PHI_DSIG2COS2
1716 d2phi_d2(2,3) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1717 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig2
1718 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig2)
1719c
1720 ! THIRD LINE OF THE MATRIX
1721 ! D2PHI_D2COS2SIG1
1722 d2phi_d2(3,1) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1723 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig1
1724 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig1)
1725 ! D2PHI_D2COS2SIG2
1726 d2phi_d2(3,2) = (-one/denom)*((df2_dmu*(d2f1_d2mucos2+d2f1_d2mu*dmu_dcos2) -
1727 . df1_dmu*(d2f2_d2mucos2+d2f2_d2mu*dmu_dcos2))*dmu_dsig2
1728 . + (df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dsig2)
1729 ! D2PHI_D2COS2
1730 d2phi_d2(3,3) = (-one/denom)*(df2_dmu*(d2f1_d2cos2 + two*d2f1_d2mucos2*dmu_dcos2 + d2f1_d2mu*(dmu_dcos2**2))
1731 . - df1_dmu*(d2f2_d2cos2 + two*d2f2_d2mucos2*dmu_dcos2 + d2f2_d2mu*(dmu_dcos2**2))
1732 . + two*(df2_dmu*df1_dcos2 - df1_dmu*df2_dcos2)*dphi_dcos2)
1733c
1734 ! Filling the T matrix
1735 t(1,1) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig1-dphi_dsig2) +
1736 . dphi_dcos2*(-three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1737 t(1,2) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig2 - dphi_dsig1) +
1738 . dphi_dcos2*(three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1739 t(1,3) = (one/(f1-f2))*(sin(two*thet)*cos(two*thet))*(dphi_dsig2 - dphi_dsig1) +
1740 . dphi_dcos2*(one/(f1-f2)**2)*(four*sin(two*thet)*(cos(two*thet))**2 -
1741 . two*(sin(two*thet)**3))
1742 t(2,1) = t(1,2)
1743 t(2,2) = (half/(f1-f2))*((sin(two*thet))**2)*(dphi_dsig1-dphi_dsig2) +
1744 . dphi_dcos2*(-three/(f1-f2)**2)*cos(two*thet)*(sin(two*thet))**2
1745 t(2,3) = (one/(f1-f2))*(sin(two*thet)*cos(two*thet))*(dphi_dsig1-dphi_dsig2) +
1746 . dphi_dcos2*(one/(f1-f2)**2)*(-four*sin(two*thet)*(cos(two*thet))**2 +
1747 . two*(sin(two*thet)**3))
1748 t(3,1) = t(1,3)
1749 t(3,2) = t(2,3)
1750 t(3,3) = (two/(f1-f2))*((cos(two*thet))**2)*(dphi_dsig1-dphi_dsig2)+
1751 . dphi_dcos2*(one/(f1-f2)**2)*(eight*(sin(two*thet)**2)*cos(two*thet) -
1752 . four*(cos(two*thet)**3))
1753
1754 ! Filling the rotation D matrix
1755 d(1,1) = half*(one+cos(two*thet))
1756 d(1,2) = half*(one-cos(two*thet))
1757 d(1,3) = sin(two*thet)
1758 d(2,1) = half*(one-cos(two*thet))
1759 d(2,2) = half*(one+cos(two*thet))
1760 d(2,3) = -sin(two*thet)
1761 d(3,1) = (sin(two*thet)**2)/(f1-f2)
1762 d(3,2) = -(sin(two*thet)**2)/(f1-f2)
1763 d(3,3) = -(two*sin(two*thet)*cos(two*thet))/(f1-f2)
1764
1765 ! Assembling the Hessian Matrix Dt*DPHI_DSIG*D + T
1766 d2phi_d2 = matmul(d2phi_d2,d)
1767 d2phi_d2 = matmul(transpose(d),d2phi_d2)
1768 d2phi_d2(1:3,1:3) = d2phi_d2(1:3,1:3) + t(1:3,1:3)
1769c
1770 ! Computation of eigen values
1771 wr = zero
1772 wi = zero
1773 vl = zero
1774 vr = zero
1775 work = zero
1776#ifndef WITHOUT_LINALG
1777 CALL dgeev('N','n',3,D2PHI_D2,3,WR,WI,VL,3,VR,3,WORK,102,INFO)
1778#else
1779 WRITE(6,*) "Error: Blas/Lapack required for MAT110"
1780#endif
1781
1782c
1783 ! Storing the minimal value
1784 IF (MINVAL(WR)<MINEIG) THEN
1785 MINEIG = MINVAL(WR)
1786 MINEIGMU = MU
1787 MINEIGTHET = THET*(180.0D0/PI)
1788 ENDIF
1789 ! Updating the MU parameter
1790 MU = MU + ONE/TWENTY
1791 ENDDO
1792 ! Updating the theta angle
1793 THET = THET + (PI/TWO)/60.0D0
1794 ENDDO
1795c
1796 ! Convexity check and error message
1797 IF (MINEIG<-TOL) THEN
1798 CALL ANCMSG(MSGID=1807,MSGTYPE=MSGERROR,
1799 . ANMODE=ANINFO_BLIND_1,I1=MAT_ID,C1=TITR,
1800 . R1=MINEIG,R2=MINEIGMU,R3=MINEIGTHET)
1801 ENDIF
1802c
1803 ENDIF
1804c
1805c--------------------------
1806c Filling buffer tables
1807c--------------------------
1808 ! Number of material parameter
1809 NUPARAM = 29
1810.OR..OR. IF ((Icrit == 1)(Icrit == 2)(Icrit == 3)) THEN
1811 NUPARAM = NUPARAM + 17*NANGLE
1812 ELSE
1813 NUPARAM = NUPARAM + 11*NANGLE
1814 ENDIF
1815 ! Number of function and temp variable
1816 IF (TAB_YLD > 0) THEN
1817 NTABL = 1
1818 ITABLE(1) = TAB_YLD ! Yield function table = f(epsp,epsdot)
1819 NVARTMP = 1
1820 IF (TAB_TEMP > 0) THEN
1821 NTABL = 2
1822 ITABLE(2) = TAB_TEMP ! Temperature scale factor table = f(epsp,T)
1823 NVARTMP = 3
1824 ENDIF
1825 ELSE
1826 NTABL = 0
1827 NVARTMP = 0
1828 ENDIF
1829 ! Number of user variable
1830 NUVAR = 1
1831 IF (Ires == 1) THEN
1832 NUVAR = 2
1833 ENDIF
1834c
1835 ! Material parameters
1836 UPARAM(1) = YOUNG
1837 UPARAM(2) = NU
1838 UPARAM(3) = G
1839 UPARAM(4) = TWO*G
1840 UPARAM(5) = NU/(ONE-NU)
1841 UPARAM(6) = A11
1842 UPARAM(7) = A11*NU
1843 UPARAM(8) = YLD0
1844 UPARAM(9) = DSIGM
1845 UPARAM(10) = BETA
1846 UPARAM(11) = OMEGA
1847 UPARAM(12) = NEXP
1848 UPARAM(13) = EPS0
1849 UPARAM(14) = SIGSTAR
1850 UPARAM(15) = DG0
1851 UPARAM(16) = Deps0
1852 UPARAM(17) = MEXP
1853 UPARAM(18) = KBOLTZ
1854 IF (XSCALE /= ZERO) THEN
1855 UPARAM(19) = ONE/XSCALE
1856 ELSE
1857 UPARAM(19) = ZERO
1858 ENDIF
1859 UPARAM(20) = YSCALE
1860 UPARAM(21) = FISOKIN
1861 UPARAM(22) = ASRATE
1862 UPARAM(23) = Ivflag
1863 UPARAM(24) = TINI
1864 UPARAM(25) = fBI
1865 UPARAM(26) = NANGLE
1866 UPARAM(27) = Icrit
1867 UPARAM(28) = Ismooth
1868 UPARAM(29) = Ires
1869.OR..OR. IF ((Icrit == 1)(Icrit == 2)(Icrit == 3)) THEN
1870 DO J = 1,NANGLE
1871 UPARAM(30+17*(J-1)) = HIPS(J,1)
1872 UPARAM(31+17*(J-1)) = HIPS(J,2)
1873 UPARAM(32+17*(J-1)) = HIUN(J,1)
1874 UPARAM(33+17*(J-1)) = HIUN(J,2)
1875 UPARAM(34+17*(J-1)) = HISH(J,1)
1876 UPARAM(35+17*(J-1)) = HISH(J,2)
1877 UPARAM(36+17*(J-1)) = Q_fSH(J,1)
1878 UPARAM(37+17*(J-1)) = Q_fSH(J,2)
1879 UPARAM(38+17*(J-1)) = Q_fUN(J,1)
1880 UPARAM(39+17*(J-1)) = Q_fPS(J,1)
1881 UPARAM(40+17*(J-1)) = Q_fPS(J,2)
1882 UPARAM(41+17*(J-1)) = Q_HISH(J,1)
1883 UPARAM(42+17*(J-1)) = Q_HISH(J,2)
1884 UPARAM(43+17*(J-1)) = Q_HIUN(J,1)
1885 UPARAM(44+17*(J-1)) = Q_HIUN(J,2)
1886 UPARAM(45+17*(J-1)) = Q_HIPS(J,1)
1887 UPARAM(46+17*(J-1)) = Q_HIPS(J,2)
1888 ENDDO
1889 ELSE
1890 DO J = 1,NANGLE
1891 UPARAM(30+11*(J-1)) = HIPS(J,1)
1892 UPARAM(31+11*(J-1)) = HIPS(J,2)
1893 UPARAM(32+11*(J-1)) = HISH(J,1)
1894 UPARAM(33+11*(J-1)) = HISH(J,2)
1895 UPARAM(34+11*(J-1)) = Q_fUN(J,1)
1896 UPARAM(35+11*(J-1)) = Q_HISH(J,1)
1897 UPARAM(36+11*(J-1)) = Q_HISH(J,2)
1898 UPARAM(37+11*(J-1)) = Q_HIPS(J,1)
1899 UPARAM(38+11*(J-1)) = Q_HIPS(J,2)
1900 UPARAM(39+11*(J-1)) = Q_WSH(J)
1901 UPARAM(40+11*(J-1)) = Q_WPS(J)
1902 ENDDO
1903 ENDIF
1904c
1905 ! PARMAT table
1906 PARMAT(1) = BULK
1907 PARMAT(2) = YOUNG
1908 PARMAT(3) = NU
1909 PARMAT(4) = ISRATE
1910 PARMAT(5) = ASRATE
1911c
1912 ! PM table
1913 PM(1) = RHO0
1914 PM(89) = RHO0
1915 PM(27) = SQRT(A11/RHO0) ! sound speed estimation
1916 PM(100)= BULK
1917c
1918 ! MTAG variable activation
1919 MTAG%G_PLA = 1
1920 MTAG%L_PLA = 1
1921 MTAG%L_EPSD = 1
1922 MTAG%G_EPSD = 1
1923 MTAG%L_SEQ = 1
1924 MTAG%G_SEQ = 1
1925 MTAG%L_TEMP = 1
1926 MTAG%G_TEMP = 1
1927 IF (FISOKIN > ZERO) THEN
1928 MTAG%L_SIGA = 3
1929 ENDIF
1930c
1931 CALL INIT_MAT_KEYWORD(MATPARAM,"ELASTO_PLASTIC")
1932 CALL INIT_MAT_KEYWORD(MATPARAM,"INCREMENTAL" )
1933 CALL INIT_MAT_KEYWORD(MATPARAM,"LARGE_STRAIN" )
1934 CALL INIT_MAT_KEYWORD(MATPARAM,"ORTHOTROPIC" )
1935c
1936 ! Properties compatibility
1937 CALL INIT_MAT_KEYWORD(MATPARAM,"SHELL_ORTHOTROPIC")
1938c
1939c--------------------------
1940c Parameters printout
1941c--------------------------
1942 WRITE(IOUT,1000) TRIM(TITR),MAT_ID,ILAW
1943 IF (Icrit == 1) THEN
1944 WRITE(IOUT,1100)
1945 ELSEIF (Icrit == 2) THEN
1946 WRITE(IOUT,1125)
1947 ELSEIF (Icrit == 3) THEN
1948 WRITE(IOUT,1135)
1949 ELSE
1950 WRITE(IOUT,1150)
1951 ENDIF
1952 IF (IS_ENCRYPTED) THEN
1953 WRITE(IOUT,'(5x,a,//)')'confidential data'
1954 ELSE
1955 WRITE(IOUT,1200) RHO0
1956 WRITE(IOUT,1300) YOUNG,NU
1957 WRITE(IOUT,1350) Ires
1958 IF (ITABLE(1)>0) THEN
1959 WRITE(IOUT,1500) ITABLE(1),XSCALE,YSCALE,Ismooth,EPS0,FISOKIN,
1960 . Ivflag,ITABLE(2),TINI
1961 ELSE
1962 WRITE(IOUT,1400) YLD0,DSIGM,BETA,OMEGA,NEXP,EPS0,FISOKIN
1963 IF (SIGSTAR>ZERO) THEN
1964 WRITE(IOUT,1600) SIGSTAR,DG0,Deps0,MEXP,TINI,ASRATE,Ivflag
1965 ENDIF
1966 ENDIF
1967.OR..OR. IF ((Icrit == 1)(Icrit == 2)(Icrit == 3)) THEN
1968 DO J = 1,NANGLE
1969 WRITE(IOUT,1700) THETA(J),fSH(J,1),fSH(J,2),HISH(J,1),HISH(J,2),
1970 . fUN(J,1),fUN(J,2),rLank(J),HIUN(J,1),HIUN(J,2),
1971 . fPS(J,1),fPS(J,2),HIPS(J,1),HIPS(J,2),fBI,fBI
1972 ENDDO
1973 ELSE
1974 DO J = 1,NANGLE
1975 WRITE(IOUT,1750) THETA(J),HISH(J,1),HISH(J,2),WSH(J),
1976 . fUN(J,1),fUN(J,2),rLank(J),HIPS(J,1),
1977 . HIPS(J,2),WPS(J),fBI,fBI
1978 ENDDO
1979 ENDIF
1980 ENDIF
1981c-----------------------------------------------------------------------
1982 1000 FORMAT(/
1983 & 5X,A,/,
1984 & 5X,'material number. . . . . . . . . . . . =',I10/,
1985 & 5X,'material law . . . . . . . . . . . . . =',I10/)
1986 1100 FORMAT
1987 &(5X,'material model : vegter elasto-plastic',/,
1988 & 5X,'--------------------------------------',/)
1989 1125 FORMAT
1990 &(5X,'material model : vegter standard elasto-plastic',/,
1991 & 5X,'-----------------------------------------------',/)
1992 1135 FORMAT
1993 &(5X,'material model : vegter 2017 elasto-plastic',/,
1994 & 5X,'-------------------------------------------',/)
1995 1150 FORMAT
1996 &(5X,'material model : simplified vegter lite elasto-plastic',/,
1997 & 5X,'------------------------------------------------------',/)
1998 1200 FORMAT(
1999 & 5X,'initial density . . . . . . . . . . . .=',1PG20.13/)
2000 1300 FORMAT(
2001 & 5X,'young modulus . . . . . . . . . . . . .=',1PG20.13/
2002 & 5X,'poisson ratio . . . . . . . . . . . . .=',1PG20.13/)
2003 1350 FORMAT(
2004 & 5X,'RETURN mapping algorithm flag . . . . .=',I3/
2005 & 5X,' ires=1 nice explicit'/
2006 & 5X,' ires=2 newton-iteration implicit(cutting plane)'/)
2007 1400 FORMAT(
2008 & 5X,'initial yield stress . . . . . . . . .=',1PG20.13/
2009 & 5X,'stress hardening increment. . . . . . .=',1PG20.13/
2010 & 5X,'small strain hardening PARAMETER . . .=',1PG20.13/
2011 & 5X,'large strain hardening PARAMETER . . .=',1PG20.13/
2012 & 5X,'hardening exponent . . . . . . . . . .=',1PG20.13/
2013 & 5X,'initial plastic strain . . . . . . . .=',1PG20.13/
2014 & 5X,'kinematic hardening factor . . . . . .=',1PG20.13/)
2015 1500 FORMAT(
2016 & 5X,'tabulated yield - strain rate table id.=',I10/
2017 & 5X,'tabulated yield x factor . . . . . . .=',1PG20.13/
2018 & 5X,'tabulated yield y factor . . . . . . .=',1PG20.13/
2019 & 5X,'table interpolation flag . . . . . . .=',I10/
2020 & 5X,' ismooth=1 linear interpolation'/
2021 & 5X,' ismooth=2 logarithmic interpolation base 10'/
2022 & 5X,' ismooth=3 logarithmic interpolation base n'/
2023 & 5X,'initial plastic strain . . . . . . . .=',1PG20.13/
2024 & 5X,'kinematic hardening factor . . . . . .=',1PG20.13/
2025 & 5X,'strain-rate choice flag . . . . . . . .=',I10/
2026 & 5X,' vp=1 equivalent plastic strain rate'/
2027 & 5X,' vp=2 total strain rate(default)'/
2028 & 5X,' vp=3 deviatoric strain rate'/
2029 & 5X,'tabulated yield - temperature table id .=',I10/
2030 & 5X,'initial(reference) temperature. . . . .=',1PG20.13/)
2031 1600 FORMAT(
2032 & 5X,'limit dynamic flow stress . . . . . . .=',1PG20.13/
2033 & 5X,'maximum activation enthalpy . . . . . .=',1PG20.13/
2034 & 5X,'thermal activation limit strain-rate .=',1PG20.13/
2035 & 5X,'strain-rate exponent . . . . . . . . .=',1PG20.13/
2036 & 5X,'initial temperature . . . . . . . . . .=',1PG20.13/
2037 & 5X,'strain-rate cutting frequency . . . . .=',1PG20.13/
2038 & 5X,'strain-rate choice flag . . . . . . . .=',I10/
2039 & 5X,' vp=1 equivalent plastic strain rate'/
2040 & 5X,' vp=2 total strain rate(default)'/
2041 & 5X,' vp=3 deviatoric strain rate'/)
2042 1700 FORMAT(
2043 & 5X,'|| DATA for angle(deg) . . . . . . . .=',1PG20.13/
2044 & 5X,'shear reference point . . . . . . . . . '/
2045 & 5X,' first coordinate fsh1. . . . . .=',1PG20.13/
2046 & 5X,' second coordinate fsh2. . . . . .=',1PG20.13/
2047 & 5X,'shear/uniaxial hinge point. . . . . . . '/
2048 & 5X,' first coordinate hish1 . . . . .=',1PG20.13/
2049 & 5X,' second coordinate hish2 . . . . .=',1PG20.13/
2050 & 5X,'uniaxial reference point. . . . . . . . '/
2051 & 5X,' first coordinate fun1. . . . . .=',1PG20.13/
2052 & 5X,' second coordinate fun2. . . . . .=',1PG20.13/
2053 & 5X,' lankford coefficient . . . . . .=',1PG20.13/
2054 & 5X,'uniaxial/plane-strain hinge point . . . '/
2055 & 5X,' first coordinate hiun1 . . . . .=',1PG20.13/
2056 & 5X,' second coordinate hiun2 . . . . .=',1PG20.13/
2057 & 5X,'plane-strain reference point . . . . . '/
2058 & 5X,' first coordinate fps1. . . . . .=',1PG20.13/
2059 & 5X,' second coordinate fps2. . . . . .=',1PG20.13/
2060 & 5X,'plane-strain/biaxial hinge point. . . . '/
2061 & 5X,' first coordinate hips1 . . . . .=',1PG20.13/
2062 & 5X,' second coordinate hips2 . . . . .=',1PG20.13/
2063 & 5X,'biaxial reference point . . . . . . . . '/
2064 & 5X,' first coordinate fbi1. . . . . .=',1PG20.13/
2065 & 5X,' second coordinate fbi2. . . . . .=',1PG20.13/)
2066 1750 FORMAT(
2067 & 5X,'|| DATA for angle(deg) . . . . . . . .=',1PG20.13/
2068 & 5X,'shear hinge point . . . . . . . . . . . '/
2069 & 5X,' first coordinate hish1 . . . . .=',1PG20.13/
2070 & 5X,' second coordinate hish2 . . . . .=',1PG20.13/
2071 & 5X,' shear weight factor . . . . . . .=',1PG20.13/
2072 & 5X,'uniaxial reference point. . . . . . . . '/
2073 & 5X,' first coordinate fun1. . . . . .=',1PG20.13/
2074 & 5X,' second coordinate fun2. . . . . .=',1PG20.13/
2075 & 5X,' lankford coefficient . . . . . .=',1PG20.13/
2076 & 5X,'plane-strain hinge point . . . . . . . '/
2077 & 5X,' first coordinate hips1 . . . . .=',1PG20.13/
2078 & 5X,' second coordinate hips2 . . . . .=',1PG20.13/
2079 & 5X,' plane-strain weight factor. . . .=',1PG20.13/
2080 & 5X,'biaxial reference point . . . . . . . . '/
2081 & 5X,' first coordinate fbi1. . . . . .=',1PG20.13/
2082 & 5X,' second coordinate fbi2. . . . . .=',1PG20.13/)
2083c-----------------------------------------------------------------------
2084 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
the
end diagonal values have been computed in the(sparse) matrix id.SOL
dgeev
subroutine dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition dgeev.f:192
dgesv
subroutine dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)
DGESV computes the solution to system of linear equations A * X = B for GE matrices
Definition dgesv.f:122
hm_get_floatv
subroutine hm_get_floatv(name, rval, is_available, lsubmodel, unitab)
Definition hm_get_floatv.F:418
hm_get_floatv_dim
subroutine hm_get_floatv_dim(name, dim_fac, is_available, lsubmodel, unitab)
Definition hm_get_floatv_dim.F:161
hm_get_intv
subroutine hm_get_intv(name, ival, is_available, lsubmodel)
Definition hm_get_intv.F:528
hm_option_is_encrypted
subroutine hm_option_is_encrypted(is_encrypted)
Definition hm_option_is_encrypted.F:288
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
elbuftag_mod
Definition elbuftag_mod.F:230
message_mod
Definition message_mod.F:1257
names_and_titles_mod
Definition names_and_titles_mod.F:1002
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1008
submodel_mod
Definition submodel_mod.F:717
unitab_mod
Definition unitab_mod.F:517
second
real function second()
SECOND Using ETIME
Definition second_EXT_ETIME.f:35
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:895