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!     This 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 defplas(elconloc,elas,emec,ithermal,icmd,
20     &  beta,stre,ckl,vj)
21!
22!     calculates stiffness and stresses for the deformation plasticity
23!     material law
24!
25!     icmd=3: calculates stress at mechanical strain
26!     else: calculates stress and stiffness matrix at mechanical strain
27!
28      implicit none
29!
30      integer cauchy
31!
32      integer ithermal(*),icmd,i,j,k,l,m,n,ii,istart,iend,nt,kk(84)
33!
34      real*8 elconloc(*),elas(*),emec(*),beta(*),s(6),al,
35     &  ee,un,s0,xn,stre(*),eq,c0,c1,c2,c3,dkl(3,3),ekl(3,3),
36     &  q,dq,pp,el(6),ckl(3,3),vj
37!
38      kk=(/1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3,
39     &  1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3,3,3,1,3,
40     &  1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3,1,2,2,3,1,3,2,3,
41     &  2,3,2,3/)
42!
43      cauchy=1
44!
45      istart=1
46      iend=1
47!
48!     determining linear elastic material constants
49!
50      ee=elconloc(1)
51      un=elconloc(2)
52      s0=elconloc(3)
53      xn=elconloc(4)
54      al=elconloc(5)
55!
56      do i=1,6
57         el(i)=emec(i)
58      enddo
59!
60!     major loop
61!
62      do ii=istart,iend
63!
64         c0=(el(1)+el(2)+el(3))/3.d0
65!
66         el(1)=el(1)-c0
67         el(2)=el(2)-c0
68         el(3)=el(3)-c0
69!
70!        equivalent deviatoric strain
71!
72         eq=dsqrt(2.d0/3.d0*(el(1)*el(1)+el(2)*el(2)+
73     &                       el(3)*el(3)+2.d0*(el(4)*el(4)+
74     &                       el(5)*el(5)+el(6)*el(6))))
75!
76!        initial value of the Mises equivalent stress (q)
77!
78         c1=3.d0*ee*eq/(2.d0*(1.d0+un))
79!
80         if(c1.le.s0) then
81            q=c1
82         else
83            q=(s0**(xn-1)*ee*eq/al)**(1.d0/xn)
84         endif
85!
86!        determining the Mises equivalent stress q
87!
88         c1=2.d0*(1.d0+un)/3.d0
89         do
90            c2=al*(q/s0)**(xn-1.d0)
91            dq=(ee*eq-(c1+c2)*q)/(c1+xn*c2)
92            if((dabs(dq).lt.q*1.d-4).or.(dabs(dq).lt.1.d-10)) exit
93            q=q+dq
94         enddo
95!
96         if(icmd.ne.3) then
97!
98!           calculating the tangent stiffness matrix
99!
100!           initialization of the Delta Dirac function
101!
102            do i=1,3
103               do j=1,3
104                  dkl(i,j)=0.d0
105               enddo
106            enddo
107            do i=1,3
108               dkl(i,i)=1.d0
109            enddo
110!
111            ekl(1,1)=el(1)
112            ekl(2,2)=el(2)
113            ekl(3,3)=el(3)
114            ekl(1,2)=el(4)
115            ekl(1,3)=el(5)
116            ekl(2,3)=el(6)
117            ekl(2,1)=ekl(1,2)
118            ekl(3,1)=ekl(1,3)
119            ekl(3,2)=ekl(2,3)
120!
121            if(eq.lt.1.d-10) then
122               c1=ee/(1.d0+un)
123               c2=0.d0
124            else
125               c1=2.d0/(3.d0*eq)
126               c2=c1*(1.d0/eq-1.d0/(eq+(xn-1.d0)*c2*q/ee))
127               c1=c1*q
128            endif
129            c3=(ee/(1.d0-2.d0*un)-c1)/3.d0
130!
131            nt=0
132            do i=1,21
133               k=kk(nt+1)
134               l=kk(nt+2)
135               m=kk(nt+3)
136               n=kk(nt+4)
137               nt=nt+4
138               elas(i)=c1*((dkl(k,m)*dkl(l,n)+dkl(k,n)*dkl(l,m))/2.d0
139     &                     -c2*ekl(k,l)*ekl(m,n))
140     &                 +c3*dkl(k,l)*dkl(m,n)
141            enddo
142!
143!           conversion of the stiffness matrix from spatial coordinates
144!           coordinates into material coordinates
145!
146            call stiff2mat(elas,ckl,vj,cauchy)
147!
148         endif
149!
150!        calculating the stress
151!
152         pp=-ee*c0/(1.d0-2.d0*un)
153!
154         if(eq.lt.1.d-10) then
155            c1=0.d0
156         else
157            c1=2.d0*q/(3.d0*eq)
158         endif
159!
160         do i=1,6
161            s(i)=el(i)*c1
162         enddo
163         do i=1,3
164            s(i)=s(i)-pp
165         enddo
166!
167         do i=1,6
168            stre(i)=s(i)
169         enddo
170!
171!        converting the stress into the material frame of
172!        reference
173!
174         call str2mat(stre,ckl,vj,cauchy)
175!
176      enddo
177!
178      return
179      end
180