Polynomial Class

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.)

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)

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) &
+  isScalar = (.not. isKeyValue(line) .and. &
-                                                            .and. .not.isFlow(line))
+              .not. isKey(line) .and. &
+              .not. isListItem(line) .and. &
+              .not. isFlow(line))
 
 end function isScalar
 

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"

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'
 
 

src/phase.f90

@@ -8,6 +8,7 @@ module phase
   use constants
   use math
   use rotations
+  use polynomials
   use IO
   use config
   use material

src/phase_mechanical_elastic.f90

@@ -1,15 +1,13 @@
 submodule(phase:mechanical) elastic
 
   type :: tParameters
-    real(pReal),dimension(3) :: &
+    type(tPolynomial) :: &
-      C_11 = 0.0_pReal, &
+      C_11, &
-      C_12 = 0.0_pReal, &
+      C_12, &
-      C_13 = 0.0_pReal, &
+      C_13, &
-      C_33 = 0.0_pReal, &
+      C_33, &
-      C_44 = 0.0_pReal, &
+      C_44, &
-      C_66 = 0.0_pReal
+      C_66
-    real(pReal) :: &
-      T_ref
   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 = polynomial(elastic%asDict(),'C_11','T')
-
+      prm%C_12 = polynomial(elastic%asDict(),'C_12','T')
-      prm%C_11(1) = elastic%get_asFloat('C_11')
+      prm%C_44 = polynomial(elastic%asDict(),'C_44','T')
-      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)
  
       if (any(phase_lattice(ph) == ['hP','tI'])) then
-        prm%C_13(1) = elastic%get_asFloat('C_13')
+        prm%C_13 = polynomial(elastic%asDict(),'C_13','T')
-        prm%C_13(2) = elastic%get_asFloat('C_13,T',  defaultVal=0.0_pReal)
+        prm%C_33 = polynomial(elastic%asDict(),'C_33','T')
-        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)
       end if
 
-      if (phase_lattice(ph) == 'tI') then
+      if (phase_lattice(ph) == 'tI') &
-        prm%C_66(1) = elastic%get_asFloat('C_66')
+        prm%C_66 = polynomial(elastic%asDict(),'C_66','T')
-        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
 
     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) &
+    C66(1,1) = prm%C_11%at(T)
-             + prm%C_11(2)*(T - prm%T_ref) &
+    C66(1,2) = prm%C_12%at(T)
-             + prm%C_11(3)*(T - prm%T_ref)**2
+    C66(4,4) = prm%C_44%at(T)
-
-    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
-
 
     if (any(phase_lattice(ph) == ['hP','tI'])) then
-      C66(1,3) = prm%C_13(1) &
+      C66(1,3) = prm%C_13%at(T)
-               + prm%C_13(2)*(T - prm%T_ref) &
+      C66(3,3) = prm%C_33%at(T)
-               + 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
-
     end if
 
-    if (phase_lattice(ph) == 'tI') then
+    if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66%at(T)
-      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
 
     C66 = lattice_symmetrize_C66(C66,phase_lattice(ph))
 

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