1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! The program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine umat_single_crystal_creep(amat,iel,iint,kode,elconloc, 20 & emec,emec0,beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time, 21 & ttime,icmd,ielas, 22 & mi,nstate_,xstateini,xstate,stre,stiff,iorien,pgauss, 23 & orab,pnewdt) 24! 25! calculates stiffness and stresses for a user defined material 26! law 27! 28! icmd=3: calculates stress at mechanical strain 29! else: calculates stress at mechanical strain and the stiffness 30! matrix 31! 32! INPUT: 33! 34! amat material name 35! iel element number 36! iint integration point number 37! 38! kode material type (-100-#of constants entered 39! under *USER MATERIAL): can be used for materials 40! with varying number of constants 41! 42! elconloc(21) user defined constants defined by the keyword 43! card *USER MATERIAL (max. 21, actual # = 44! -kode-100), interpolated for the 45! actual temperature t1l 46! 47! emec(6) Lagrange mechanical strain tensor (component order: 48! 11,22,33,12,13,23) at the end of the increment 49! (thermal strains are subtracted) 50! emec0(6) Lagrange mechanical strain tensor at the start of the 51! increment (thermal strains are subtracted) 52! beta(6) residual stress tensor (the stress entered under 53! the keyword *INITIAL CONDITIONS,TYPE=STRESS) 54! 55! xokl(3,3) deformation gradient at the start of the increment 56! voj Jacobian at the start of the increment 57! xkl(3,3) deformation gradient at the end of the increment 58! vj Jacobian at the end of the increment 59! 60! ithermal 0: no thermal effects are taken into account 61! >0: thermal effects are taken into account (triggered 62! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE) 63! t1l temperature at the end of the increment 64! dtime time length of the increment 65! time step time at the end of the current increment 66! ttime total time at the start of the current step 67! 68! icmd not equal to 3: calculate stress and stiffness 69! 3: calculate only stress 70! ielas 0: no elastic iteration: irreversible effects 71! are allowed 72! 1: elastic iteration, i.e. no irreversible 73! deformation allowed 74! 75! mi(1) max. # of integration points per element in the 76! model 77! nstate_ max. # of state variables in the model 78! 79! xstateini(nstate_,mi(1),# of elements) 80! state variables at the start of the increment 81! xstate(nstate_,mi(1),# of elements) 82! state variables at the end of the increment 83! 84! stre(6) Piola-Kirchhoff stress of the second kind 85! at the start of the increment 86! 87! iorien number of the local coordinate axis system 88! in the integration point at stake (takes the value 89! 0 if no local system applies) 90! pgauss(3) global coordinates of the integration point 91! orab(7,*) description of all local coordinate systems. 92! If a local coordinate system applies the global 93! tensors can be obtained by premultiplying the local 94! tensors with skl(3,3). skl is determined by calling 95! the subroutine transformatrix: 96! call transformatrix(orab(1,iorien),pgauss,skl) 97! 98! 99! OUTPUT: 100! 101! xstate(nstate_,mi(1),# of elements) 102! updated state variables at the end of the increment 103! stre(6) Piola-Kirchhoff stress of the second kind at the 104! end of the increment 105! stiff(21): consistent tangent stiffness matrix in the material 106! frame of reference at the end of the increment. In 107! other words: the derivative of the PK2 stress with 108! respect to the Lagrangian strain tensor. The matrix 109! is supposed to be symmetric, only the upper half is 110! to be given in the same order as for a fully 111! anisotropic elastic material (*ELASTIC,TYPE=ANISO). 112! Notice that the matrix is an integral part of the 113! fourth order material tensor, i.e. the Voigt notation 114! is not used. 115! 116 implicit none 117! 118 integer convergence 119! 120 character*80 amat 121! 122 integer ithermal(*),icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien 123! 124 integer i,j,k,l,ipiv(18),info,neq,lda,ldb, 125 & nrhs,iplas,icounter 126! 127 real*8 ep0(6),ep(6),pnewdt, 128 & dg(18),ddg(18),xm(6,18),ck(18),cn(18), 129 & stri(6),htri(18),sg(18),r(6),xmc(6,18), 130 & t(6),gl(18,18),gr(18,18),ee(6),c1111,c1122,c1212,dd, 131 & skl(3,3),xmtran(3,3),ddsdde(6,6),xx(6,18) 132! 133 real*8 elconloc(*),stiff(21),emec(6),emec0(6),beta(6),stre(6), 134 & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*), 135 & elas(21),time,ttime 136! 137 real*8 xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*) 138! 139! crystallographic slip planes: 140! 141! 1. n=1,1,1 l=1,-1,0 142! 2. n=1,1,1 l=1,0,-1 143! 3. n=1,1,1 l=0,1,-1 144! 4. n=1,-1,1 l=0,1,1 145! 5. n=1,-1,1 l=1,0,-1 146! 6. n=1,-1,1 l=1,1,0 147! 7. n=1,-1,-1 l=0,1,-1 148! 8. n=1,-1,-1 l=1,0,1 149! 9. n=1,-1,-1 l=1,1,0 150! 10. n=1,1,-1 l=0,1,1 151! 11. n=1,1,-1 l=1,0,1 152! 12. n=1,1,-1 l=1,-1,0 153! 13. n=1,0,0 l=0,1,1 154! 14. n=1,0,0 l=0,1,-1 155! 15. n=0,1,0 l=1,0,1 156! 16. n=0,1,0 l=1,0,-1 157! 17. n=0,0,1 l=1,1,0 158! 18. n=0,0,1 l=1,-1,0 159! 160 xm=reshape(( 161 & /0.4082482904639d+00,-0.4082482904639d+00, 0.0000000000000d+00, 162 & 0.0000000000000d+00, 0.2041241452319d+00,-0.2041241452319d+00, 163 & 0.4082482904639d+00, 0.0000000000000d+00,-0.4082482904639d+00, 164 & 0.2041241452319d+00, 0.0000000000000d+00,-0.2041241452319d+00, 165 & 0.0000000000000d+00, 0.4082482904639d+00,-0.4082482904639d+00, 166 & 0.2041241452319d+00,-0.2041241452319d+00, 0.0000000000000d+00, 167 & 0.0000000000000d+00,-0.4082482904639d+00, 0.4082482904639d+00, 168 & 0.2041241452319d+00, 0.2041241452319d+00, 0.0000000000000d+00, 169 & 0.4082482904639d+00, 0.0000000000000d+00,-0.4082482904639d+00, 170 & -0.2041241452319d+00, 0.0000000000000d+00, 0.2041241452319d+00, 171 & 0.4082482904639d+00,-0.4082482904639d+00, 0.0000000000000d+00, 172 & 0.0000000000000d+00, 0.2041241452319d+00, 0.2041241452319d+00, 173 & 0.0000000000000d+00,-0.4082482904639d+00, 0.4082482904639d+00, 174 & 0.2041241452319d+00,-0.2041241452319d+00, 0.0000000000000d+00, 175 & 0.4082482904639d+00, 0.0000000000000d+00,-0.4082482904639d+00, 176 & -0.2041241452319d+00, 0.0000000000000d+00,-0.2041241452319d+00, 177 & 0.4082482904639d+00,-0.4082482904639d+00, 0.0000000000000d+00, 178 & 0.0000000000000d+00,-0.2041241452319d+00,-0.2041241452319d+00, 179 & 0.0000000000000d+00, 0.4082482904639d+00,-0.4082482904639d+00, 180 & 0.2041241452319d+00, 0.2041241452319d+00, 0.0000000000000d+00, 181 & 0.4082482904639d+00, 0.0000000000000d+00,-0.4082482904639d+00, 182 & 0.2041241452319d+00, 0.0000000000000d+00, 0.2041241452319d+00, 183 & 0.4082482904639d+00,-0.4082482904639d+00, 0.0000000000000d+00, 184 & 0.0000000000000d+00,-0.2041241452319d+00, 0.2041241452319d+00, 185 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 186 & 0.3535533905933d+00, 0.3535533905933d+00, 0.0000000000000d+00, 187 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 188 & 0.3535533905933d+00,-0.3535533905933d+00, 0.0000000000000d+00, 189 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 190 & 0.3535533905933d+00, 0.0000000000000d+00, 0.3535533905933d+00, 191 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 192 & 0.3535533905933d+00, 0.0000000000000d+00,-0.3535533905933d+00, 193 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 194 & 0.0000000000000d+00, 0.3535533905933d+00, 0.3535533905933d+00, 195 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 196 & 0.0000000000000d+00, 0.3535533905933d+00,-0.3535533905933d+00/ 197 & ),(/6,18/)) 198! 199 xx=reshape(( 200 & /0.4082482904639d+00,-0.4082482904639d+00, 0.0000000000000d+00, 201 & 0.0000000000000d+00, 0.2041241452319d+00,-0.2041241452319d+00, 202 & 0.4082482904639d+00, 0.0000000000000d+00,-0.4082482904639d+00, 203 & 0.2041241452319d+00, 0.0000000000000d+00,-0.2041241452319d+00, 204 & 0.0000000000000d+00, 0.4082482904639d+00,-0.4082482904639d+00, 205 & 0.2041241452319d+00,-0.2041241452319d+00, 0.0000000000000d+00, 206 & 0.0000000000000d+00,-0.4082482904639d+00, 0.4082482904639d+00, 207 & 0.2041241452319d+00, 0.2041241452319d+00, 0.0000000000000d+00, 208 & 0.4082482904639d+00, 0.0000000000000d+00,-0.4082482904639d+00, 209 & -0.2041241452319d+00, 0.0000000000000d+00, 0.2041241452319d+00, 210 & 0.4082482904639d+00,-0.4082482904639d+00, 0.0000000000000d+00, 211 & 0.0000000000000d+00, 0.2041241452319d+00, 0.2041241452319d+00, 212 & 0.0000000000000d+00,-0.4082482904639d+00, 0.4082482904639d+00, 213 & 0.2041241452319d+00,-0.2041241452319d+00, 0.0000000000000d+00, 214 & 0.4082482904639d+00, 0.0000000000000d+00,-0.4082482904639d+00, 215 & -0.2041241452319d+00, 0.0000000000000d+00,-0.2041241452319d+00, 216 & 0.4082482904639d+00,-0.4082482904639d+00, 0.0000000000000d+00, 217 & 0.0000000000000d+00,-0.2041241452319d+00,-0.2041241452319d+00, 218 & 0.0000000000000d+00, 0.4082482904639d+00,-0.4082482904639d+00, 219 & 0.2041241452319d+00, 0.2041241452319d+00, 0.0000000000000d+00, 220 & 0.4082482904639d+00, 0.0000000000000d+00,-0.4082482904639d+00, 221 & 0.2041241452319d+00, 0.0000000000000d+00, 0.2041241452319d+00, 222 & 0.4082482904639d+00,-0.4082482904639d+00, 0.0000000000000d+00, 223 & 0.0000000000000d+00,-0.2041241452319d+00, 0.2041241452319d+00, 224 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 225 & 0.3535533905933d+00, 0.3535533905933d+00, 0.0000000000000d+00, 226 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 227 & 0.3535533905933d+00,-0.3535533905933d+00, 0.0000000000000d+00, 228 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 229 & 0.3535533905933d+00, 0.0000000000000d+00, 0.3535533905933d+00, 230 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 231 & 0.3535533905933d+00, 0.0000000000000d+00,-0.3535533905933d+00, 232 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 233 & 0.0000000000000d+00, 0.3535533905933d+00, 0.3535533905933d+00, 234 & 0.0000000000000d+00, 0.0000000000000d+00, 0.0000000000000d+00, 235 & 0.0000000000000d+00, 0.3535533905933d+00,-0.3535533905933d+00/ 236 & ),(/6,18/)) 237! 238c pnewdt=-1.d0 239! 240! elastic constants 241! 242 c1111=elconloc(1) 243 c1122=elconloc(2) 244 c1212=elconloc(3) 245! 246 if(iorien.gt.0) then 247 call transformatrix(orab(1,iorien),pgauss,skl) 248 do k=1,18 249 do i=1,3 250 do j=i,3 251 xmtran(i,j)=skl(i,1)*skl(j,1)*xx(1,k)+ 252 & skl(i,2)*skl(j,2)*xx(2,k)+ 253 & skl(i,3)*skl(j,3)*xx(3,k)+ 254 & (skl(i,1)*skl(j,2)+ 255 & skl(i,2)*skl(j,1))*xx(4,k)+ 256 & (skl(i,1)*skl(j,3)+ 257 & skl(i,3)*skl(j,1))*xx(5,k)+ 258 & (skl(i,2)*skl(j,3)+ 259 & skl(i,3)*skl(j,2))*xx(6,k) 260 enddo 261 enddo 262 xm(1,k)=xmtran(1,1) 263 xm(2,k)=xmtran(2,2) 264 xm(3,k)=xmtran(3,3) 265 xm(4,k)=xmtran(1,2) 266 xm(5,k)=xmtran(1,3) 267 xm(6,k)=xmtran(2,3) 268 enddo 269! 270 elas( 1)= 271 & skl(1,1)*skl(1,1)*skl(1,1)*skl(1,1)*c1111+ 272 & skl(1,1)*skl(1,1)*skl(1,2)*skl(1,2)*c1122+ 273 & skl(1,1)*skl(1,1)*skl(1,3)*skl(1,3)*c1122+ 274 & skl(1,1)*skl(1,2)*skl(1,1)*skl(1,2)*c1212+ 275 & skl(1,1)*skl(1,2)*skl(1,2)*skl(1,1)*c1212+ 276 & skl(1,1)*skl(1,3)*skl(1,1)*skl(1,3)*c1212+ 277 & skl(1,1)*skl(1,3)*skl(1,3)*skl(1,1)*c1212+ 278 & skl(1,2)*skl(1,1)*skl(1,1)*skl(1,2)*c1212+ 279 & skl(1,2)*skl(1,1)*skl(1,2)*skl(1,1)*c1212+ 280 & skl(1,2)*skl(1,2)*skl(1,1)*skl(1,1)*c1122+ 281 & skl(1,2)*skl(1,2)*skl(1,2)*skl(1,2)*c1111+ 282 & skl(1,2)*skl(1,2)*skl(1,3)*skl(1,3)*c1122+ 283 & skl(1,2)*skl(1,3)*skl(1,2)*skl(1,3)*c1212+ 284 & skl(1,2)*skl(1,3)*skl(1,3)*skl(1,2)*c1212+ 285 & skl(1,3)*skl(1,1)*skl(1,1)*skl(1,3)*c1212+ 286 & skl(1,3)*skl(1,1)*skl(1,3)*skl(1,1)*c1212+ 287 & skl(1,3)*skl(1,2)*skl(1,2)*skl(1,3)*c1212+ 288 & skl(1,3)*skl(1,2)*skl(1,3)*skl(1,2)*c1212+ 289 & skl(1,3)*skl(1,3)*skl(1,1)*skl(1,1)*c1122+ 290 & skl(1,3)*skl(1,3)*skl(1,2)*skl(1,2)*c1122+ 291 & skl(1,3)*skl(1,3)*skl(1,3)*skl(1,3)*c1111 292 elas( 2)= 293 & skl(1,1)*skl(1,1)*skl(2,1)*skl(2,1)*c1111+ 294 & skl(1,1)*skl(1,1)*skl(2,2)*skl(2,2)*c1122+ 295 & skl(1,1)*skl(1,1)*skl(2,3)*skl(2,3)*c1122+ 296 & skl(1,1)*skl(1,2)*skl(2,1)*skl(2,2)*c1212+ 297 & skl(1,1)*skl(1,2)*skl(2,2)*skl(2,1)*c1212+ 298 & skl(1,1)*skl(1,3)*skl(2,1)*skl(2,3)*c1212+ 299 & skl(1,1)*skl(1,3)*skl(2,3)*skl(2,1)*c1212+ 300 & skl(1,2)*skl(1,1)*skl(2,1)*skl(2,2)*c1212+ 301 & skl(1,2)*skl(1,1)*skl(2,2)*skl(2,1)*c1212+ 302 & skl(1,2)*skl(1,2)*skl(2,1)*skl(2,1)*c1122+ 303 & skl(1,2)*skl(1,2)*skl(2,2)*skl(2,2)*c1111+ 304 & skl(1,2)*skl(1,2)*skl(2,3)*skl(2,3)*c1122+ 305 & skl(1,2)*skl(1,3)*skl(2,2)*skl(2,3)*c1212+ 306 & skl(1,2)*skl(1,3)*skl(2,3)*skl(2,2)*c1212+ 307 & skl(1,3)*skl(1,1)*skl(2,1)*skl(2,3)*c1212+ 308 & skl(1,3)*skl(1,1)*skl(2,3)*skl(2,1)*c1212+ 309 & skl(1,3)*skl(1,2)*skl(2,2)*skl(2,3)*c1212+ 310 & skl(1,3)*skl(1,2)*skl(2,3)*skl(2,2)*c1212+ 311 & skl(1,3)*skl(1,3)*skl(2,1)*skl(2,1)*c1122+ 312 & skl(1,3)*skl(1,3)*skl(2,2)*skl(2,2)*c1122+ 313 & skl(1,3)*skl(1,3)*skl(2,3)*skl(2,3)*c1111 314 elas( 3)= 315 & skl(2,1)*skl(2,1)*skl(2,1)*skl(2,1)*c1111+ 316 & skl(2,1)*skl(2,1)*skl(2,2)*skl(2,2)*c1122+ 317 & skl(2,1)*skl(2,1)*skl(2,3)*skl(2,3)*c1122+ 318 & skl(2,1)*skl(2,2)*skl(2,1)*skl(2,2)*c1212+ 319 & skl(2,1)*skl(2,2)*skl(2,2)*skl(2,1)*c1212+ 320 & skl(2,1)*skl(2,3)*skl(2,1)*skl(2,3)*c1212+ 321 & skl(2,1)*skl(2,3)*skl(2,3)*skl(2,1)*c1212+ 322 & skl(2,2)*skl(2,1)*skl(2,1)*skl(2,2)*c1212+ 323 & skl(2,2)*skl(2,1)*skl(2,2)*skl(2,1)*c1212+ 324 & skl(2,2)*skl(2,2)*skl(2,1)*skl(2,1)*c1122+ 325 & skl(2,2)*skl(2,2)*skl(2,2)*skl(2,2)*c1111+ 326 & skl(2,2)*skl(2,2)*skl(2,3)*skl(2,3)*c1122+ 327 & skl(2,2)*skl(2,3)*skl(2,2)*skl(2,3)*c1212+ 328 & skl(2,2)*skl(2,3)*skl(2,3)*skl(2,2)*c1212+ 329 & skl(2,3)*skl(2,1)*skl(2,1)*skl(2,3)*c1212+ 330 & skl(2,3)*skl(2,1)*skl(2,3)*skl(2,1)*c1212+ 331 & skl(2,3)*skl(2,2)*skl(2,2)*skl(2,3)*c1212+ 332 & skl(2,3)*skl(2,2)*skl(2,3)*skl(2,2)*c1212+ 333 & skl(2,3)*skl(2,3)*skl(2,1)*skl(2,1)*c1122+ 334 & skl(2,3)*skl(2,3)*skl(2,2)*skl(2,2)*c1122+ 335 & skl(2,3)*skl(2,3)*skl(2,3)*skl(2,3)*c1111 336 elas( 4)= 337 & skl(1,1)*skl(1,1)*skl(3,1)*skl(3,1)*c1111+ 338 & skl(1,1)*skl(1,1)*skl(3,2)*skl(3,2)*c1122+ 339 & skl(1,1)*skl(1,1)*skl(3,3)*skl(3,3)*c1122+ 340 & skl(1,1)*skl(1,2)*skl(3,1)*skl(3,2)*c1212+ 341 & skl(1,1)*skl(1,2)*skl(3,2)*skl(3,1)*c1212+ 342 & skl(1,1)*skl(1,3)*skl(3,1)*skl(3,3)*c1212+ 343 & skl(1,1)*skl(1,3)*skl(3,3)*skl(3,1)*c1212+ 344 & skl(1,2)*skl(1,1)*skl(3,1)*skl(3,2)*c1212+ 345 & skl(1,2)*skl(1,1)*skl(3,2)*skl(3,1)*c1212+ 346 & skl(1,2)*skl(1,2)*skl(3,1)*skl(3,1)*c1122+ 347 & skl(1,2)*skl(1,2)*skl(3,2)*skl(3,2)*c1111+ 348 & skl(1,2)*skl(1,2)*skl(3,3)*skl(3,3)*c1122+ 349 & skl(1,2)*skl(1,3)*skl(3,2)*skl(3,3)*c1212+ 350 & skl(1,2)*skl(1,3)*skl(3,3)*skl(3,2)*c1212+ 351 & skl(1,3)*skl(1,1)*skl(3,1)*skl(3,3)*c1212+ 352 & skl(1,3)*skl(1,1)*skl(3,3)*skl(3,1)*c1212+ 353 & skl(1,3)*skl(1,2)*skl(3,2)*skl(3,3)*c1212+ 354 & skl(1,3)*skl(1,2)*skl(3,3)*skl(3,2)*c1212+ 355 & skl(1,3)*skl(1,3)*skl(3,1)*skl(3,1)*c1122+ 356 & skl(1,3)*skl(1,3)*skl(3,2)*skl(3,2)*c1122+ 357 & skl(1,3)*skl(1,3)*skl(3,3)*skl(3,3)*c1111 358 elas( 5)= 359 & skl(2,1)*skl(2,1)*skl(3,1)*skl(3,1)*c1111+ 360 & skl(2,1)*skl(2,1)*skl(3,2)*skl(3,2)*c1122+ 361 & skl(2,1)*skl(2,1)*skl(3,3)*skl(3,3)*c1122+ 362 & skl(2,1)*skl(2,2)*skl(3,1)*skl(3,2)*c1212+ 363 & skl(2,1)*skl(2,2)*skl(3,2)*skl(3,1)*c1212+ 364 & skl(2,1)*skl(2,3)*skl(3,1)*skl(3,3)*c1212+ 365 & skl(2,1)*skl(2,3)*skl(3,3)*skl(3,1)*c1212+ 366 & skl(2,2)*skl(2,1)*skl(3,1)*skl(3,2)*c1212+ 367 & skl(2,2)*skl(2,1)*skl(3,2)*skl(3,1)*c1212+ 368 & skl(2,2)*skl(2,2)*skl(3,1)*skl(3,1)*c1122+ 369 & skl(2,2)*skl(2,2)*skl(3,2)*skl(3,2)*c1111+ 370 & skl(2,2)*skl(2,2)*skl(3,3)*skl(3,3)*c1122+ 371 & skl(2,2)*skl(2,3)*skl(3,2)*skl(3,3)*c1212+ 372 & skl(2,2)*skl(2,3)*skl(3,3)*skl(3,2)*c1212+ 373 & skl(2,3)*skl(2,1)*skl(3,1)*skl(3,3)*c1212+ 374 & skl(2,3)*skl(2,1)*skl(3,3)*skl(3,1)*c1212+ 375 & skl(2,3)*skl(2,2)*skl(3,2)*skl(3,3)*c1212+ 376 & skl(2,3)*skl(2,2)*skl(3,3)*skl(3,2)*c1212+ 377 & skl(2,3)*skl(2,3)*skl(3,1)*skl(3,1)*c1122+ 378 & skl(2,3)*skl(2,3)*skl(3,2)*skl(3,2)*c1122+ 379 & skl(2,3)*skl(2,3)*skl(3,3)*skl(3,3)*c1111 380 elas( 6)= 381 & skl(3,1)*skl(3,1)*skl(3,1)*skl(3,1)*c1111+ 382 & skl(3,1)*skl(3,1)*skl(3,2)*skl(3,2)*c1122+ 383 & skl(3,1)*skl(3,1)*skl(3,3)*skl(3,3)*c1122+ 384 & skl(3,1)*skl(3,2)*skl(3,1)*skl(3,2)*c1212+ 385 & skl(3,1)*skl(3,2)*skl(3,2)*skl(3,1)*c1212+ 386 & skl(3,1)*skl(3,3)*skl(3,1)*skl(3,3)*c1212+ 387 & skl(3,1)*skl(3,3)*skl(3,3)*skl(3,1)*c1212+ 388 & skl(3,2)*skl(3,1)*skl(3,1)*skl(3,2)*c1212+ 389 & skl(3,2)*skl(3,1)*skl(3,2)*skl(3,1)*c1212+ 390 & skl(3,2)*skl(3,2)*skl(3,1)*skl(3,1)*c1122+ 391 & skl(3,2)*skl(3,2)*skl(3,2)*skl(3,2)*c1111+ 392 & skl(3,2)*skl(3,2)*skl(3,3)*skl(3,3)*c1122+ 393 & skl(3,2)*skl(3,3)*skl(3,2)*skl(3,3)*c1212+ 394 & skl(3,2)*skl(3,3)*skl(3,3)*skl(3,2)*c1212+ 395 & skl(3,3)*skl(3,1)*skl(3,1)*skl(3,3)*c1212+ 396 & skl(3,3)*skl(3,1)*skl(3,3)*skl(3,1)*c1212+ 397 & skl(3,3)*skl(3,2)*skl(3,2)*skl(3,3)*c1212+ 398 & skl(3,3)*skl(3,2)*skl(3,3)*skl(3,2)*c1212+ 399 & skl(3,3)*skl(3,3)*skl(3,1)*skl(3,1)*c1122+ 400 & skl(3,3)*skl(3,3)*skl(3,2)*skl(3,2)*c1122+ 401 & skl(3,3)*skl(3,3)*skl(3,3)*skl(3,3)*c1111 402 elas( 7)= 403 & skl(1,1)*skl(1,1)*skl(1,1)*skl(2,1)*c1111+ 404 & skl(1,1)*skl(1,1)*skl(1,2)*skl(2,2)*c1122+ 405 & skl(1,1)*skl(1,1)*skl(1,3)*skl(2,3)*c1122+ 406 & skl(1,1)*skl(1,2)*skl(1,1)*skl(2,2)*c1212+ 407 & skl(1,1)*skl(1,2)*skl(1,2)*skl(2,1)*c1212+ 408 & skl(1,1)*skl(1,3)*skl(1,1)*skl(2,3)*c1212+ 409 & skl(1,1)*skl(1,3)*skl(1,3)*skl(2,1)*c1212+ 410 & skl(1,2)*skl(1,1)*skl(1,1)*skl(2,2)*c1212+ 411 & skl(1,2)*skl(1,1)*skl(1,2)*skl(2,1)*c1212+ 412 & skl(1,2)*skl(1,2)*skl(1,1)*skl(2,1)*c1122+ 413 & skl(1,2)*skl(1,2)*skl(1,2)*skl(2,2)*c1111+ 414 & skl(1,2)*skl(1,2)*skl(1,3)*skl(2,3)*c1122+ 415 & skl(1,2)*skl(1,3)*skl(1,2)*skl(2,3)*c1212+ 416 & skl(1,2)*skl(1,3)*skl(1,3)*skl(2,2)*c1212+ 417 & skl(1,3)*skl(1,1)*skl(1,1)*skl(2,3)*c1212+ 418 & skl(1,3)*skl(1,1)*skl(1,3)*skl(2,1)*c1212+ 419 & skl(1,3)*skl(1,2)*skl(1,2)*skl(2,3)*c1212+ 420 & skl(1,3)*skl(1,2)*skl(1,3)*skl(2,2)*c1212+ 421 & skl(1,3)*skl(1,3)*skl(1,1)*skl(2,1)*c1122+ 422 & skl(1,3)*skl(1,3)*skl(1,2)*skl(2,2)*c1122+ 423 & skl(1,3)*skl(1,3)*skl(1,3)*skl(2,3)*c1111 424 elas( 8)= 425 & skl(2,1)*skl(2,1)*skl(1,1)*skl(2,1)*c1111+ 426 & skl(2,1)*skl(2,1)*skl(1,2)*skl(2,2)*c1122+ 427 & skl(2,1)*skl(2,1)*skl(1,3)*skl(2,3)*c1122+ 428 & skl(2,1)*skl(2,2)*skl(1,1)*skl(2,2)*c1212+ 429 & skl(2,1)*skl(2,2)*skl(1,2)*skl(2,1)*c1212+ 430 & skl(2,1)*skl(2,3)*skl(1,1)*skl(2,3)*c1212+ 431 & skl(2,1)*skl(2,3)*skl(1,3)*skl(2,1)*c1212+ 432 & skl(2,2)*skl(2,1)*skl(1,1)*skl(2,2)*c1212+ 433 & skl(2,2)*skl(2,1)*skl(1,2)*skl(2,1)*c1212+ 434 & skl(2,2)*skl(2,2)*skl(1,1)*skl(2,1)*c1122+ 435 & skl(2,2)*skl(2,2)*skl(1,2)*skl(2,2)*c1111+ 436 & skl(2,2)*skl(2,2)*skl(1,3)*skl(2,3)*c1122+ 437 & skl(2,2)*skl(2,3)*skl(1,2)*skl(2,3)*c1212+ 438 & skl(2,2)*skl(2,3)*skl(1,3)*skl(2,2)*c1212+ 439 & skl(2,3)*skl(2,1)*skl(1,1)*skl(2,3)*c1212+ 440 & skl(2,3)*skl(2,1)*skl(1,3)*skl(2,1)*c1212+ 441 & skl(2,3)*skl(2,2)*skl(1,2)*skl(2,3)*c1212+ 442 & skl(2,3)*skl(2,2)*skl(1,3)*skl(2,2)*c1212+ 443 & skl(2,3)*skl(2,3)*skl(1,1)*skl(2,1)*c1122+ 444 & skl(2,3)*skl(2,3)*skl(1,2)*skl(2,2)*c1122+ 445 & skl(2,3)*skl(2,3)*skl(1,3)*skl(2,3)*c1111 446 elas( 9)= 447 & skl(3,1)*skl(3,1)*skl(1,1)*skl(2,1)*c1111+ 448 & skl(3,1)*skl(3,1)*skl(1,2)*skl(2,2)*c1122+ 449 & skl(3,1)*skl(3,1)*skl(1,3)*skl(2,3)*c1122+ 450 & skl(3,1)*skl(3,2)*skl(1,1)*skl(2,2)*c1212+ 451 & skl(3,1)*skl(3,2)*skl(1,2)*skl(2,1)*c1212+ 452 & skl(3,1)*skl(3,3)*skl(1,1)*skl(2,3)*c1212+ 453 & skl(3,1)*skl(3,3)*skl(1,3)*skl(2,1)*c1212+ 454 & skl(3,2)*skl(3,1)*skl(1,1)*skl(2,2)*c1212+ 455 & skl(3,2)*skl(3,1)*skl(1,2)*skl(2,1)*c1212+ 456 & skl(3,2)*skl(3,2)*skl(1,1)*skl(2,1)*c1122+ 457 & skl(3,2)*skl(3,2)*skl(1,2)*skl(2,2)*c1111+ 458 & skl(3,2)*skl(3,2)*skl(1,3)*skl(2,3)*c1122+ 459 & skl(3,2)*skl(3,3)*skl(1,2)*skl(2,3)*c1212+ 460 & skl(3,2)*skl(3,3)*skl(1,3)*skl(2,2)*c1212+ 461 & skl(3,3)*skl(3,1)*skl(1,1)*skl(2,3)*c1212+ 462 & skl(3,3)*skl(3,1)*skl(1,3)*skl(2,1)*c1212+ 463 & skl(3,3)*skl(3,2)*skl(1,2)*skl(2,3)*c1212+ 464 & skl(3,3)*skl(3,2)*skl(1,3)*skl(2,2)*c1212+ 465 & skl(3,3)*skl(3,3)*skl(1,1)*skl(2,1)*c1122+ 466 & skl(3,3)*skl(3,3)*skl(1,2)*skl(2,2)*c1122+ 467 & skl(3,3)*skl(3,3)*skl(1,3)*skl(2,3)*c1111 468 elas(10)= 469 & skl(1,1)*skl(2,1)*skl(1,1)*skl(2,1)*c1111+ 470 & skl(1,1)*skl(2,1)*skl(1,2)*skl(2,2)*c1122+ 471 & skl(1,1)*skl(2,1)*skl(1,3)*skl(2,3)*c1122+ 472 & skl(1,1)*skl(2,2)*skl(1,1)*skl(2,2)*c1212+ 473 & skl(1,1)*skl(2,2)*skl(1,2)*skl(2,1)*c1212+ 474 & skl(1,1)*skl(2,3)*skl(1,1)*skl(2,3)*c1212+ 475 & skl(1,1)*skl(2,3)*skl(1,3)*skl(2,1)*c1212+ 476 & skl(1,2)*skl(2,1)*skl(1,1)*skl(2,2)*c1212+ 477 & skl(1,2)*skl(2,1)*skl(1,2)*skl(2,1)*c1212+ 478 & skl(1,2)*skl(2,2)*skl(1,1)*skl(2,1)*c1122+ 479 & skl(1,2)*skl(2,2)*skl(1,2)*skl(2,2)*c1111+ 480 & skl(1,2)*skl(2,2)*skl(1,3)*skl(2,3)*c1122+ 481 & skl(1,2)*skl(2,3)*skl(1,2)*skl(2,3)*c1212+ 482 & skl(1,2)*skl(2,3)*skl(1,3)*skl(2,2)*c1212+ 483 & skl(1,3)*skl(2,1)*skl(1,1)*skl(2,3)*c1212+ 484 & skl(1,3)*skl(2,1)*skl(1,3)*skl(2,1)*c1212+ 485 & skl(1,3)*skl(2,2)*skl(1,2)*skl(2,3)*c1212+ 486 & skl(1,3)*skl(2,2)*skl(1,3)*skl(2,2)*c1212+ 487 & skl(1,3)*skl(2,3)*skl(1,1)*skl(2,1)*c1122+ 488 & skl(1,3)*skl(2,3)*skl(1,2)*skl(2,2)*c1122+ 489 & skl(1,3)*skl(2,3)*skl(1,3)*skl(2,3)*c1111 490 elas(11)= 491 & skl(1,1)*skl(1,1)*skl(1,1)*skl(3,1)*c1111+ 492 & skl(1,1)*skl(1,1)*skl(1,2)*skl(3,2)*c1122+ 493 & skl(1,1)*skl(1,1)*skl(1,3)*skl(3,3)*c1122+ 494 & skl(1,1)*skl(1,2)*skl(1,1)*skl(3,2)*c1212+ 495 & skl(1,1)*skl(1,2)*skl(1,2)*skl(3,1)*c1212+ 496 & skl(1,1)*skl(1,3)*skl(1,1)*skl(3,3)*c1212+ 497 & skl(1,1)*skl(1,3)*skl(1,3)*skl(3,1)*c1212+ 498 & skl(1,2)*skl(1,1)*skl(1,1)*skl(3,2)*c1212+ 499 & skl(1,2)*skl(1,1)*skl(1,2)*skl(3,1)*c1212+ 500 & skl(1,2)*skl(1,2)*skl(1,1)*skl(3,1)*c1122+ 501 & skl(1,2)*skl(1,2)*skl(1,2)*skl(3,2)*c1111+ 502 & skl(1,2)*skl(1,2)*skl(1,3)*skl(3,3)*c1122+ 503 & skl(1,2)*skl(1,3)*skl(1,2)*skl(3,3)*c1212+ 504 & skl(1,2)*skl(1,3)*skl(1,3)*skl(3,2)*c1212+ 505 & skl(1,3)*skl(1,1)*skl(1,1)*skl(3,3)*c1212+ 506 & skl(1,3)*skl(1,1)*skl(1,3)*skl(3,1)*c1212+ 507 & skl(1,3)*skl(1,2)*skl(1,2)*skl(3,3)*c1212+ 508 & skl(1,3)*skl(1,2)*skl(1,3)*skl(3,2)*c1212+ 509 & skl(1,3)*skl(1,3)*skl(1,1)*skl(3,1)*c1122+ 510 & skl(1,3)*skl(1,3)*skl(1,2)*skl(3,2)*c1122+ 511 & skl(1,3)*skl(1,3)*skl(1,3)*skl(3,3)*c1111 512 elas(12)= 513 & skl(2,1)*skl(2,1)*skl(1,1)*skl(3,1)*c1111+ 514 & skl(2,1)*skl(2,1)*skl(1,2)*skl(3,2)*c1122+ 515 & skl(2,1)*skl(2,1)*skl(1,3)*skl(3,3)*c1122+ 516 & skl(2,1)*skl(2,2)*skl(1,1)*skl(3,2)*c1212+ 517 & skl(2,1)*skl(2,2)*skl(1,2)*skl(3,1)*c1212+ 518 & skl(2,1)*skl(2,3)*skl(1,1)*skl(3,3)*c1212+ 519 & skl(2,1)*skl(2,3)*skl(1,3)*skl(3,1)*c1212+ 520 & skl(2,2)*skl(2,1)*skl(1,1)*skl(3,2)*c1212+ 521 & skl(2,2)*skl(2,1)*skl(1,2)*skl(3,1)*c1212+ 522 & skl(2,2)*skl(2,2)*skl(1,1)*skl(3,1)*c1122+ 523 & skl(2,2)*skl(2,2)*skl(1,2)*skl(3,2)*c1111+ 524 & skl(2,2)*skl(2,2)*skl(1,3)*skl(3,3)*c1122+ 525 & skl(2,2)*skl(2,3)*skl(1,2)*skl(3,3)*c1212+ 526 & skl(2,2)*skl(2,3)*skl(1,3)*skl(3,2)*c1212+ 527 & skl(2,3)*skl(2,1)*skl(1,1)*skl(3,3)*c1212+ 528 & skl(2,3)*skl(2,1)*skl(1,3)*skl(3,1)*c1212+ 529 & skl(2,3)*skl(2,2)*skl(1,2)*skl(3,3)*c1212+ 530 & skl(2,3)*skl(2,2)*skl(1,3)*skl(3,2)*c1212+ 531 & skl(2,3)*skl(2,3)*skl(1,1)*skl(3,1)*c1122+ 532 & skl(2,3)*skl(2,3)*skl(1,2)*skl(3,2)*c1122+ 533 & skl(2,3)*skl(2,3)*skl(1,3)*skl(3,3)*c1111 534 elas(13)= 535 & skl(3,1)*skl(3,1)*skl(1,1)*skl(3,1)*c1111+ 536 & skl(3,1)*skl(3,1)*skl(1,2)*skl(3,2)*c1122+ 537 & skl(3,1)*skl(3,1)*skl(1,3)*skl(3,3)*c1122+ 538 & skl(3,1)*skl(3,2)*skl(1,1)*skl(3,2)*c1212+ 539 & skl(3,1)*skl(3,2)*skl(1,2)*skl(3,1)*c1212+ 540 & skl(3,1)*skl(3,3)*skl(1,1)*skl(3,3)*c1212+ 541 & skl(3,1)*skl(3,3)*skl(1,3)*skl(3,1)*c1212+ 542 & skl(3,2)*skl(3,1)*skl(1,1)*skl(3,2)*c1212+ 543 & skl(3,2)*skl(3,1)*skl(1,2)*skl(3,1)*c1212+ 544 & skl(3,2)*skl(3,2)*skl(1,1)*skl(3,1)*c1122+ 545 & skl(3,2)*skl(3,2)*skl(1,2)*skl(3,2)*c1111+ 546 & skl(3,2)*skl(3,2)*skl(1,3)*skl(3,3)*c1122+ 547 & skl(3,2)*skl(3,3)*skl(1,2)*skl(3,3)*c1212+ 548 & skl(3,2)*skl(3,3)*skl(1,3)*skl(3,2)*c1212+ 549 & skl(3,3)*skl(3,1)*skl(1,1)*skl(3,3)*c1212+ 550 & skl(3,3)*skl(3,1)*skl(1,3)*skl(3,1)*c1212+ 551 & skl(3,3)*skl(3,2)*skl(1,2)*skl(3,3)*c1212+ 552 & skl(3,3)*skl(3,2)*skl(1,3)*skl(3,2)*c1212+ 553 & skl(3,3)*skl(3,3)*skl(1,1)*skl(3,1)*c1122+ 554 & skl(3,3)*skl(3,3)*skl(1,2)*skl(3,2)*c1122+ 555 & skl(3,3)*skl(3,3)*skl(1,3)*skl(3,3)*c1111 556 elas(14)= 557 & skl(1,1)*skl(2,1)*skl(1,1)*skl(3,1)*c1111+ 558 & skl(1,1)*skl(2,1)*skl(1,2)*skl(3,2)*c1122+ 559 & skl(1,1)*skl(2,1)*skl(1,3)*skl(3,3)*c1122+ 560 & skl(1,1)*skl(2,2)*skl(1,1)*skl(3,2)*c1212+ 561 & skl(1,1)*skl(2,2)*skl(1,2)*skl(3,1)*c1212+ 562 & skl(1,1)*skl(2,3)*skl(1,1)*skl(3,3)*c1212+ 563 & skl(1,1)*skl(2,3)*skl(1,3)*skl(3,1)*c1212+ 564 & skl(1,2)*skl(2,1)*skl(1,1)*skl(3,2)*c1212+ 565 & skl(1,2)*skl(2,1)*skl(1,2)*skl(3,1)*c1212+ 566 & skl(1,2)*skl(2,2)*skl(1,1)*skl(3,1)*c1122+ 567 & skl(1,2)*skl(2,2)*skl(1,2)*skl(3,2)*c1111+ 568 & skl(1,2)*skl(2,2)*skl(1,3)*skl(3,3)*c1122+ 569 & skl(1,2)*skl(2,3)*skl(1,2)*skl(3,3)*c1212+ 570 & skl(1,2)*skl(2,3)*skl(1,3)*skl(3,2)*c1212+ 571 & skl(1,3)*skl(2,1)*skl(1,1)*skl(3,3)*c1212+ 572 & skl(1,3)*skl(2,1)*skl(1,3)*skl(3,1)*c1212+ 573 & skl(1,3)*skl(2,2)*skl(1,2)*skl(3,3)*c1212+ 574 & skl(1,3)*skl(2,2)*skl(1,3)*skl(3,2)*c1212+ 575 & skl(1,3)*skl(2,3)*skl(1,1)*skl(3,1)*c1122+ 576 & skl(1,3)*skl(2,3)*skl(1,2)*skl(3,2)*c1122+ 577 & skl(1,3)*skl(2,3)*skl(1,3)*skl(3,3)*c1111 578 elas(15)= 579 & skl(1,1)*skl(3,1)*skl(1,1)*skl(3,1)*c1111+ 580 & skl(1,1)*skl(3,1)*skl(1,2)*skl(3,2)*c1122+ 581 & skl(1,1)*skl(3,1)*skl(1,3)*skl(3,3)*c1122+ 582 & skl(1,1)*skl(3,2)*skl(1,1)*skl(3,2)*c1212+ 583 & skl(1,1)*skl(3,2)*skl(1,2)*skl(3,1)*c1212+ 584 & skl(1,1)*skl(3,3)*skl(1,1)*skl(3,3)*c1212+ 585 & skl(1,1)*skl(3,3)*skl(1,3)*skl(3,1)*c1212+ 586 & skl(1,2)*skl(3,1)*skl(1,1)*skl(3,2)*c1212+ 587 & skl(1,2)*skl(3,1)*skl(1,2)*skl(3,1)*c1212+ 588 & skl(1,2)*skl(3,2)*skl(1,1)*skl(3,1)*c1122+ 589 & skl(1,2)*skl(3,2)*skl(1,2)*skl(3,2)*c1111+ 590 & skl(1,2)*skl(3,2)*skl(1,3)*skl(3,3)*c1122+ 591 & skl(1,2)*skl(3,3)*skl(1,2)*skl(3,3)*c1212+ 592 & skl(1,2)*skl(3,3)*skl(1,3)*skl(3,2)*c1212+ 593 & skl(1,3)*skl(3,1)*skl(1,1)*skl(3,3)*c1212+ 594 & skl(1,3)*skl(3,1)*skl(1,3)*skl(3,1)*c1212+ 595 & skl(1,3)*skl(3,2)*skl(1,2)*skl(3,3)*c1212+ 596 & skl(1,3)*skl(3,2)*skl(1,3)*skl(3,2)*c1212+ 597 & skl(1,3)*skl(3,3)*skl(1,1)*skl(3,1)*c1122+ 598 & skl(1,3)*skl(3,3)*skl(1,2)*skl(3,2)*c1122+ 599 & skl(1,3)*skl(3,3)*skl(1,3)*skl(3,3)*c1111 600 elas(16)= 601 & skl(1,1)*skl(1,1)*skl(2,1)*skl(3,1)*c1111+ 602 & skl(1,1)*skl(1,1)*skl(2,2)*skl(3,2)*c1122+ 603 & skl(1,1)*skl(1,1)*skl(2,3)*skl(3,3)*c1122+ 604 & skl(1,1)*skl(1,2)*skl(2,1)*skl(3,2)*c1212+ 605 & skl(1,1)*skl(1,2)*skl(2,2)*skl(3,1)*c1212+ 606 & skl(1,1)*skl(1,3)*skl(2,1)*skl(3,3)*c1212+ 607 & skl(1,1)*skl(1,3)*skl(2,3)*skl(3,1)*c1212+ 608 & skl(1,2)*skl(1,1)*skl(2,1)*skl(3,2)*c1212+ 609 & skl(1,2)*skl(1,1)*skl(2,2)*skl(3,1)*c1212+ 610 & skl(1,2)*skl(1,2)*skl(2,1)*skl(3,1)*c1122+ 611 & skl(1,2)*skl(1,2)*skl(2,2)*skl(3,2)*c1111+ 612 & skl(1,2)*skl(1,2)*skl(2,3)*skl(3,3)*c1122+ 613 & skl(1,2)*skl(1,3)*skl(2,2)*skl(3,3)*c1212+ 614 & skl(1,2)*skl(1,3)*skl(2,3)*skl(3,2)*c1212+ 615 & skl(1,3)*skl(1,1)*skl(2,1)*skl(3,3)*c1212+ 616 & skl(1,3)*skl(1,1)*skl(2,3)*skl(3,1)*c1212+ 617 & skl(1,3)*skl(1,2)*skl(2,2)*skl(3,3)*c1212+ 618 & skl(1,3)*skl(1,2)*skl(2,3)*skl(3,2)*c1212+ 619 & skl(1,3)*skl(1,3)*skl(2,1)*skl(3,1)*c1122+ 620 & skl(1,3)*skl(1,3)*skl(2,2)*skl(3,2)*c1122+ 621 & skl(1,3)*skl(1,3)*skl(2,3)*skl(3,3)*c1111 622 elas(17)= 623 & skl(2,1)*skl(2,1)*skl(2,1)*skl(3,1)*c1111+ 624 & skl(2,1)*skl(2,1)*skl(2,2)*skl(3,2)*c1122+ 625 & skl(2,1)*skl(2,1)*skl(2,3)*skl(3,3)*c1122+ 626 & skl(2,1)*skl(2,2)*skl(2,1)*skl(3,2)*c1212+ 627 & skl(2,1)*skl(2,2)*skl(2,2)*skl(3,1)*c1212+ 628 & skl(2,1)*skl(2,3)*skl(2,1)*skl(3,3)*c1212+ 629 & skl(2,1)*skl(2,3)*skl(2,3)*skl(3,1)*c1212+ 630 & skl(2,2)*skl(2,1)*skl(2,1)*skl(3,2)*c1212+ 631 & skl(2,2)*skl(2,1)*skl(2,2)*skl(3,1)*c1212+ 632 & skl(2,2)*skl(2,2)*skl(2,1)*skl(3,1)*c1122+ 633 & skl(2,2)*skl(2,2)*skl(2,2)*skl(3,2)*c1111+ 634 & skl(2,2)*skl(2,2)*skl(2,3)*skl(3,3)*c1122+ 635 & skl(2,2)*skl(2,3)*skl(2,2)*skl(3,3)*c1212+ 636 & skl(2,2)*skl(2,3)*skl(2,3)*skl(3,2)*c1212+ 637 & skl(2,3)*skl(2,1)*skl(2,1)*skl(3,3)*c1212+ 638 & skl(2,3)*skl(2,1)*skl(2,3)*skl(3,1)*c1212+ 639 & skl(2,3)*skl(2,2)*skl(2,2)*skl(3,3)*c1212+ 640 & skl(2,3)*skl(2,2)*skl(2,3)*skl(3,2)*c1212+ 641 & skl(2,3)*skl(2,3)*skl(2,1)*skl(3,1)*c1122+ 642 & skl(2,3)*skl(2,3)*skl(2,2)*skl(3,2)*c1122+ 643 & skl(2,3)*skl(2,3)*skl(2,3)*skl(3,3)*c1111 644 elas(18)= 645 & skl(3,1)*skl(3,1)*skl(2,1)*skl(3,1)*c1111+ 646 & skl(3,1)*skl(3,1)*skl(2,2)*skl(3,2)*c1122+ 647 & skl(3,1)*skl(3,1)*skl(2,3)*skl(3,3)*c1122+ 648 & skl(3,1)*skl(3,2)*skl(2,1)*skl(3,2)*c1212+ 649 & skl(3,1)*skl(3,2)*skl(2,2)*skl(3,1)*c1212+ 650 & skl(3,1)*skl(3,3)*skl(2,1)*skl(3,3)*c1212+ 651 & skl(3,1)*skl(3,3)*skl(2,3)*skl(3,1)*c1212+ 652 & skl(3,2)*skl(3,1)*skl(2,1)*skl(3,2)*c1212+ 653 & skl(3,2)*skl(3,1)*skl(2,2)*skl(3,1)*c1212+ 654 & skl(3,2)*skl(3,2)*skl(2,1)*skl(3,1)*c1122+ 655 & skl(3,2)*skl(3,2)*skl(2,2)*skl(3,2)*c1111+ 656 & skl(3,2)*skl(3,2)*skl(2,3)*skl(3,3)*c1122+ 657 & skl(3,2)*skl(3,3)*skl(2,2)*skl(3,3)*c1212+ 658 & skl(3,2)*skl(3,3)*skl(2,3)*skl(3,2)*c1212+ 659 & skl(3,3)*skl(3,1)*skl(2,1)*skl(3,3)*c1212+ 660 & skl(3,3)*skl(3,1)*skl(2,3)*skl(3,1)*c1212+ 661 & skl(3,3)*skl(3,2)*skl(2,2)*skl(3,3)*c1212+ 662 & skl(3,3)*skl(3,2)*skl(2,3)*skl(3,2)*c1212+ 663 & skl(3,3)*skl(3,3)*skl(2,1)*skl(3,1)*c1122+ 664 & skl(3,3)*skl(3,3)*skl(2,2)*skl(3,2)*c1122+ 665 & skl(3,3)*skl(3,3)*skl(2,3)*skl(3,3)*c1111 666 elas(19)= 667 & skl(1,1)*skl(2,1)*skl(2,1)*skl(3,1)*c1111+ 668 & skl(1,1)*skl(2,1)*skl(2,2)*skl(3,2)*c1122+ 669 & skl(1,1)*skl(2,1)*skl(2,3)*skl(3,3)*c1122+ 670 & skl(1,1)*skl(2,2)*skl(2,1)*skl(3,2)*c1212+ 671 & skl(1,1)*skl(2,2)*skl(2,2)*skl(3,1)*c1212+ 672 & skl(1,1)*skl(2,3)*skl(2,1)*skl(3,3)*c1212+ 673 & skl(1,1)*skl(2,3)*skl(2,3)*skl(3,1)*c1212+ 674 & skl(1,2)*skl(2,1)*skl(2,1)*skl(3,2)*c1212+ 675 & skl(1,2)*skl(2,1)*skl(2,2)*skl(3,1)*c1212+ 676 & skl(1,2)*skl(2,2)*skl(2,1)*skl(3,1)*c1122+ 677 & skl(1,2)*skl(2,2)*skl(2,2)*skl(3,2)*c1111+ 678 & skl(1,2)*skl(2,2)*skl(2,3)*skl(3,3)*c1122+ 679 & skl(1,2)*skl(2,3)*skl(2,2)*skl(3,3)*c1212+ 680 & skl(1,2)*skl(2,3)*skl(2,3)*skl(3,2)*c1212+ 681 & skl(1,3)*skl(2,1)*skl(2,1)*skl(3,3)*c1212+ 682 & skl(1,3)*skl(2,1)*skl(2,3)*skl(3,1)*c1212+ 683 & skl(1,3)*skl(2,2)*skl(2,2)*skl(3,3)*c1212+ 684 & skl(1,3)*skl(2,2)*skl(2,3)*skl(3,2)*c1212+ 685 & skl(1,3)*skl(2,3)*skl(2,1)*skl(3,1)*c1122+ 686 & skl(1,3)*skl(2,3)*skl(2,2)*skl(3,2)*c1122+ 687 & skl(1,3)*skl(2,3)*skl(2,3)*skl(3,3)*c1111 688 elas(20)= 689 & skl(1,1)*skl(3,1)*skl(2,1)*skl(3,1)*c1111+ 690 & skl(1,1)*skl(3,1)*skl(2,2)*skl(3,2)*c1122+ 691 & skl(1,1)*skl(3,1)*skl(2,3)*skl(3,3)*c1122+ 692 & skl(1,1)*skl(3,2)*skl(2,1)*skl(3,2)*c1212+ 693 & skl(1,1)*skl(3,2)*skl(2,2)*skl(3,1)*c1212+ 694 & skl(1,1)*skl(3,3)*skl(2,1)*skl(3,3)*c1212+ 695 & skl(1,1)*skl(3,3)*skl(2,3)*skl(3,1)*c1212+ 696 & skl(1,2)*skl(3,1)*skl(2,1)*skl(3,2)*c1212+ 697 & skl(1,2)*skl(3,1)*skl(2,2)*skl(3,1)*c1212+ 698 & skl(1,2)*skl(3,2)*skl(2,1)*skl(3,1)*c1122+ 699 & skl(1,2)*skl(3,2)*skl(2,2)*skl(3,2)*c1111+ 700 & skl(1,2)*skl(3,2)*skl(2,3)*skl(3,3)*c1122+ 701 & skl(1,2)*skl(3,3)*skl(2,2)*skl(3,3)*c1212+ 702 & skl(1,2)*skl(3,3)*skl(2,3)*skl(3,2)*c1212+ 703 & skl(1,3)*skl(3,1)*skl(2,1)*skl(3,3)*c1212+ 704 & skl(1,3)*skl(3,1)*skl(2,3)*skl(3,1)*c1212+ 705 & skl(1,3)*skl(3,2)*skl(2,2)*skl(3,3)*c1212+ 706 & skl(1,3)*skl(3,2)*skl(2,3)*skl(3,2)*c1212+ 707 & skl(1,3)*skl(3,3)*skl(2,1)*skl(3,1)*c1122+ 708 & skl(1,3)*skl(3,3)*skl(2,2)*skl(3,2)*c1122+ 709 & skl(1,3)*skl(3,3)*skl(2,3)*skl(3,3)*c1111 710 elas(21)= 711 & skl(2,1)*skl(3,1)*skl(2,1)*skl(3,1)*c1111+ 712 & skl(2,1)*skl(3,1)*skl(2,2)*skl(3,2)*c1122+ 713 & skl(2,1)*skl(3,1)*skl(2,3)*skl(3,3)*c1122+ 714 & skl(2,1)*skl(3,2)*skl(2,1)*skl(3,2)*c1212+ 715 & skl(2,1)*skl(3,2)*skl(2,2)*skl(3,1)*c1212+ 716 & skl(2,1)*skl(3,3)*skl(2,1)*skl(3,3)*c1212+ 717 & skl(2,1)*skl(3,3)*skl(2,3)*skl(3,1)*c1212+ 718 & skl(2,2)*skl(3,1)*skl(2,1)*skl(3,2)*c1212+ 719 & skl(2,2)*skl(3,1)*skl(2,2)*skl(3,1)*c1212+ 720 & skl(2,2)*skl(3,2)*skl(2,1)*skl(3,1)*c1122+ 721 & skl(2,2)*skl(3,2)*skl(2,2)*skl(3,2)*c1111+ 722 & skl(2,2)*skl(3,2)*skl(2,3)*skl(3,3)*c1122+ 723 & skl(2,2)*skl(3,3)*skl(2,2)*skl(3,3)*c1212+ 724 & skl(2,2)*skl(3,3)*skl(2,3)*skl(3,2)*c1212+ 725 & skl(2,3)*skl(3,1)*skl(2,1)*skl(3,3)*c1212+ 726 & skl(2,3)*skl(3,1)*skl(2,3)*skl(3,1)*c1212+ 727 & skl(2,3)*skl(3,2)*skl(2,2)*skl(3,3)*c1212+ 728 & skl(2,3)*skl(3,2)*skl(2,3)*skl(3,2)*c1212+ 729 & skl(2,3)*skl(3,3)*skl(2,1)*skl(3,1)*c1122+ 730 & skl(2,3)*skl(3,3)*skl(2,2)*skl(3,2)*c1122+ 731 & skl(2,3)*skl(3,3)*skl(2,3)*skl(3,3)*c1111 732 endif 733! 734 do i=1,6 735 ep0(i)=xstateini(i,iint,iel) 736 enddo 737! 738! elastic strains 739! 740 do i=1,6 741 ee(i)=emec(i)-ep0(i) 742 enddo 743! 744! (visco)plastic constants: octahedral slip system 745! 746 do i=1,12 747 ck(i)=elconloc(4) 748 cn(i)=elconloc(5) 749 enddo 750! 751! (visco)plastic constants: cubic slip system 752! 753 do i=13,18 754 ck(i)=elconloc(6) 755 cn(i)=elconloc(7) 756 enddo 757! 758! global trial stress tensor 759! 760 if(iorien.gt.0) then 761 stri(1)=elas(1)*ee(1)+elas(2)*ee(2)+elas(4)*ee(3)+ 762 & 2.d0*(elas(7)*ee(4)+elas(11)*ee(5)+elas(16)*ee(6)) 763 & -beta(1) 764 stri(2)=elas(2)*ee(1)+elas(3)*ee(2)+elas(5)*ee(3)+ 765 & 2.d0*(elas(8)*ee(4)+elas(12)*ee(5)+elas(17)*ee(6)) 766 & -beta(2) 767 stri(3)=elas(4)*ee(1)+elas(5)*ee(2)+elas(6)*ee(3)+ 768 & 2.d0*(elas(9)*ee(4)+elas(13)*ee(5)+elas(18)*ee(6)) 769 & -beta(3) 770 stri(4)=elas(7)*ee(1)+elas(8)*ee(2)+elas(9)*ee(3)+ 771 & 2.d0*(elas(10)*ee(4)+elas(14)*ee(5)+elas(19)*ee(6)) 772 & -beta(4) 773 stri(5)=elas(11)*ee(1)+elas(12)*ee(2)+elas(13)*ee(3)+ 774 & 2.d0*(elas(14)*ee(4)+elas(15)*ee(5)+elas(20)*ee(6)) 775 & -beta(5) 776 stri(6)=elas(16)*ee(1)+elas(17)*ee(2)+elas(18)*ee(3)+ 777 & 2.d0*(elas(19)*ee(4)+elas(20)*ee(5)+elas(21)*ee(6)) 778 & -beta(6) 779 else 780 stri(1)=c1111*ee(1)+c1122*(ee(2)+ee(3))-beta(1) 781 stri(2)=c1111*ee(2)+c1122*(ee(1)+ee(3))-beta(2) 782 stri(3)=c1111*ee(3)+c1122*(ee(1)+ee(2))-beta(3) 783 stri(4)=2.d0*c1212*ee(4)-beta(4) 784 stri(5)=2.d0*c1212*ee(5)-beta(5) 785 stri(6)=2.d0*c1212*ee(6)-beta(6) 786 endif 787! 788! stress radius in each slip plane 789! 790 do i=1,18 791 sg(i)=xm(1,i)*stri(1)+xm(2,i)*stri(2)+xm(3,i)*stri(3)+ 792 & 2.d0*(xm(4,i)*stri(4)+xm(5,i)*stri(5)+xm(6,i)*stri(6)) 793 enddo 794! 795! evaluation of the yield surface 796! 797 do i=1,18 798 htri(i)=dabs(sg(i)) 799 enddo 800! 801 if(ielas.eq.1) then 802! 803! elastic stress 804! 805 do i=1,6 806 stre(i)=stri(i) 807 enddo 808! 809! elastic stiffness 810! 811 if(icmd.ne.3) then 812 if(iorien.gt.0) then 813 do i=1,21 814 stiff(i)=elas(i) 815 enddo 816 else 817 stiff(1)=c1111 818 stiff(2)=c1122 819 stiff(3)=c1111 820 stiff(4)=c1122 821 stiff(5)=c1122 822 stiff(6)=c1111 823 stiff(7)=0.d0 824 stiff(8)=0.d0 825 stiff(9)=0.d0 826 stiff(10)=c1212 827 stiff(11)=0.d0 828 stiff(12)=0.d0 829 stiff(13)=0.d0 830 stiff(14)=0.d0 831 stiff(15)=c1212 832 stiff(16)=0.d0 833 stiff(17)=0.d0 834 stiff(18)=0.d0 835 stiff(19)=0.d0 836 stiff(20)=0.d0 837 stiff(21)=c1212 838 endif 839! 840! updating the state variables 841! 842 do i=1,6 843 xstate(i,iint,iel)=ep0(i) 844 enddo 845 do i=7,24 846 xstate(i,iint,iel)=xstateini(i,iint,iel) 847 enddo 848 endif 849! 850 return 851 endif 852! 853! plastic deformation 854! 855 nrhs=1 856 lda=18 857 ldb=18 858! 859! initializing the state variables 860! 861 do i=1,6 862 ep(i)=ep0(i) 863 enddo 864 do i=1,18 865 dg(i)=0.d0 866c dg(i)=dtime*htri(i)**cn(i)/ck(i)**cn(i) 867 enddo 868! 869! major loop 870! 871 icounter=0 872 do 873 icounter=icounter+1 874 if(icounter.gt.100) then 875 pnewdt=0.25d0 876 return 877c write(*,*) 878c & '*ERROR in umat_single_crystal_creep: no convergence' 879c call exit(201) 880 endif 881! 882! elastic strains 883! 884 do i=1,6 885 ee(i)=emec(i)-ep(i) 886c write(*,*) 'umat_single_crystal_creep ee',icounter,i,ee(i) 887 enddo 888! 889! global trial stress tensor 890! 891 if(iorien.gt.0) then 892 stri(1)=elas(1)*ee(1)+elas(2)*ee(2)+elas(4)*ee(3)+ 893 & 2.d0*(elas(7)*ee(4)+elas(11)*ee(5)+elas(16)*ee(6)) 894 & -beta(1) 895 stri(2)=elas(2)*ee(1)+elas(3)*ee(2)+elas(5)*ee(3)+ 896 & 2.d0*(elas(8)*ee(4)+elas(12)*ee(5)+elas(17)*ee(6)) 897 & -beta(2) 898 stri(3)=elas(4)*ee(1)+elas(5)*ee(2)+elas(6)*ee(3)+ 899 & 2.d0*(elas(9)*ee(4)+elas(13)*ee(5)+elas(18)*ee(6)) 900 & -beta(3) 901 stri(4)=elas(7)*ee(1)+elas(8)*ee(2)+elas(9)*ee(3)+ 902 & 2.d0*(elas(10)*ee(4)+elas(14)*ee(5)+elas(19)*ee(6)) 903 & -beta(4) 904 stri(5)=elas(11)*ee(1)+elas(12)*ee(2)+elas(13)*ee(3)+ 905 & 2.d0*(elas(14)*ee(4)+elas(15)*ee(5)+elas(20)*ee(6)) 906 & -beta(5) 907 stri(6)=elas(16)*ee(1)+elas(17)*ee(2)+elas(18)*ee(3)+ 908 & 2.d0*(elas(19)*ee(4)+elas(20)*ee(5)+elas(21)*ee(6)) 909 & -beta(6) 910 else 911 stri(1)=c1111*ee(1)+c1122*(ee(2)+ee(3))-beta(1) 912 stri(2)=c1111*ee(2)+c1122*(ee(1)+ee(3))-beta(2) 913 stri(3)=c1111*ee(3)+c1122*(ee(1)+ee(2))-beta(3) 914 stri(4)=2.d0*c1212*ee(4)-beta(4) 915 stri(5)=2.d0*c1212*ee(5)-beta(5) 916 stri(6)=2.d0*c1212*ee(6)-beta(6) 917 endif 918 do i=1,6 919c write(*,*) 'umat_single_crystal_creep stri',icounter,i,stri(i) 920 enddo 921! 922! stress radius in each slip plane 923! 924 do i=1,18 925 sg(i)=xm(1,i)*stri(1)+xm(2,i)*stri(2)+xm(3,i)*stri(3)+ 926 & 2.d0*(xm(4,i)*stri(4)+xm(5,i)*stri(5)+xm(6,i)*stri(6)) 927c write(*,*) 'umat_single_crystal_creep sg',icounter,i,sg(i) 928 enddo 929! 930! evaluation of the yield surface 931! 932 do i=1,18 933c write(*,*) 'umat_single_crystal_creep ck,dtime,cn,dg',icounter, 934c & i,ck(i),dtime,cn(i),dg(i) 935 htri(i)=dabs(sg(i))-ck(i)*(dg(i)/dtime)**(1.d0/cn(i)) 936c write(*,*) 'umat_single_crystal_creep htri',icounter,i,htri(i) 937 enddo 938! 939! replace sg(i) by sgn(sg(i)) 940! 941 do i=1,18 942 if(sg(i).lt.0.d0) then 943 sg(i)=-1.d0 944 else 945 sg(i)=1.d0 946 endif 947 enddo 948! 949! determining the residual matrix 950! 951 do i=1,6 952 r(i)=ep0(i)-ep(i) 953 enddo 954 do i=1,18 955 do j=1,6 956 r(j)=r(j)+xm(j,i)*sg(i)*dg(i) 957 enddo 958 enddo 959! 960! check convergence 961! 962 if(icounter.gt.1) then 963 convergence=1 964 do i=1,18 965c if(dabs(htri(i)).gt.1.d-10) then 966 if(dabs(htri(i)).gt.1.d-3) then 967c write(*,*) 'divergence htri dg ', 968c & icounter,i,htri(i),dg(i) 969 convergence=0 970 exit 971 endif 972 enddo 973 if(convergence.eq.1) exit 974 endif 975c do i=1,18 976c if(dabs(htri(i)).gt.1.d-5) then 977c write(*,*) 'htri ',i,htri(i) 978c convergence=0 979c exit 980c endif 981c enddo 982c if(convergence.eq.1) then 983c dd=0.d0 984c do i=1,6 985c dd=dd+r(i)*r(i) 986c enddo 987c dd=sqrt(dd) 988c if(dd.gt.1.d-10) then 989c convergence=0 990c write(*,*) 'dd ',dd 991c else 992c exit 993c endif 994c endif 995! 996! compute xmc=c:xm 997! 998 do i=1,18 999 if(iorien.gt.0) then 1000 xmc(1,i)=elas(1)*xm(1,i)+elas(2)*xm(2,i)+ 1001 & elas(4)*xm(3,i)+2.d0*(elas(7)*xm(4,i)+ 1002 & elas(11)*xm(5,i)+elas(16)*xm(6,i)) 1003 xmc(2,i)=elas(2)*xm(1,i)+elas(3)*xm(2,i)+ 1004 & elas(5)*xm(3,i)+2.d0*(elas(8)*xm(4,i)+ 1005 & elas(12)*xm(5,i)+elas(17)*xm(6,i)) 1006 xmc(3,i)=elas(4)*xm(1,i)+elas(5)*xm(2,i)+ 1007 & elas(6)*xm(3,i)+2.d0*(elas(9)*xm(4,i)+ 1008 & elas(13)*xm(5,i)+elas(18)*xm(6,i)) 1009 xmc(4,i)=elas(7)*xm(1,i)+elas(8)*xm(2,i)+ 1010 & elas(9)*xm(3,i)+2.d0*(elas(10)*xm(4,i)+ 1011 & elas(14)*xm(5,i)+elas(19)*xm(6,i)) 1012 xmc(5,i)=elas(11)*xm(1,i)+elas(12)*xm(2,i)+ 1013 & elas(13)*xm(3,i)+2.d0*(elas(14)*xm(4,i)+ 1014 & elas(15)*xm(5,i)+elas(20)*xm(6,i)) 1015 xmc(6,i)=elas(16)*xm(1,i)+elas(17)*xm(2,i)+ 1016 & elas(18)*xm(3,i)+2.d0*(elas(19)*xm(4,i)+ 1017 & elas(20)*xm(5,i)+elas(21)*xm(6,i)) 1018 else 1019 xmc(1,i)=c1111*xm(1,i)+c1122*(xm(2,i)+xm(3,i)) 1020 xmc(2,i)=c1111*xm(2,i)+c1122*(xm(1,i)+xm(3,i)) 1021 xmc(3,i)=c1111*xm(3,i)+c1122*(xm(1,i)+xm(2,i)) 1022 xmc(4,i)=2.d0*c1212*xm(4,i) 1023 xmc(5,i)=2.d0*c1212*xm(5,i) 1024 xmc(6,i)=2.d0*c1212*xm(6,i) 1025 endif 1026 enddo 1027! 1028 neq=18 1029! 1030! filling the LHS 1031! 1032 do i=1,18 1033 do j=1,18 1034 if(i.ne.j) then 1035 gl(i,j)=(xm(1,i)*xmc(1,j)+ 1036 & xm(2,i)*xmc(2,j)+xm(3,i)*xmc(3,j)+2.d0* 1037 & (xm(4,i)*xmc(4,j)+xm(5,i)*xmc(5,j)+ 1038 & xm(6,i)*xmc(6,j))) 1039 & *sg(i)*sg(j) 1040 else 1041 gl(i,j)=(xm(1,i)*xmc(1,j)+ 1042 & xm(2,i)*xmc(2,j)+xm(3,i)*xmc(3,j)+2.d0* 1043 & (xm(4,i)*xmc(4,j)+xm(5,i)*xmc(5,j)+ 1044 & xm(6,i)*xmc(6,j))) 1045 endif 1046 enddo 1047 if(dg(i).gt.1.d-30) then 1048 gl(i,i)=gl(i,i)+ 1049 & (dg(i)/dtime)**(1.d0/cn(i)-1.d0)*ck(i)/ 1050 & (cn(i)*dtime) 1051 else 1052! 1053! for gamma ein default of 1.d-30 is taken to 1054! obtain a finite gradient 1055! 1056 gl(i,i)=gl(i,i)+ 1057 & (1.d-30/dtime)**(1.d0/cn(i)-1.d0)*ck(i)/ 1058 & (cn(i)*dtime) 1059 endif 1060 enddo 1061! 1062! filling the RHS 1063! 1064 do i=1,18 1065! 1066 gr(i,1)=htri(i) 1067cccccc & -ck(i)*(dg(i)/dtime)**(1.d0/cn(i)) 1068! 1069 do j=1,6 1070 t(j)=xmc(j,i)*sg(i) 1071 enddo 1072 do j=1,6 1073 gr(i,1)=gr(i,1)-t(j)*r(j) 1074 enddo 1075 gr(i,1)=gr(i,1) 1076 & -t(4)*r(4)-t(5)*r(5)-t(6)*r(6) 1077 enddo 1078! 1079! solve gl*ddg=gr 1080! 1081 call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info) 1082 if(info.ne.0) then 1083 write(*,*) '*ERROR in umat_single_crystal_creep:' 1084 write(*,*) ' linear equation solver exited' 1085 write(*,*) ' with error: info = ',info 1086 call exit(201) 1087 endif 1088! 1089 do i=1,18 1090 ddg(i)=gr(i,1) 1091 enddo 1092! 1093! updating the plastic strain 1094! 1095 do j=1,6 1096 ep(j)=ep(j)+r(j) 1097 do i=1,18 1098 ep(j)=ep(j)+xm(j,i)*sg(i)*ddg(i) 1099 enddo 1100c write(*,*) 'umat_single_crystal_creep ep',icounter,j,ep(j) 1101 enddo 1102! 1103! updating the consistency parameter 1104! 1105 do i=1,18 1106 dg(i)=max(dg(i)+ddg(i),0.d0) 1107c write(*,*) 'umat_single_crystal_creep ddg',icounter,i,ddg(i) 1108 enddo 1109c write(*,*) 'umat_single_crystal_creep ',icounter,dg(2) 1110! 1111! end of major loop 1112! 1113 enddo 1114! 1115! inversion of G 1116! 1117 do i=1,neq 1118 do j=1,neq 1119 gr(i,j)=0.d0 1120 enddo 1121 gr(i,i)=1.d0 1122 enddo 1123 nrhs=neq 1124 call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb,info) 1125 if(info.ne.0) then 1126 write(*,*) '*ERROR in sc.f: linear equation solver' 1127 write(*,*) ' exited with error: info = ',info 1128 call exit(201) 1129 endif 1130! 1131! storing the stress 1132! 1133 do i=1,6 1134 stre(i)=stri(i) 1135 enddo 1136! 1137! calculating the tangent stiffness matrix 1138! 1139 if(icmd.ne.3) then 1140 if(iorien.gt.0) then 1141 ddsdde(1,1)=elas(1) 1142 ddsdde(1,2)=elas(2) 1143 ddsdde(1,3)=elas(4) 1144 ddsdde(1,4)=elas(7) 1145 ddsdde(1,5)=elas(11) 1146 ddsdde(1,6)=elas(16) 1147c ddsdde(2,1)=elas(2) 1148 ddsdde(2,2)=elas(3) 1149 ddsdde(2,3)=elas(5) 1150 ddsdde(2,4)=elas(8) 1151 ddsdde(2,5)=elas(12) 1152 ddsdde(2,6)=elas(17) 1153c ddsdde(3,1)=elas(4) 1154c ddsdde(3,2)=elas(5) 1155 ddsdde(3,3)=elas(6) 1156 ddsdde(3,4)=elas(9) 1157 ddsdde(3,5)=elas(13) 1158 ddsdde(3,6)=elas(18) 1159c ddsdde(4,1)=elas(7) 1160c ddsdde(4,2)=elas(8) 1161c ddsdde(4,3)=elas(9) 1162 ddsdde(4,4)=elas(10) 1163 ddsdde(4,5)=elas(14) 1164 ddsdde(4,6)=elas(19) 1165c ddsdde(5,1)=elas(11) 1166c ddsdde(5,2)=elas(12) 1167c ddsdde(5,3)=elas(13) 1168c ddsdde(5,4)=elas(14) 1169 ddsdde(5,5)=elas(15) 1170 ddsdde(5,6)=elas(20) 1171c ddsdde(6,1)=elas(16) 1172c ddsdde(6,2)=elas(17) 1173c ddsdde(6,3)=elas(18) 1174c ddsdde(6,4)=elas(19) 1175c ddsdde(6,5)=elas(20) 1176 ddsdde(6,6)=elas(21) 1177 else 1178 do i=1,6 1179c do j=1,6 1180 do j=i,6 1181 ddsdde(i,j)=0.d0 1182 enddo 1183 enddo 1184 do i=1,3 1185 ddsdde(i,i)=c1111 1186 enddo 1187 do i=1,3 1188 do j=i+1,3 1189 ddsdde(i,j)=c1122 1190 enddo 1191c do j=1,i-1 1192c ddsdde(i,j)=c1122 1193c enddo 1194 ddsdde(i+3,i+3)=c1212 1195 enddo 1196 endif 1197 do i=1,18 1198 do j=1,18 1199 do k=1,6 1200c do l=1,6 1201 do l=k,6 1202 ddsdde(k,l)=ddsdde(k,l)- 1203 & gr(i,j)*xmc(k,i)*sg(i)*xmc(l,j)*sg(j) 1204 enddo 1205 enddo 1206 enddo 1207 enddo 1208! 1209! symmatrizing the stiffness matrix: to be changed: the 1210! matrix is symmetric! 1211! 1212 stiff(1)=ddsdde(1,1) 1213 stiff(2)=ddsdde(1,2) 1214 stiff(3)=ddsdde(2,2) 1215 stiff(4)=ddsdde(1,3) 1216 stiff(5)=ddsdde(2,3) 1217 stiff(6)=ddsdde(3,3) 1218 stiff(7)=ddsdde(1,4) 1219 stiff(8)=ddsdde(2,4) 1220 stiff(9)=ddsdde(3,4) 1221 stiff(10)=ddsdde(4,4) 1222 stiff(11)=ddsdde(1,5) 1223 stiff(12)=ddsdde(2,5) 1224 stiff(13)=ddsdde(3,5) 1225 stiff(14)=ddsdde(4,5) 1226 stiff(15)=ddsdde(5,5) 1227 stiff(16)=ddsdde(1,6) 1228 stiff(17)=ddsdde(2,6) 1229 stiff(18)=ddsdde(3,6) 1230 stiff(19)=ddsdde(4,6) 1231 stiff(20)=ddsdde(5,6) 1232 stiff(21)=ddsdde(6,6) 1233! 1234 endif 1235! 1236! updating the state variables 1237! 1238 do i=1,6 1239 xstate(i,iint,iel)=ep(i) 1240 enddo 1241 do i=7,24 1242 xstate(i,iint,iel)=xstateini(i,iint,iel)+dg(i-6) 1243 enddo 1244! 1245 return 1246 end 1247