FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
calMatMatrix.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
6 module m_matmatrix
7  use hecmw_util
8  use mmaterial
9  use mmechgauss
10  use m_elasticlinear
11  use mhyperelastic
12  use m_elastoplastic
13  use mviscoelastic
14  use mcreep
15  use muelastic
16  use mumat
17 
18  implicit none
19 
20 contains
21 
23  integer function getnlgeomflag( gauss )
24  type( tgaussstatus ), intent(in) :: gauss
25  getnlgeomflag = gauss%pMaterial%nlgeom_flag
26  end function
27 
29  subroutine matlmatrix( gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp )
30  type( tgaussstatus ), intent(in) :: gauss
31  integer, intent(in) :: sectType
32  real(kind=kreal), intent(out) :: matrix(:,:)
33  real(kind=kreal), intent(in) :: time
34  real(kind=kreal), intent(in) :: dtime
35  real(kind=kreal), intent(in) :: cdsys(3,3)
36  real(kind=kreal), intent(in), optional :: temperature
37  integer(kind=kint), intent(in), optional :: isEp
38 
39  integer :: i
40  integer :: flag
41  real(kind=kreal) :: cijkl(3,3,3,3)
42  type( tmaterial ), pointer :: matl
43  matl=>gauss%pMaterial
44 
45  flag = 0
46  if( present(isep) )then
47  if( isep == 1 )flag = 1
48  endif
49 
50  if( matl%mtype==userelastic ) then
51  call uelasticmatrix( matl%variables(101:), gauss%strain, matrix )
52  elseif( isviscoelastic(matl%mtype) ) then
53  if( present(temperature) ) then
54  call calviscoelasticmatrix( matl, secttype, dtime, matrix, temperature )
55  else
56  call calviscoelasticmatrix( matl, secttype, dtime, matrix )
57  endif
58  elseif( iselastic(matl%mtype) .or. flag==1 ) then
59  if(flag==1)then
61  else
62  i = getelastictype(gauss%pMaterial%mtype)
63  endif
64 
65  if( i==0 ) then
66  if( present(temperature) ) then
67  call calelasticmatrix( matl, secttype, matrix, temperature )
68  else
69  call calelasticmatrix( matl, secttype, matrix )
70  endif
71  elseif( i==1 ) then
72  if( present(temperature) ) then
73  call calelasticmatrix_ortho( gauss%pMaterial, secttype, cdsys, matrix, temperature )
74  else
75  call calelasticmatrix_ortho( gauss%pMaterial, secttype, cdsys, matrix )
76  endif
77  else
78  print *, "Elasticity type", matl%mtype, "not supported"
79  stop
80  endif
81 
82  elseif( matl%mtype==neohooke .or. matl%mtype==mooneyrivlin ) then
83  call calelasticmooneyrivlin( matl, secttype, cijkl, gauss%strain )
84  call mat_c2d( cijkl, matrix, secttype )
85  elseif( matl%mtype==arrudaboyce ) then
86  call calelasticarrudaboyce( matl, secttype, cijkl, gauss%strain )
87  call mat_c2d( cijkl, matrix, secttype )
88  elseif( matl%mtype==mooneyrivlin_aniso ) then
89  call calelasticmooneyrivlinaniso( matl, secttype, cijkl, gauss%strain, cdsys )
90  call mat_c2d( cijkl, matrix, secttype )
91  elseif( matl%mtype==userhyperelastic ) then
92  call uelasticmatrix( matl%variables(101:), gauss%strain, matrix )
93  elseif( iselastoplastic(matl%mtype) ) then
94  if( present( temperature ) ) then
95  call calelastoplasticmatrix( matl, secttype, gauss%stress, &
96  gauss%istatus(1), gauss%fstatus, gauss%plstrain, matrix, temperature )
97  else
98  call calelastoplasticmatrix( matl, secttype, gauss%stress, &
99  gauss%istatus(1), gauss%fstatus, gauss%plstrain, matrix )
100  endif
101  elseif( matl%mtype==usermaterial ) then
102  call umatlmatrix( matl%name, matl%variables(101:), gauss%strain, &
103  gauss%stress, gauss%fstatus, matrix, dtime, time )
104  elseif( matl%mtype==norton ) then
105  if( present( temperature ) ) then
106  call iso_creep( matl, secttype, gauss%stress, gauss%strain, gauss%fstatus, &
107  gauss%plstrain, dtime, time, matrix, temperature )
108  else
109  call iso_creep( matl, secttype, gauss%stress, gauss%strain, gauss%fstatus, &
110  gauss%plstrain, dtime, time, matrix )
111  endif
112  else
113  stop "Material type not supported!"
114  endif
115 
116  end subroutine
117 
118  !
120  subroutine stressupdate( gauss, sectType, strain, stress, cdsys, time, dtime, temp, tempn )
121  type( tgaussstatus ), intent(inout) :: gauss
122  integer, intent(in) :: sectType
123  real(kind=kreal), intent(in) :: strain(6)
124  real(kind=kreal), intent(out) :: stress(6)
125  real(kind=kreal), intent(in) :: cdsys(3,3)
126  real(kind=kreal), intent(in), optional :: time
127  real(kind=kreal), intent(in), optional :: dtime
128  real(kind=kreal), optional :: temp
129  real(kind=kreal), optional :: tempn
130 
131  if( gauss%pMaterial%mtype==neohooke .or. gauss%pMaterial%mtype==mooneyrivlin ) then
132  call calupdateelasticmooneyrivlin( gauss%pMaterial, secttype, strain, stress )
133  elseif( gauss%pMaterial%mtype==arrudaboyce ) then ! Arruda-Boyce Hyperelastic material
134  call calupdateelasticarrudaboyce( gauss%pMaterial, secttype, strain, stress )
135  elseif( gauss%pMaterial%mtype==mooneyrivlin_aniso ) then
136  call calupdateelasticmooneyrivlinaniso( gauss%pMaterial, secttype, strain, stress, cdsys )
137  elseif( gauss%pMaterial%mtype==userhyperelastic .or. gauss%pMaterial%mtype==userelastic ) then ! user-defined
138  call uelasticupdate( gauss%pMaterial%variables(101:), strain, stress )
139  elseif( isviscoelastic( gauss%pMaterial%mtype) ) then
140  if( .not. present(dtime) ) stop "error in viscoelastic update!"
141  if( present(temp) .and. present(tempn) ) then
142  call updateviscoelastic( gauss%pMaterial, secttype, strain, stress, gauss%fstatus, dtime, temp, tempn )
143  else
144  call updateviscoelastic( gauss%pMaterial, secttype, strain, stress, gauss%fstatus, dtime )
145  endif
146  elseif ( gauss%pMaterial%mtype==norton ) then
147  if( .not. present(dtime) ) stop "error in viscoelastic update!"
148  if( present(temp) ) then
149  call update_iso_creep( gauss%pMaterial, secttype, strain, stress, gauss%fstatus, gauss%plstrain, dtime, time, temp )
150  else
151  call update_iso_creep( gauss%pMaterial, secttype, strain, stress, gauss%fstatus, gauss%plstrain, dtime, time )
152  endif
153  elseif ( gauss%pMaterial%mtype==usermaterial) then ! user-defined
154  call uupdate( gauss%pMaterial%name, gauss%pMaterial%variables(101:), &
155  strain, stress, gauss%fstatus, dtime, time )
156  end if
157 
158  end subroutine stressupdate
159 
161  subroutine mat_c2d( cijkl, dij, itype )
162  real(kind=kreal), intent(in) :: cijkl(3,3,3,3)
163  real(kind=kreal), intent(out) :: dij(6,6)
164  integer, intent(in) :: itype
165 
166  dij(:,:) = 0.d0
167  select case( itype )
168  case( d3 )
169  dij(1,1) = cijkl(1,1,1,1) ! ---
170  dij(1,2) = cijkl(1,1,2,2)
171  dij(1,3) = cijkl(1,1,3,3)
172  dij(1,4) = cijkl(1,1,1,2)
173  dij(1,5) = cijkl(1,1,2,3)
174  dij(1,6) = cijkl(1,1,3,1)
175  dij(2,1) = cijkl(2,2,1,1) ! ---
176  dij(2,2) = cijkl(2,2,2,2)
177  dij(2,3) = cijkl(2,2,3,3)
178  dij(2,4) = cijkl(2,2,1,2)
179  dij(2,5) = cijkl(2,2,2,3)
180  dij(2,6) = cijkl(2,2,3,1)
181  dij(3,1) = cijkl(3,3,1,1) ! ---
182  dij(3,2) = cijkl(3,3,2,2)
183  dij(3,3) = cijkl(3,3,3,3)
184  dij(3,4) = cijkl(3,3,1,2)
185  dij(3,5) = cijkl(3,3,2,3)
186  dij(3,6) = cijkl(3,3,3,1)
187  dij(4,1) = cijkl(1,2,1,1) ! ---
188  dij(4,2) = cijkl(1,2,2,2)
189  dij(4,3) = cijkl(1,2,3,3)
190  dij(4,4) = cijkl(1,2,1,2)
191  dij(4,5) = cijkl(1,2,2,3)
192  dij(4,6) = cijkl(1,2,3,1)
193  dij(5,1) = cijkl(2,3,1,1) ! ---
194  dij(5,2) = cijkl(2,3,2,2)
195  dij(5,3) = cijkl(2,3,3,3)
196  dij(5,4) = cijkl(2,3,1,2)
197  dij(5,5) = cijkl(2,3,2,3)
198  dij(5,6) = cijkl(2,3,3,1)
199  dij(6,1) = cijkl(3,1,1,1) ! ---
200  dij(6,2) = cijkl(3,1,2,2)
201  dij(6,3) = cijkl(3,1,3,3)
202  dij(6,4) = cijkl(3,1,1,2)
203  dij(6,5) = cijkl(3,1,2,3)
204  dij(6,6) = cijkl(3,1,3,1)
205  !
206  case( planestress, planestrain )
207  dij(1,1) = cijkl(1,1,1,1) ! ---
208  dij(1,2) = cijkl(1,1,2,2)
209  dij(1,3) = cijkl(1,1,1,2)
210  dij(2,1) = cijkl(2,2,1,1) ! ---
211  dij(2,2) = cijkl(2,2,2,2)
212  dij(2,3) = cijkl(2,2,1,2)
213  dij(3,1) = cijkl(1,2,1,1) ! ---
214  dij(3,2) = cijkl(1,2,2,2)
215  dij(3,3) = cijkl(1,2,1,2)
216  case( axissymetric )
217  dij(1,1) = cijkl(1,1,1,1)
218  dij(1,2) = cijkl(1,1,2,2)
219  dij(1,3) = cijkl(1,1,1,2)
220  dij(1,4) = cijkl(1,1,3,3)
221  dij(2,1) = cijkl(2,2,1,1)
222  dij(2,2) = cijkl(2,2,2,2)
223  dij(2,3) = cijkl(2,2,1,2)
224  dij(2,4) = cijkl(2,2,3,3)
225  dij(3,1) = cijkl(1,2,1,1)
226  dij(3,2) = cijkl(1,2,2,2)
227  dij(3,3) = cijkl(1,2,1,2)
228  dij(3,4) = cijkl(1,2,3,3)
229  dij(4,1) = cijkl(3,3,1,1)
230  dij(4,2) = cijkl(3,3,2,2)
231  dij(4,3) = cijkl(3,3,1,2)
232  dij(4,4) = cijkl(3,3,3,3)
233  case( shell )
234  end select
235 
236  end subroutine mat_c2d
237 
238  ! (Gaku Hashimoto, The University of Tokyo, 2012/11/15) <
239  !####################################################################
240  subroutine matlmatrix_shell &
241  (gauss, secttype, d, &
242  e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
243  alpha, n_layer)
244  !####################################################################
245 
246  type(tgaussstatus), intent(in) :: gauss
247  integer, intent(in) :: sectType, n_layer
248  real(kind = kreal), intent(out) :: d(:, :)
249  real(kind = kreal), intent(in) :: e1_hat(3), e2_hat(3), e3_hat(3)
250  real(kind = kreal), intent(in) :: cg1(3), cg2(3), cg3(3)
251  real(kind = kreal), intent(out) :: alpha
252 
253  !--------------------------------------------------------------------
254 
255  real(kind = kreal) :: c(3, 3, 3, 3)
256  type(tmaterial), pointer :: matl
257 
258  !--------------------------------------------------------------------
259 
260  matl => gauss%pMaterial
261 
262  !--------------------------------------------------------------------
263 
264  if( iselastic(matl%mtype) ) then
265  call linearelastic_shell &
266  (matl, secttype, c, &
267  e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
268  alpha, n_layer)
269 
270  call mat_c2d_shell(c, d, secttype)
271  else
272  stop "Material type not supported!"
273  end if
274 
275  !--------------------------------------------------------------------
276 
277  return
278 
279  !####################################################################
280  end subroutine matlmatrix_shell
281  !####################################################################
282  ! > (Gaku Hashimoto, The University of Tokyo, 2012/11/15)
283 
284 
285  ! (Gaku Hashimoto, The University of Tokyo, 2012/11/15) <
286  !####################################################################
287  subroutine mat_c2d_shell(c, D, itype)
288  !####################################################################
289 
290  real(kind = kreal), intent(in) :: c(:, :, :, :)
291  real(kind = kreal), intent(out) :: d(:, :)
292  integer, intent(in) :: itype
293 
294  !--------------------------------------------------------------------
295 
296  integer :: index_i(5), index_j(5), &
297  index_k(5), index_l(5)
298  integer :: i, j, k, l
299  integer :: is, js
300 
301  !--------------------------------------------------------------------
302 
303  index_i(1) = 1
304  index_i(2) = 2
305  index_i(3) = 1
306  index_i(4) = 2
307  index_i(5) = 3
308 
309  index_j(1) = 1
310  index_j(2) = 2
311  index_j(3) = 2
312  index_j(4) = 3
313  index_j(5) = 1
314 
315  index_k(1) = 1
316  index_k(2) = 2
317  index_k(3) = 1
318  index_k(4) = 2
319  index_k(5) = 3
320 
321  index_l(1) = 1
322  index_l(2) = 2
323  index_l(3) = 2
324  index_l(4) = 3
325  index_l(5) = 1
326 
327  !--------------------------------------------------------------------
328 
329  d(:, :) = 0.0d0
330 
331  !--------------------------------------------------------------------
332 
333  select case( itype )
334  case( shell )
335 
336  do js = 1, 5
337 
338  do is = 1, 5
339 
340  i = index_i(is)
341  j = index_j(is)
342  k = index_k(js)
343  l = index_l(js)
344 
345  d(is, js) = c(i, j, k, l)
346 
347  end do
348 
349  end do
350 
351  end select
352 
353  !--------------------------------------------------------------------
354 
355  return
356 
357  !####################################################################
358  end subroutine mat_c2d_shell
359  !####################################################################
360  ! > (Gaku Hashimoto, The University of Tokyo, 2012/11/15)
361 
362 end module m_matmatrix
I/O and Utility.
Definition: hecmw_util_f.F90:7
This module provides functions for elastic material.
subroutine linearelastic_shell(matl, sectType, c, e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, alpha, n_layer)
subroutine calelasticmatrix_ortho(matl, sectType, bij, DMAT, temp)
Calculate orthotropic elastic matrix.
subroutine calelasticmatrix(matl, sectType, D, temp)
Calculate isotropic elastic matrix.
This module provide functions for elastoplastic calculation.
subroutine calelastoplasticmatrix(matl, sectType, stress, istat, extval, plstrain, D, temperature)
This subroutine calculates elastoplastic constitutive relation.
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp)
Calculate constituive matrix.
subroutine stressupdate(gauss, sectType, strain, stress, cdsys, time, dtime, temp, tempn)
Update strain and stress for elastic and hyperelastic materials.
subroutine mat_c2d_shell(c, D, itype)
integer function getnlgeomflag(gauss)
Fetch the nlgeom flag of the material.
subroutine matlmatrix_shell(gauss, sectType, D, e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, alpha, n_layer)
subroutine mat_c2d(cijkl, dij, itype)
Transfer rank 4 constituive matrix to rank 2 form.
This module provides functions for creep calculation.
Definition: creep.f90:6
subroutine iso_creep(matl, sectType, stress, strain, extval, plstrain, dtime, ttime, stiffness, temp)
This subroutine calculates stiffness for elastically isotropic materials with isotropic creep.
Definition: creep.f90:19
subroutine update_iso_creep(matl, sectType, strain, stress, extval, plstrain, dtime, ttime, temp)
This subroutine calculates stresses and creep status for an elastically isotropic material with isotr...
Definition: creep.f90:122
This module provides functions for hyperelastic calculation.
Definition: Hyperelastic.f90:6
subroutine calelasticmooneyrivlin(matl, sectType, cijkl, strain)
This subroutine provides elastic tangent coefficient for Mooney-Rivlin hyperelastic material.
subroutine calupdateelasticmooneyrivlinaniso(matl, sectType, strain, stress, cdsys)
This subroutine provides to update stress and strain for anisotropic Mooney-Rivlin material.
subroutine calupdateelasticmooneyrivlin(matl, sectType, strain, stress)
This subroutine provides to update stress and strain for Mooney-Rivlin material.
subroutine calupdateelasticarrudaboyce(matl, sectType, dstrain, dstress)
This subroutine provides to update stress and strain for Arrude-Royce material.
subroutine calelasticarrudaboyce(matl, sectType, cijkl, strain)
This subroutine provides elastic tangent coefficient for Arruda-Boyce hyperelastic material.
subroutine calelasticmooneyrivlinaniso(matl, sectType, cijkl, strain, cdsys)
This subroutine provides elastic tangent coefficient for anisotropic Mooney-Rivlin hyperelastic mater...
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter mooneyrivlin
Definition: material.f90:65
integer(kind=kint), parameter mooneyrivlin_aniso
Definition: material.f90:68
integer(kind=kint), parameter planestress
Definition: material.f90:77
integer function getelastictype(mtype)
Get elastic type.
Definition: material.f90:282
integer(kind=kint), parameter d3
Definition: material.f90:76
integer(kind=kint), parameter planestrain
Definition: material.f90:78
integer(kind=kint), parameter elastic
Definition: material.f90:58
integer(kind=kint), parameter arrudaboyce
Definition: material.f90:66
integer(kind=kint), parameter shell
Definition: material.f90:80
integer(kind=kint), parameter norton
Definition: material.f90:71
integer(kind=kint), parameter axissymetric
Definition: material.f90:79
integer(kind=kint), parameter userelastic
Definition: material.f90:60
integer(kind=kint), parameter neohooke
Definition: material.f90:64
integer(kind=kint), parameter usermaterial
Definition: material.f90:56
logical function iselastic(mtype)
If it is an elastic material?
Definition: material.f90:327
integer(kind=kint), parameter userhyperelastic
Definition: material.f90:67
logical function isviscoelastic(mtype)
If it is an viscoelastic material?
Definition: material.f90:354
logical function iselastoplastic(mtype)
If it is an elastoplastic material?
Definition: material.f90:336
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
This module provides interface for elastic or hyperelastic calculation.
Definition: uelastic.f90:7
subroutine uelasticmatrix(matl, strain, D)
This subroutine calculates constitutive relation.
Definition: uelastic.f90:13
subroutine uelasticupdate(matl, strain, stress)
This subroutine calculate updated strain and stress.
Definition: uelastic.f90:41
This subroutine read in used-defined material properties tangent.
Definition: umat.f90:7
subroutine uupdate(mname, matl, strain, stress, fstat, dtime, ttime, temperature)
This subroutine calculate strain and stress increment.
Definition: umat.f90:31
subroutine umatlmatrix(mname, matl, strain, stress, fstat, D, dtime, ttime, temperature)
This subroutine calculates constitutive matrix.
Definition: umat.f90:17
This module provides functions for viscoelastic calculation.
Definition: Viscoelastic.f90:6
subroutine updateviscoelastic(matl, sectType, eps, sig, vsig, dt, temp, tempn)
This subroutine provides to update stress for viscoelastic material.
subroutine calviscoelasticmatrix(matl, sectType, dt, D, temp)
This subroutine calculates tangent moduli for isotropic viscoelastic material.
Stucture to management all material relates data.
Definition: material.f90:144
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13