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 umat_lin_iso_el(amat,iel,iint,kode,elconloc,emec,emec0,
20     &        beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time,ttime,
21     &        icmd,ielas,mi,
22     &        nstate_,xstateini,xstate,stre,stiff,iorien,pgauss,orab)
23!
24!     calculates stiffness and stresses for a user defined material
25!     law
26!
27!     icmd=3: calculates stress at mechanical strain
28!     else: calculates stress at mechanical strain and the stiffness
29!           matrix
30!
31!     INPUT:
32!
33!     amat               material name
34!     iel                element number
35!     iint               integration point number
36!
37!     kode               material type (-100-#of constants entered
38!                        under *USER MATERIAL): can be used for materials
39!                        with varying number of constants
40!
41!     elconloc(21)       user defined constants defined by the keyword
42!                        card *USER MATERIAL (max. 21, actual # =
43!                        -kode-100), interpolated for the
44!                        actual temperature t1l
45!
46!     emec(6)            Lagrange mechanical strain tensor (component order:
47!                        11,22,33,12,13,23) at the end of the increment
48!                        (thermal strains are subtracted)
49!     emec0(6)           Lagrange mechanical strain tensor at the start of the
50!                        increment (thermal strains are subtracted)
51!     beta(6)            residual stress tensor (the stress entered under
52!                        the keyword *INITIAL CONDITIONS,TYPE=STRESS)
53!
54!     xokl(3,3)          deformation gradient at the start of the increment
55!     voj                Jacobian at the start of the increment
56!     xkl(3,3)           deformation gradient at the end of the increment
57!     vj                 Jacobian at the end of the increment
58!
59!     ithermal           0: no thermal effects are taken into account
60!                        >0: thermal effects are taken into account (triggered
61!                        by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE)
62!     t1l                temperature at the end of the increment
63!     dtime              time length of the increment
64!     time               step time at the end of the current increment
65!     ttime              total time at the start of the current step
66!
67!     icmd               not equal to 3: calculate stress and stiffness
68!                        3: calculate only stress
69!     ielas              0: no elastic iteration: irreversible effects
70!                        are allowed
71!                        1: elastic iteration, i.e. no irreversible
72!                           deformation allowed
73!
74!     mi(1)              max. # of integration points per element in the
75!                        model
76!     nstate_            max. # of state variables in the model
77!
78!     xstateini(nstate_,mi(1),# of elements)
79!                        state variables at the start of the increment
80!     xstate(nstate_,mi(1),# of elements)
81!                        state variables at the end of the increment
82!
83!     stre(6)            Piola-Kirchhoff stress of the second kind
84!                        at the start of the increment
85!
86!     iorien             number of the local coordinate axis system
87!                        in the integration point at stake (takes the value
88!                        0 if no local system applies)
89!     pgauss(3)          global coordinates of the integration point
90!     orab(7,*)          description of all local coordinate systems.
91!                        If a local coordinate system applies the global
92!                        tensors can be obtained by premultiplying the local
93!                        tensors with skl(3,3). skl is  determined by calling
94!                        the subroutine transformatrix:
95!                        call transformatrix(orab(1,iorien),pgauss,skl)
96!
97!
98!     OUTPUT:
99!
100!     xstate(nstate_,mi(1),# of elements)
101!                        updated state variables at the end of the increment
102!     stre(6)            Piola-Kirchhoff stress of the second kind at the
103!                        end of the increment
104!     stiff(21):         consistent tangent stiffness matrix in the material
105!                        frame of reference at the end of the increment. In
106!                        other words: the derivative of the PK2 stress with
107!                        respect to the Lagrangian strain tensor. The matrix
108!                        is supposed to be symmetric, only the upper half is
109!                        to be given in the same order as for a fully
110!                        anisotropic elastic material (*ELASTIC,TYPE=ANISO).
111!                        Notice that the matrix is an integral part of the
112!                        fourth order material tensor, i.e. the Voigt notation
113!                        is not used.
114!
115      implicit none
116!
117      character*80 amat
118!
119      integer ithermal(*),icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien,
120     &     i
121!
122      real*8 elconloc(*),stiff(21),emec(6),emec0(6),beta(6),stre(6),
123     &  vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*),
124     &  time,ttime
125!
126      real*8 xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
127!
128      real*8 e,un,al,um,am1,am2
129!
130!     insert here code to calculate the stresses
131!
132      e=elconloc(1)
133      un=elconloc(2)
134      al=un*e/(1.d0+un)/(1.d0-2.d0*un)
135      um=e/2.d0/(1.d0+un)
136      am1=al+2.d0*um
137      am2=2.d0*um
138!
139      stre(1)=am1*emec(1)+al*(emec(2)+emec(3))-beta(1)
140      stre(2)=am1*emec(2)+al*(emec(1)+emec(3))-beta(2)
141      stre(3)=am1*emec(3)+al*(emec(1)+emec(2))-beta(3)
142      stre(4)=am2*emec(4)-beta(4)
143      stre(5)=am2*emec(5)-beta(5)
144      stre(6)=am2*emec(6)-beta(6)
145c      write(*,*) 'stre',(stre(i),i=1,6)
146!
147      if(icmd.ne.3) then
148!
149!        insert here code to calculate the stiffness matrix
150!
151         stiff(1)=al+2.d0*um
152         stiff(2)=al
153         stiff(3)=al+2.d0*um
154         stiff(4)=al
155         stiff(5)=al
156         stiff(6)=al+2.d0*um
157         stiff(7)=0.d0
158         stiff(8)=0.d0
159         stiff(9)=0.d0
160         stiff(10)=um
161         stiff(11)=0.d0
162         stiff(12)=0.d0
163         stiff(13)=0.d0
164         stiff(14)=0.d0
165         stiff(15)=um
166         stiff(16)=0.d0
167         stiff(17)=0.d0
168         stiff(18)=0.d0
169         stiff(19)=0.d0
170         stiff(20)=0.d0
171         stiff(21)=um
172c      write(*,*) 'stiff',(stiff(i),i=1,21)
173      endif
174!
175      return
176      end
177