Merge branch 'polynomial' into 'development'
Polynomial Class See merge request damask/DAMASK!512
This commit is contained in:
Franz Roters
commit
30bf386418
|
|
@ -14,6 +14,7 @@ module CPFEM
|
|||
use config
|
||||
use math
|
||||
use rotations
|
||||
use polynomials
|
||||
use lattice
|
||||
use material
|
||||
use phase
|
||||
|
|
@ -83,6 +84,7 @@ subroutine CPFEM_initAll
|
|||
call config_init
|
||||
call math_init
|
||||
call rotations_init
|
||||
call polynomials_init
|
||||
call lattice_init
|
||||
call discretization_marc_init
|
||||
call material_init(.false.)
|
||||
|
|
|
|||
|
|
@ -16,6 +16,7 @@ module CPFEM2
|
|||
use config
|
||||
use math
|
||||
use rotations
|
||||
use polynomials
|
||||
use lattice
|
||||
use material
|
||||
use phase
|
||||
|
|
@ -57,6 +58,7 @@ subroutine CPFEM_initAll
|
|||
call config_init
|
||||
call math_init
|
||||
call rotations_init
|
||||
call polynomials_init
|
||||
call lattice_init
|
||||
#if defined(MESH)
|
||||
call discretization_mesh_init(restart=interface_restartInc>0)
|
||||
|
|
|
|||
|
|
@ -191,8 +191,10 @@ logical function isScalar(line)
|
|||
|
||||
character(len=*), intent(in) :: line
|
||||
|
||||
isScalar = (.not.isKeyValue(line) .and. .not.isKey(line) .and. .not.isListItem(line) &
|
||||
.and. .not.isFlow(line))
|
||||
isScalar = (.not. isKeyValue(line) .and. &
|
||||
.not. isKey(line) .and. &
|
||||
.not. isListItem(line) .and. &
|
||||
.not. isFlow(line))
|
||||
|
||||
end function isScalar
|
||||
|
||||
|
|
|
|||
|
|
@ -14,6 +14,7 @@
|
|||
#include "LAPACK_interface.f90"
|
||||
#include "math.f90"
|
||||
#include "rotations.f90"
|
||||
#include "polynomials.f90"
|
||||
#include "lattice.f90"
|
||||
#include "element.f90"
|
||||
#include "geometry_plastic_nonlocal.f90"
|
||||
|
|
|
|||
|
|
@ -66,7 +66,7 @@ subroutine material_init(restart)
|
|||
print'(/,1x,a)', '<<<+- material init -+>>>'; flush(IO_STDOUT)
|
||||
|
||||
|
||||
call parse
|
||||
call parse()
|
||||
print'(/,1x,a)', 'parsed material.yaml'
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -8,6 +8,7 @@ module phase
|
|||
use constants
|
||||
use math
|
||||
use rotations
|
||||
use polynomials
|
||||
use IO
|
||||
use config
|
||||
use material
|
||||
|
|
|
|||
|
|
@ -1,15 +1,13 @@
|
|||
submodule(phase:mechanical) elastic
|
||||
|
||||
type :: tParameters
|
||||
real(pReal),dimension(3) :: &
|
||||
C_11 = 0.0_pReal, &
|
||||
C_12 = 0.0_pReal, &
|
||||
C_13 = 0.0_pReal, &
|
||||
C_33 = 0.0_pReal, &
|
||||
C_44 = 0.0_pReal, &
|
||||
C_66 = 0.0_pReal
|
||||
real(pReal) :: &
|
||||
T_ref
|
||||
type(tPolynomial) :: &
|
||||
C_11, &
|
||||
C_12, &
|
||||
C_13, &
|
||||
C_33, &
|
||||
C_44, &
|
||||
C_66
|
||||
end type tParameters
|
||||
|
||||
type(tParameters), allocatable, dimension(:) :: param
|
||||
|
|
@ -47,35 +45,17 @@ module subroutine elastic_init(phases)
|
|||
|
||||
associate(prm => param(ph))
|
||||
|
||||
prm%T_ref = elastic%get_asFloat('T_ref', defaultVal=T_ROOM)
|
||||
|
||||
prm%C_11(1) = elastic%get_asFloat('C_11')
|
||||
prm%C_11(2) = elastic%get_asFloat('C_11,T', defaultVal=0.0_pReal)
|
||||
prm%C_11(3) = elastic%get_asFloat('C_11,T^2',defaultVal=0.0_pReal)
|
||||
|
||||
prm%C_12(1) = elastic%get_asFloat('C_12')
|
||||
prm%C_12(2) = elastic%get_asFloat('C_12,T', defaultVal=0.0_pReal)
|
||||
prm%C_12(3) = elastic%get_asFloat('C_12,T^2',defaultVal=0.0_pReal)
|
||||
|
||||
prm%C_44(1) = elastic%get_asFloat('C_44')
|
||||
prm%C_44(2) = elastic%get_asFloat('C_44,T', defaultVal=0.0_pReal)
|
||||
prm%C_44(3) = elastic%get_asFloat('C_44,T^2',defaultVal=0.0_pReal)
|
||||
prm%C_11 = polynomial(elastic%asDict(),'C_11','T')
|
||||
prm%C_12 = polynomial(elastic%asDict(),'C_12','T')
|
||||
prm%C_44 = polynomial(elastic%asDict(),'C_44','T')
|
||||
|
||||
if (any(phase_lattice(ph) == ['hP','tI'])) then
|
||||
prm%C_13(1) = elastic%get_asFloat('C_13')
|
||||
prm%C_13(2) = elastic%get_asFloat('C_13,T', defaultVal=0.0_pReal)
|
||||
prm%C_13(3) = elastic%get_asFloat('C_13,T^2',defaultVal=0.0_pReal)
|
||||
|
||||
prm%C_33(1) = elastic%get_asFloat('C_33')
|
||||
prm%C_33(2) = elastic%get_asFloat('C_33,T', defaultVal=0.0_pReal)
|
||||
prm%C_33(3) = elastic%get_asFloat('C_33,T^2',defaultVal=0.0_pReal)
|
||||
prm%C_13 = polynomial(elastic%asDict(),'C_13','T')
|
||||
prm%C_33 = polynomial(elastic%asDict(),'C_33','T')
|
||||
end if
|
||||
|
||||
if (phase_lattice(ph) == 'tI') then
|
||||
prm%C_66(1) = elastic%get_asFloat('C_66')
|
||||
prm%C_66(2) = elastic%get_asFloat('C_66,T', defaultVal=0.0_pReal)
|
||||
prm%C_66(3) = elastic%get_asFloat('C_66,T^2',defaultVal=0.0_pReal)
|
||||
end if
|
||||
if (phase_lattice(ph) == 'tI') &
|
||||
prm%C_66 = polynomial(elastic%asDict(),'C_66','T')
|
||||
|
||||
end associate
|
||||
end do
|
||||
|
|
@ -97,38 +77,20 @@ pure module function elastic_C66(ph,en) result(C66)
|
|||
|
||||
|
||||
associate(prm => param(ph))
|
||||
|
||||
C66 = 0.0_pReal
|
||||
T = thermal_T(ph,en)
|
||||
|
||||
C66(1,1) = prm%C_11(1) &
|
||||
+ prm%C_11(2)*(T - prm%T_ref) &
|
||||
+ prm%C_11(3)*(T - prm%T_ref)**2
|
||||
|
||||
C66(1,2) = prm%C_12(1) &
|
||||
+ prm%C_12(2)*(T - prm%T_ref) &
|
||||
+ prm%C_12(3)*(T - prm%T_ref)**2
|
||||
|
||||
C66(4,4) = prm%C_44(1) &
|
||||
+ prm%C_44(2)*(T - prm%T_ref) &
|
||||
+ prm%C_44(3)*(T - prm%T_ref)**2
|
||||
|
||||
C66(1,1) = prm%C_11%at(T)
|
||||
C66(1,2) = prm%C_12%at(T)
|
||||
C66(4,4) = prm%C_44%at(T)
|
||||
|
||||
if (any(phase_lattice(ph) == ['hP','tI'])) then
|
||||
C66(1,3) = prm%C_13(1) &
|
||||
+ prm%C_13(2)*(T - prm%T_ref) &
|
||||
+ prm%C_13(3)*(T - prm%T_ref)**2
|
||||
|
||||
C66(3,3) = prm%C_33(1) &
|
||||
+ prm%C_33(2)*(T - prm%T_ref) &
|
||||
+ prm%C_33(3)*(T - prm%T_ref)**2
|
||||
|
||||
C66(1,3) = prm%C_13%at(T)
|
||||
C66(3,3) = prm%C_33%at(T)
|
||||
end if
|
||||
|
||||
if (phase_lattice(ph) == 'tI') then
|
||||
C66(6,6) = prm%C_66(1) &
|
||||
+ prm%C_66(2)*(T - prm%T_ref) &
|
||||
+ prm%C_66(3)*(T - prm%T_ref)**2
|
||||
end if
|
||||
if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66%at(T)
|
||||
|
||||
C66 = lattice_symmetrize_C66(C66,phase_lattice(ph))
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,179 @@
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @author Martin Diehl, KU Leuven
|
||||
!> @brief Polynomial representation for variable data
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module polynomials
|
||||
use prec
|
||||
use IO
|
||||
use YAML_parse
|
||||
use YAML_types
|
||||
|
||||
implicit none
|
||||
private
|
||||
|
||||
type, public :: tPolynomial
|
||||
real(pReal), dimension(:), allocatable :: coef
|
||||
real(pReal) :: x_ref
|
||||
contains
|
||||
procedure, public :: at => eval
|
||||
procedure, public :: der1_at => eval_der1
|
||||
end type tPolynomial
|
||||
|
||||
interface polynomial
|
||||
module procedure polynomial_from_dict
|
||||
module procedure polynomial_from_coef
|
||||
end interface polynomial
|
||||
|
||||
public :: &
|
||||
polynomial, &
|
||||
polynomials_init
|
||||
|
||||
contains
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Run self-test.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine polynomials_init()
|
||||
|
||||
print'(/,1x,a)', '<<<+- polynomials init -+>>>'; flush(IO_STDOUT)
|
||||
|
||||
call selfTest()
|
||||
|
||||
end subroutine polynomials_init
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Initialize a Polynomial from Coefficients.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function polynomial_from_coef(coef,x_ref) result(p)
|
||||
|
||||
real(pReal), dimension(:), intent(in) :: coef
|
||||
real(pReal), intent(in) :: x_ref
|
||||
type(tPolynomial) :: p
|
||||
|
||||
|
||||
allocate(p%coef(0:size(coef)-1),source=coef) ! should be zero based
|
||||
p%x_ref = x_ref
|
||||
|
||||
end function polynomial_from_coef
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Initialize a Polynomial from a Dictionary with Coefficients.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function polynomial_from_dict(dict,y,x) result(p)
|
||||
|
||||
type(tDict), intent(in) :: dict
|
||||
character(len=*), intent(in) :: y, x
|
||||
type(tPolynomial) :: p
|
||||
|
||||
real(pReal), dimension(:), allocatable :: coef
|
||||
real(pReal) :: x_ref
|
||||
|
||||
|
||||
allocate(coef(1),source=dict%get_asFloat(y))
|
||||
|
||||
if (dict%contains(y//','//x)) then
|
||||
x_ref = dict%get_asFloat(x//'_ref')
|
||||
coef = [coef,dict%get_asFloat(y//','//x)]
|
||||
if (dict%contains(y//','//x//'^2')) then
|
||||
coef = [coef,dict%get_asFloat(y//','//x//'^2')]
|
||||
end if
|
||||
else
|
||||
x_ref = huge(0.0_pReal) ! Simplify debugging
|
||||
end if
|
||||
|
||||
p = Polynomial(coef,x_ref)
|
||||
|
||||
end function polynomial_from_dict
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Evaluate a Polynomial.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function eval(self,x) result(y)
|
||||
|
||||
class(tPolynomial), intent(in) :: self
|
||||
real(pReal), intent(in) :: x
|
||||
real(pReal) :: y
|
||||
|
||||
integer :: i
|
||||
|
||||
|
||||
y = self%coef(0)
|
||||
do i = 1, ubound(self%coef,1)
|
||||
y = y + self%coef(i) * (x-self%x_ref)**i
|
||||
enddo
|
||||
|
||||
end function eval
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Evaluate a first derivative of Polynomial.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function eval_der1(self,x) result(y)
|
||||
|
||||
class(tPolynomial), intent(in) :: self
|
||||
real(pReal), intent(in) :: x
|
||||
real(pReal) :: y
|
||||
|
||||
integer :: i
|
||||
|
||||
|
||||
y = 0.0_pReal
|
||||
do i = 1, ubound(self%coef,1)
|
||||
y = y + real(i,pReal)*self%coef(i) * (x-self%x_ref)**(i-1)
|
||||
enddo
|
||||
|
||||
end function eval_der1
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Check correctness of polynomical functionality.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine selfTest
|
||||
|
||||
type(tPolynomial) :: p1, p2
|
||||
real(pReal), dimension(3) :: coef
|
||||
real(pReal) :: x_ref, x
|
||||
class(tNode), pointer :: dict
|
||||
character(len=pStringLen), dimension(3) :: coef_s
|
||||
character(len=pStringLen) :: x_ref_s, x_s, YAML_s
|
||||
|
||||
call random_number(coef)
|
||||
call random_number(x_ref)
|
||||
call random_number(x)
|
||||
|
||||
coef = coef*10_pReal -0.5_pReal
|
||||
x_ref = x_ref*10_pReal -0.5_pReal
|
||||
x = x*10_pReal -0.5_pReal
|
||||
|
||||
p1 = polynomial(coef,x_ref)
|
||||
if (dNeq(p1%at(x_ref),coef(1))) error stop 'polynomial: @ref'
|
||||
|
||||
write(coef_s(1),*) coef(1)
|
||||
write(coef_s(2),*) coef(2)
|
||||
write(coef_s(3),*) coef(3)
|
||||
write(x_ref_s,*) x_ref
|
||||
write(x_s,*) x
|
||||
YAML_s = 'C: '//trim(adjustl(coef_s(1)))//IO_EOL//&
|
||||
'C,T: '//trim(adjustl(coef_s(2)))//IO_EOL//&
|
||||
'C,T^2: '//trim(adjustl(coef_s(3)))//IO_EOL//&
|
||||
'T_ref: '//trim(adjustl(x_ref_s))//IO_EOL
|
||||
Dict => YAML_parse_str(trim(YAML_s))
|
||||
p2 = polynomial(dict%asDict(),'C','T')
|
||||
if (dNeq(p1%at(x),p2%at(x),1.0e-12_pReal)) error stop 'polynomials: init'
|
||||
|
||||
p1 = polynomial(coef*[0.0_pReal,1.0_pReal,0.0_pReal],x_ref)
|
||||
if (dNeq(p1%at(x_ref+x),-p1%at(x_ref-x),1.0e-10_pReal)) error stop 'polynomials: eval(odd)'
|
||||
if (dNeq(p1%der1_at(x),p1%der1_at(5.0_pReal*x),1.0e-10_pReal)) error stop 'polynomials: eval_der(odd)'
|
||||
|
||||
p1 = polynomial(coef*[0.0_pReal,0.0_pReal,1.0_pReal],x_ref)
|
||||
if (dNeq(p1%at(x_ref+x),p1%at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval(even)'
|
||||
if (dNeq(p1%der1_at(x_ref+x),-p1%der1_at(x_ref-x),1e-10_pReal)) error stop 'polynomials: eval_der(even)'
|
||||
|
||||
|
||||
end subroutine selfTest
|
||||
|
||||
end module polynomials
|
||||
Loading…
Reference in New Issue