From 39aa24369547616fddca666b56fa9248ffe411e4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 31 Jan 2022 14:05:15 +0000 Subject: [PATCH] Polynomial Class --- src/CPFEM.f90 | 2 + src/CPFEM2.f90 | 2 + src/YAML_parse.f90 | 6 +- src/commercialFEM_fileList.f90 | 1 + src/material.f90 | 2 +- src/phase.f90 | 1 + src/phase_mechanical_elastic.f90 | 80 ++++---------- src/polynomials.f90 | 179 +++++++++++++++++++++++++++++++ 8 files changed, 211 insertions(+), 62 deletions(-) create mode 100644 src/polynomials.f90 diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 8ee26ffe7..0f7590782 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -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.) diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 647befd30..ed8fb611b 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -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) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index 9ebde963b..8a3264cff 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -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 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index e67149dea..8dbbea706 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -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" diff --git a/src/material.f90 b/src/material.f90 index 1e3a4b4ec..31eeef5e1 100644 --- a/src/material.f90 +++ b/src/material.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' diff --git a/src/phase.f90 b/src/phase.f90 index 6035b4491..df7557a0b 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -8,6 +8,7 @@ module phase use constants use math use rotations + use polynomials use IO use config use material diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 8972367d5..c57849b1f 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -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)) diff --git a/src/polynomials.f90 b/src/polynomials.f90 new file mode 100644 index 000000000..49a46e2e4 --- /dev/null +++ b/src/polynomials.f90 @@ -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