1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5!> This module summarizes all infomation of material properties
6module mMaterial
7  use hecmw_util
8  use m_table
9  use Table_DICTS
10  implicit none
11
12  ! Following algorithm type
13  integer(kind=kint), parameter :: INFINITESIMAL = 0
14  integer(kind=kint), parameter :: TOTALLAG  = 1
15  integer(kind=kint), parameter :: UPDATELAG = 2
16
17  ! Following material types. All material type number consists with integer of six digits.
18  !   First digit: Indicates physical type
19  !                1: mechanical deformation analysis
20  !                2: heat conduct analysis
21  !                ......
22  !   Second digit:
23  !     Mechanical analysis
24  !     1: Elastic
25  !     2: Elastoplastic
26  !     3: Hyperelastic
27  !     4: Viscoelastic
28  !     5: Viscoplastic
29  !     Heat conductiovity
30  !     ......
31  !   Third digit:
32  !     For elastic or elastoplastic deformation, elastic
33  !       0: isotropic                      ie. 110000
34  !       1: with transever anisotropity        111000
35  !     For hyperelastic deformation
36  !       0: Neo-Hooke                          130000
37  !       1: Mooney-Rivlin                      131000
38  !       2: Arruda-Boyce                       132000
39  !   Fourth digit
40  !     For elastoplastic problem,yield function
41  !       0: isotropic (Mises)
42  !       1: Mohr-Coulomb
43  !       2: Drucker-Prager
44  !   Fifth digit:
45  !     For elastoplastic deformation, hardening law
46  !       0: Linear hardening           i.e.  120000
47  !       1: Multilinear hardening            120010
48  !       2: Swift
49  !       3: Ramberg-Osgood
50  !       4: linear kinematic
51  !       5: combined (linear kinematic + linear isotropic)
52  !   Six digit:
53  !     For visco-elastoplastic deformation, visco law
54  !       0: Norton                     i.e.  150000
55  !       1: Striab                           150001
56  integer(kind=kint), parameter :: USERMATERIAL          = 100000
57
58  integer(kind=kint), parameter :: ELASTIC               = 110000
59  integer(kind=kint), parameter :: MN_ORTHOELASTIC       = 111000
60  integer(kind=kint), parameter :: USERELASTIC           = 112000
61
62  integer(kind=kint), parameter :: EPLASTIC              = 120000
63
64  integer(kind=kint), parameter :: NEOHOOKE              = 130000
65  integer(kind=kint), parameter :: MOONEYRIVLIN          = 131000
66  integer(kind=kint), parameter :: ARRUDABOYCE           = 132000
67  integer(kind=kint), parameter :: USERHYPERELASTIC      = 133000
68  integer(kind=kint), parameter :: MOONEYRIVLIN_ANISO    = 134000
69
70  integer(kind=kint), parameter :: VISCOELASTIC          = 140000
71  integer(kind=kint), parameter :: NORTON                = 150000
72
73  integer(kind=kint), parameter :: INCOMP_NEWTONIAN      = 160000
74
75  ! Following section type
76  integer(kind=kint), parameter :: D3            = -1
77  integer(kind=kint), parameter :: PlaneStress   = 1
78  integer(kind=kint), parameter :: PlaneStrain   = 0
79  integer(kind=kint), parameter :: AxisSymetric  = 2
80  integer(kind=kint), parameter :: Shell         = 3
81
82  ! Material constants are saved in an array of size 100 and their physical meaning
83  ! are conrresponds to their position in the array
84  integer(kind=kint), parameter :: M_YOUNGS  = 1
85  integer(kind=kint), parameter :: M_POISSON = 2
86  integer(kind=kint), parameter :: M_DENSITY = 3
87  integer(kind=kint), parameter :: M_THICK   = 4
88
89  ! following plastic constitutive parameter
90  integer(kind=kint), parameter :: M_PLCONST1 = 5
91  integer(kind=kint), parameter :: M_PLCONST2 = 6
92  integer(kind=kint), parameter :: M_PLCONST3 = 7
93  integer(kind=kint), parameter :: M_PLCONST4 = 8
94  integer(kind=kint), parameter :: M_PLCONST5 = 9
95  integer(kind=kint), parameter :: M_KINEHARD = 10
96
97  integer(kind=kint), parameter :: M_EXAPNSION = 20
98
99  integer(kind=kint), parameter :: M_ALPHA_OVER_MU = 21
100
101  integer(kind=kint), parameter :: M_BEAM_RADIUS = 22
102  integer(kind=kint), parameter :: M_BEAM_ANGLE1 = 23
103  integer(kind=kint), parameter :: M_BEAM_ANGLE2 = 24
104  integer(kind=kint), parameter :: M_BEAM_ANGLE3 = 25
105  integer(kind=kint), parameter :: M_BEAM_ANGLE4 = 26
106  integer(kind=kint), parameter :: M_BEAM_ANGLE5 = 27
107  integer(kind=kint), parameter :: M_BEAM_ANGLE6 = 28
108
109  integer(kind=kint), parameter :: M_VISCOCITY = 29
110
111  ! additional plastic constitutive parameter
112  integer(kind=kint), parameter :: M_PLCONST6 = 30
113  integer(kind=kint), parameter :: M_PLCONST7 = 31
114  integer(kind=kint), parameter :: M_PLCONST8 = 32
115  integer(kind=kint), parameter :: M_PLCONST9 = 33
116  integer(kind=kint), parameter :: M_PLCONST10 = 34
117
118  ! Dictionary constants
119  character(len=DICT_KEY_LENGTH) :: MC_ISOELASTIC= 'ISOELASTIC'      ! youngs modulus, poisson's ratio
120  character(len=DICT_KEY_LENGTH) :: MC_ORTHOELASTIC= 'ORTHOELASTIC'  ! ortho elastic modulus
121  character(len=DICT_KEY_LENGTH) :: MC_YIELD = 'YIELD'               ! plastic strain, yield stress
122  character(len=DICT_KEY_LENGTH) :: MC_THEMOEXP = 'THEMOEXP'         ! thermo expansion coefficient
123  character(len=DICT_KEY_LENGTH) :: MC_ORTHOEXP = 'ORTHOEXP'         ! thermo expansion coefficient
124  character(len=DICT_KEY_LENGTH) :: MC_VISCOELASTIC = 'VISCOELASTIC' ! Prony coeff only curr.
125  character(len=DICT_KEY_LENGTH) :: MC_NORTON = 'NORTON'             ! NOrton's creep law
126  character(len=DICT_KEY_LENGTH) :: MC_INCOMP_NEWTONIAN = 'INCOMP_FLUID' ! viscocity
127
128  type tshellmat
129    integer(kind=kint)         :: ortho
130    real(kind=kreal)           :: ee
131    real(kind=kreal)           :: pp
132    real(kind=kreal)           :: ee2
133    real(kind=kreal)           :: g12
134    real(kind=kreal)           :: g23
135    real(kind=kreal)           :: g31
136    real(kind=kreal)           :: angle
137    real(kind=kreal)           :: rho
138    real(kind=kreal)           :: aplha
139    real(kind=kreal)           :: alpha_over_mu
140    real(kind=kreal)           :: weight
141  end type tshellmat
142
143  !> Stucture to management all material relates data
144  type tMaterial
145    integer(kind=kint)         :: nlgeom_flag       !< type of constitutive relation
146    integer(kind=kint)         :: mtype             !< material type
147    integer(kind=kint)         :: nfstatus          !< number of status variables
148    character(len=30)          :: name              !< material name
149    real(kind=kreal)           :: variables(200)    !< material properties
150    type(tshellmat), pointer   :: shell_var(:)      !< material properties for shell
151    integer(kind=kint)         :: totallyr          !< total layer of element
152    integer(kind=kint)         :: cdsys_ID          !< ID of material coordinate system
153    integer(kind=kint)         :: n_table           !< size of table
154    real(kind=kreal), pointer  :: table(:)=>null()  !< material properties in tables
155    type(DICT_STRUCT), pointer :: dict              !< material properties in dictionaried linked list
156  end type tMaterial
157
158  type(tMaterial), allocatable :: materials(:)
159
160contains
161
162  !> Initializer
163  subroutine initMaterial( material )
164    type( tMaterial ), intent(inout) :: material
165    material%mtype = -1                  ! not defined yet
166    material%nfstatus = 0                ! Default: no status
167    material%nlgeom_flag = INFINITESIMAL ! Default: INFINITESIMAL ANALYSIS
168    material%variables =  0.d0           ! not defined yet
169    material%totallyr =  0               ! not defined yet
170
171    call dict_create( material%dict, 'INIT', DICT_NULL )
172  end subroutine
173
174  !> Finalizer
175  subroutine finalizeMaterial( material )
176    type( tMaterial ), intent(inout) :: material
177    if( associated(material%table) ) deallocate( material%table )
178    if( associated(material%dict) ) call dict_destroy( material%dict )
179  end subroutine finalizeMaterial
180
181  !> Initializer
182  subroutine initializeMatls( nm )
183    integer, intent(in) :: nm
184    integer :: i
185    if( allocated(materials) ) deallocate( materials )
186    allocate( materials( nm ) )
187    do i=1,nm
188      call initMaterial( materials(i) )
189    enddo
190  end subroutine
191
192  !> Finalizer
193  subroutine finalizeMatls()
194    integer :: i
195    if( allocated( materials ) ) then
196      do i=1,size(materials)
197        call finalizeMaterial( materials(i) )
198      enddo
199      deallocate( materials )
200    endif
201  end subroutine
202
203  !> Set value of variable(m) of material n to v
204  subroutine modifyMatl( n,m,v)
205    integer, intent(in)  :: n
206    integer, intent(in)  :: m
207    real(kind=kreal), intent(in) :: v
208
209    if( n>size(materials) .OR. m>100 ) return
210    materials(n)%variables(m) = v
211  end subroutine
212
213  !> Print out the material properties
214  subroutine printMaterial( nfile, material )
215    integer, intent(in)           :: nfile
216    type( tMaterial ), intent(in) :: material
217    integer :: i, nt
218    write( nfile, *) "Material type:",material%mtype,material%nlgeom_flag
219    do i=1,100
220      if( material%variables(i) /= 0.d0 ) write( nfile, *) i,material%variables(i)
221    enddo
222    if( associated( material%table ) ) then
223      nt = size(material%table)
224      write( nfile,* ) "--table--"
225      do i=1,nt
226        write(nfile,*) i,material%table(i)
227      enddo
228    endif
229    call print_TableData( material%dict, nfile )
230  end subroutine
231
232  !> Fetch material type
233  integer function fetchDigit( npos, cnum )
234    integer, intent(in) :: npos
235    integer, intent(in) :: cnum
236    integer :: i, idum,cdum,dd
237    fetchDigit = -1
238    cdum = cnum
239    if( npos<=0 .or. npos>6) return
240    if( cnum<100000 .or. cnum>999999 ) return
241    dd = 100000
242    do i=1,npos-1
243      idum = cdum/dd
244      cdum = cdum-idum*dd
245      dd = dd/10
246    enddo
247    fetchDigit = cdum/10**(6-npos)
248  end function
249
250  !> Modify material type
251  subroutine setDigit( npos, ival, mtype )
252    integer, intent(in)              :: npos
253    integer, intent(in)              :: ival
254    integer, intent(inout)           :: mtype
255    integer :: i, idum,cdum, cdum1, dd
256    cdum = mtype
257    if( npos<=0 .or. npos>6 ) return
258    if( ival<0 .or. ival>9 ) return
259    dd =100000
260    cdum1 = 0
261    do i=1,npos-1
262      idum = cdum/dd
263      cdum1 = cdum1+ idum*dd
264      cdum = cdum-idum*dd
265      dd=dd/10
266    enddo
267    cdum1 = cdum1 + ival*dd
268    idum = cdum/dd
269    cdum = cdum-idum*dd
270    dd=dd/10
271    do i=npos+1,6
272      idum = cdum/dd
273      cdum1 = cdum1+ idum*dd
274      cdum = cdum-idum*dd
275      dd=dd/10
276    enddo
277    mtype = cdum1
278  end subroutine
279
280  !> Get elastic type
281  integer function getElasticType( mtype )
282    integer, intent(in) :: mtype
283    integer :: itype
284    getElasticType = -1
285    itype = fetchDigit( 1, mtype )
286    if( itype/=1 ) return  ! not defomration problem
287    itype = fetchDigit( 2, mtype )
288    if( itype/=1 .and. itype/=2 ) return  ! not defomration problem
289    getElasticType = fetchDigit( 3, mtype )
290  end function
291
292  !> Get type of yield function
293  integer function getYieldFunction( mtype )
294    integer, intent(in) :: mtype
295    integer :: itype
296    getYieldFunction = -1
297    itype = fetchDigit( 1, mtype )
298    if( itype/=1 ) return  ! not defomration problem
299    itype = fetchDigit( 2, mtype )
300    if( itype/=2 ) return  ! not elstoplastic problem
301    getYieldFunction = fetchDigit( 4, mtype )
302  end function
303
304  !> Get type of hardening
305  integer function getHardenType( mtype )
306    integer, intent(in) :: mtype
307    integer :: itype
308    getHardenType = -1
309    itype = fetchDigit( 1, mtype )
310    if( itype/=1 ) return  ! not defomration problem
311    itype = fetchDigit( 2, mtype )
312    if( itype/=2 ) return  ! not elstoplastic problem
313    getHardenType = fetchDigit( 5, mtype )
314  end function
315
316  !> If it is a kinematic hardening material?
317  logical function isKinematicHarden( mtype )
318    integer, intent(in) :: mtype
319    integer :: itype
320    isKinematicHarden = .false.
321    itype = fetchDigit( 5, mtype )
322    if( itype==4 .or. itype==5 ) isKinematicHarden = .true.
323  end function
324
325  !> If it is an elastic material?
326  logical function isElastic( mtype )
327    integer, intent(in) :: mtype
328    integer :: itype
329    isElastic = .false.
330    itype = fetchDigit( 2, mtype )
331    if( itype==1 ) isElastic = .true.
332  end function
333
334  !> If it is an elastoplastic material?
335  logical function isElastoplastic( mtype )
336    integer, intent(in) :: mtype
337    integer :: itype
338    isElastoplastic = .false.
339    itype = fetchDigit( 2, mtype )
340    if( itype==2 ) isElastoplastic = .true.
341  end function
342
343  !> If it is a hyperelastic material?
344  logical function isHyperelastic( mtype )
345    integer, intent(in) :: mtype
346    integer :: itype
347    isHyperelastic = .false.
348    itype = fetchDigit( 2, mtype )
349    if( itype==3 ) isHyperelastic = .true.
350  end function
351
352  !> If it is an viscoelastic material?
353  logical function isViscoelastic( mtype )
354    integer, intent(in) :: mtype
355    integer :: itype
356    isViscoelastic = .false.
357    itype = fetchDigit( 2, mtype )
358    if( itype==4 ) isViscoelastic = .true.
359  end function
360
361  !> Set material type of elastoplastic to elastic
362  subroutine ep2e( mtype )
363    integer, intent(inout) :: mtype
364    if( .not. isElastoplastic( mtype ) ) return
365    call setDigit( 2, 1, mtype )
366  end subroutine
367
368end module
369
370
371
372