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!     Author: Sven Kassbohm
20!     Permission to use this routine was given by the author on
21!     January 17, 2015.
22!
23      subroutine umat_ciarlet_el(amat,iel,iint,kode,elconloc,emec,emec0,
24     &        beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time,ttime,
25     &        icmd,ielas,mi,
26     &        nstate_,xstateini,xstate,stre,stiff,iorien,pgauss,orab)
27!
28!     calculates stiffness and stresses for a user defined material
29!     law
30!
31!     icmd=3: calculates stress at mechanical strain
32!     else: calculates stress at mechanical strain and the stiffness
33!           matrix
34!
35!     INPUT:
36!
37!     amat               material name
38!     iel                element number
39!     iint               integration point number
40!
41!     kode               material type (-100-#of constants entered
42!                        under *USER MATERIAL): can be used for materials
43!                        with varying number of constants
44!
45!     elconloc(21)       user defined constants defined by the keyword
46!                        card *USER MATERIAL (max. 21, actual # =
47!                        -kode-100), interpolated for the
48!                        actual temperature t1l
49!
50!     emec(6)            Lagrange mechanical strain tensor (component order:
51!                        11,22,33,12,13,23) at the end of the increment
52!                        (thermal strains are subtracted)
53!     emec0(6)           Lagrange mechanical strain tensor at the start of the
54!                        increment (thermal strains are subtracted)
55!     beta(6)            residual stress tensor (the stress entered under
56!                        the keyword *INITIAL CONDITIONS,TYPE=STRESS)
57!
58!     xokl(3,3)          deformation gradient at the start of the increment
59!     voj                Jacobian at the start of the increment
60!     xkl(3,3)           deformation gradient at the end of the increment
61!     vj                 Jacobian at the end of the increment
62!
63!     ithermal           0: no thermal effects are taken into account
64!                        >0: thermal effects are taken into account (triggered
65!                        by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE)
66!     t1l                temperature at the end of the increment
67!     dtime              time length of the increment
68!     time               step time at the end of the current increment
69!     ttime              total time at the start of the current increment
70!
71!     icmd               not equal to 3: calculate stress and stiffness
72!                        3: calculate only stress
73!     ielas              0: no elastic iteration: irreversible effects
74!                        are allowed
75!                        1: elastic iteration, i.e. no irreversible
76!                           deformation allowed
77!
78!     mi(1)              max. # of integration points per element in the
79!                        model
80!     nstate_            max. # of state variables in the model
81!
82!     xstateini(nstate_,mi(1),# of elements)
83!                        state variables at the start of the increment
84!     xstate(nstate_,mi(1),# of elements)
85!                        state variables at the end of the increment
86!
87!     stre(6)            Piola-Kirchhoff stress of the second kind
88!                        at the start of the increment
89!
90!     iorien             number of the local coordinate axis system
91!                        in the integration point at stake (takes the value
92!                        0 if no local system applies)
93!     pgauss(3)          global coordinates of the integration point
94!     orab(7,*)          description of all local coordinate systems.
95!                        If a local coordinate system applies the global
96!                        tensors can be obtained by premultiplying the local
97!                        tensors with skl(3,3). skl is  determined by calling
98!                        the subroutine transformatrix:
99!                        call transformatrix(orab(1,iorien),pgauss,skl)
100!
101!
102!     OUTPUT:
103!
104!     xstate(nstate_,mi(1),# of elements)
105!                        updated state variables at the end of the increment
106!     stre(6)            Piola-Kirchhoff stress of the second kind at the
107!                        end of the increment
108!     stiff(21):         consistent tangent stiffness matrix in the material
109!                        frame of reference at the end of the increment. In
110!                        other words: the derivative of the PK2 stress with
111!                        respect to the Lagrangian strain tensor. The matrix
112!                        is supposed to be symmetric, only the upper half is
113!                        to be given in the same order as for a fully
114!                        anisotropic elastic material (*ELASTIC,TYPE=ANISO).
115!                        Notice that the matrix is an integral part of the
116!                        fourth order material tensor, i.e. the Voigt notation
117!                        is not used.
118!
119      implicit none
120!
121      character*80 amat
122      integer ithermal(*),icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien
123      real*8 elconloc(*),stiff(21),emec(6),emec0(6),beta(6),stre(6),
124     &  vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*),
125     &  time,ttime
126      real*8 xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
127!
128!     Ela-Pla-modifications:
129!
130      real*8 C1(3,3), C1i(3,3), detC1, dS2dE(3,3,3,3), S2PK1(3,3)
131      real*8 lamda, mu,un,e
132!
133!     change on 190315: input are Young's modulus and Poisson's
134!     coefficient instead of the Lame constants
135!
136      e=elconloc(1)
137      un=elconloc(2)
138      mu=e/(1.d0+un)
139      lamda=mu*un/(1.d0-2.d0*un)
140      mu=mu/2.d0
141
142c      lamda=elconloc(1)
143c      mu=elconloc(2)
144
145C     ciarlet_doc
146C
147C     Compute pullback C1(3,3) of metric tensor g.
148c     C1  = matmul(transpose(xkl),xkl)
149
150
151!
152!     change on 190315: the Cauchy tensor should be
153!     computed based on the mechanical Lagrange strain in
154!     order to discard the thermal strains; the deformation
155!     gradient still contains the thermal strain effect
156!
157!     right Cauchy-Green tensor (eloc contains the Lagrange strain,
158!     including thermal strain)
159!
160      C1(1,1)=2.d0*emec(1)+1.d0
161      C1(2,2)=2.d0*emec(2)+1.d0
162      C1(3,3)=2.d0*emec(3)+1.d0
163      C1(1,2)=2.d0*emec(4)
164      C1(2,1)=C1(1,2)
165      C1(1,3)=2.d0*emec(5)
166      C1(3,1)=C1(1,3)
167      C1(2,3)=2.d0*emec(6)
168      C1(3,2)=C1(2,3)
169C
170C     2PK-Stress at C:
171      call S2_Ciarlet(C1,lamda,mu, S2PK1,C1i,detC1)
172C     Tangent/linearization:
173      call dS2_Ciarlet_dE(lamda,mu,detC1,C1i,dS2dE)
174
175C     Dhondt-Voigt-Order, see umat_lin_iso_el.f
176      call voigt_33mat_to_6vec(S2PK1,stre)
177      call voigt_tetrad_to_matrix(dS2dE,stiff)
178
179C     ciarlet_cod
180
181      return
182      end
183
184
185
186      subroutine voigt_33mat_to_6vec(mat,vec)
187      implicit none
188C     in:
189      double precision mat(3,3)
190C     out:
191      double precision vec(6)
192C     11,22,33,12,13,23
193      vec(1)=mat(1,1)
194      vec(2)=mat(2,2)
195      vec(3)=mat(3,3)
196      vec(4)=mat(1,2)
197      vec(5)=mat(1,3)
198      vec(6)=mat(2,3)
199
200      return
201      end
202
203
204      subroutine voigt_tetrad_to_matrix(C,stiff)
205      implicit none
206C     in:
207      double precision C(3,3,3,3)
208C     out:
209      double precision stiff(21)
210
211C     indices according to
212C     1) file anisotropic.f and
213C     2) html-documentation
214C     -> Making C symmetric in the **last two** indices:
215C
216C     stiff stores:
217C     1111, 1122, 2222, 1133
218C     2233, 3333, 1112, 2212
219C
220C     1111
221      stiff(1)  = C(1,1,1,1)
222C     1122 = 2211
223      stiff(2)  = C(1,1,2,2)
224C     2222
225      stiff(3)  = C(2,2,2,2)
226C     1133 = 3311
227      stiff(4)  = C(1,1,3,3)
228C     2233 = 3322
229      stiff(5)  = C(2,2,3,3)
230C     3333
231      stiff(6)  = C(3,3,3,3)
232C     1112 = 1121 = 1211 = 2111
233      stiff(7)  = 0.5d0*(C(1,1,1,2)+C(1,1,2,1))
234C     2212 = 1222 = 2122 = 2221
235      stiff(8)  = 0.5d0*(C(2,2,1,2)+C(2,2,2,1))
236C
237C     stiff stores:
238C     3312, 1212, 1113, 2213
239C     3313, 1213, 1313, 1123
240C
241C     3312 = 1233 = 2133 = 3321
242      stiff(9)  = 0.5d0*(C(3,3,1,2)+C(3,3,2,1))
243C     1212 = 1221 = 2112 = 2121
244      stiff(10) = 0.5d0*(C(1,2,1,2)+C(1,2,2,1))
245C     1113 = 1131 = 1311 = 3111
246      stiff(11) = 0.5d0*(C(1,1,1,3)+C(1,1,3,1))
247C     2213 = 1322 = 2231 = 3122
248      stiff(12) = 0.5d0*(C(2,2,1,3)+C(2,2,3,1))
249C
250C     3313 = 1333 = 3133 = 3331
251      stiff(13) = 0.5d0*(C(3,3,1,3)+C(3,3,3,1))
252C     1213 = 1231 = 1312 = 1321 = 2113 = 2131 = 3112 = 3121:
253      stiff(14) = 0.5d0*(C(1,2,1,3)+C(1,2,3,1))
254C     1313 = 1331 = 3113 = 3131
255      stiff(15) = 0.5d0*(C(1,3,1,3)+C(1,3,3,1))
256C     1123 = 1132 = 2311 = 3211
257      stiff(16) = 0.5d0*(C(1,1,2,3)+C(1,1,3,2))
258C
259C     stiff stores:
260C     2223, 3323, 1223, 1323
261C     2323
262C
263C     2223 = 2232 = 2322 = 3222
264      stiff(17) = 0.5d0*(C(2,2,2,3)+C(2,2,3,2))
265C     3323 = 2333 = 3233 = 3332
266      stiff(18) = 0.5d0*(C(3,3,2,3)+C(3,3,3,2))
267C     1223 = 1232 = 2123 = 2132 = 2312 = 2321 = 3212 = 3221
268      stiff(19) = 0.5d0*(C(1,2,2,3)+C(1,2,3,2))
269C     1323 = 1332 = 2313 = 2331 = 3123 = 3132 = 3213 = 3231
270      stiff(20) = 0.5d0*(C(1,3,2,3)+C(1,3,3,2))
271C     2323 = 2332
272      stiff(21) = 0.5d0*(C(2,3,2,3)+C(2,3,3,2))
273
274      return
275      end
276
277      subroutine S2_Ciarlet(Cb,lamda,mu, S2,Cbi,detC)
278C     computes 2PK stresses for Ciarlet-model
279C
280      implicit none
281C     input:
282      double precision Cb(3,3),lamda,mu
283C     output:
284      double precision S2(3,3), Cbi(3,3), detC
285C     local:
286      double precision id(3,3)
287      double precision f1
288      integer I,J
289      logical OK_FLAG
290
291C     Set id:
292      do I=1,3
293         do J=1,3
294            id(I,J)=0.d0
295         enddo
296         id(I,I)=1.d0
297      enddo
298
299C     Inverse of Cb and det(C):
300      call M33INV_DET(Cb, Cbi, detC, OK_FLAG)
301C
302C
303      f1 = 0.5d0*lamda*(detC-1.d0)
304      do I=1,3
305         do J=1,3
306            S2(I,J)=f1*Cbi(I,J) + mu * ( id(I,J)-Cbi(I,J) )
307         enddo
308      enddo
309      return
310      end
311
312      subroutine dS2_Ciarlet_dE(lamda,mu,detC,Cbi,A)
313C     computes linearization of 2PK stresses with E
314C     for Ciarlet-model
315      implicit none
316C     input:
317      double precision lamda,mu,detC,Cbi(3,3)
318C     output:
319      double precision A(3,3,3,3)
320C     local:
321      double precision f1,f2
322      integer I,J,K,L
323
324      f1 = lamda*detC
325      f2 = 2.d0*mu - lamda*(detC-1.d0)
326      do I=1,3
327         do J=1,3
328            do K=1,3
329               do L=1,3
330                  A(I,J,K,L) =
331     &                 f1 * Cbi(K,L)*Cbi(I,J) +
332     &                 f2 * Cbi(I,K)*Cbi(L,J)
333               enddo
334            enddo
335         enddo
336      enddo
337      return
338      end
339
340      subroutine M33INV_DET(A, AINV, DET, OK_FLAG)
341C
342C
343C     !******************************************************************
344C     !  M33INV_DET  -  Compute the inverse of a 3x3 matrix.
345C     !
346C     !  A       = input 3x3 matrix to be inverted
347C     !  AINV    = output 3x3 inverse of matrix A
348C     !  OK_FLAG = (output) .TRUE. if the input matrix could be inverted,
349C     !            and .FALSE. if the input matrix is singular.
350C     !******************************************************************
351
352      IMPLICIT NONE
353
354      DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN)  :: A
355      DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: AINV
356      LOGICAL, INTENT(OUT) :: OK_FLAG
357
358      DOUBLE PRECISION, PARAMETER :: EPS = 1.0D-10
359      DOUBLE PRECISION :: DET
360      DOUBLE PRECISION, DIMENSION(3,3) :: COFACTOR
361
362
363      DET =   A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)
364     &      - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)
365     &      + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)
366
367      IF (ABS(DET) .LE. EPS) THEN
368         AINV = 0.0D0
369         OK_FLAG = .FALSE.
370         RETURN
371      END IF
372
373      COFACTOR(1,1) = +(A(2,2)*A(3,3)-A(2,3)*A(3,2))
374      COFACTOR(1,2) = -(A(2,1)*A(3,3)-A(2,3)*A(3,1))
375      COFACTOR(1,3) = +(A(2,1)*A(3,2)-A(2,2)*A(3,1))
376      COFACTOR(2,1) = -(A(1,2)*A(3,3)-A(1,3)*A(3,2))
377      COFACTOR(2,2) = +(A(1,1)*A(3,3)-A(1,3)*A(3,1))
378      COFACTOR(2,3) = -(A(1,1)*A(3,2)-A(1,2)*A(3,1))
379      COFACTOR(3,1) = +(A(1,2)*A(2,3)-A(1,3)*A(2,2))
380      COFACTOR(3,2) = -(A(1,1)*A(2,3)-A(1,3)*A(2,1))
381      COFACTOR(3,3) = +(A(1,1)*A(2,2)-A(1,2)*A(2,1))
382
383      AINV = TRANSPOSE(COFACTOR) / DET
384
385      OK_FLAG = .TRUE.
386
387      RETURN
388      END
389
390
391