From 025cbddd00c8cd47886a6329576db5c83740265e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 13:08:32 +0200 Subject: [PATCH 01/45] generic function and interfaces for slipSlip, twinTwin, transTrans --- src/lattice.f90 | 153 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 153 insertions(+) diff --git a/src/lattice.f90 b/src/lattice.f90 index ca1cd597a..d175e2b2d 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2117,4 +2117,157 @@ pure function lattice_qDisorientation(Q1, Q2, struct) end function lattice_qDisorientation + +!-------------------------------------------------------------------------------------------------- +!> @brief Populates reduced slip-slip interaction matrix +!> ToDo: prefix "2" needed as long as deprecated array lattice_interactionSlipSlip exists +!-------------------------------------------------------------------------------------------------- +function lattice_interactionSlipSlip2(Nslip,interactionValues,structure) + use IO, only: & + IO_error + + implicit none + integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family + real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values slip-slip + character(len=*), intent(in) :: structure !< lattice structure + real(pReal), dimension(sum(Nslip),sum(Nslip)) :: lattice_interactionSlipSlip2 + + integer(pInt), dimension(:), allocatable :: NslipMax + integer(pInt), dimension(:,:), allocatable :: interactionSlipSlip + + select case(structure) + case('fcc') + interactionSlipSlip = lattice_fcc_interactionSlipSlip + NslipMax = lattice_fcc_Nslip + case('bcc') + interactionSlipSlip = lattice_bcc_interactionSlipSlip + NslipMax = lattice_bcc_Nslip + case('hex','hexagonal') !ToDo: "No alias policy": long or short? + interactionSlipSlip = lattice_hex_interactionTwinTwin + NslipMax = lattice_hex_Nslip + case('bct') + interactionSlipSlip = lattice_bct_interactionSlipSlip + NslipMax = lattice_bct_Nslip + case default + write(6,*) 'mist' + end select + + !if (size(Ntwin) > count(Ntwin > 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + !if (any(NtwinMax(1:size(Ntwin)) -Ntwin < 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + + lattice_interactionSlipSlip2 = & + buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionSlipSlip) + +end function lattice_interactionSlipSlip2 + + +!-------------------------------------------------------------------------------------------------- +!> @brief Populates reduced twin-twin interaction matrix +!> ToDo: prefix "2" needed as long as deprecated array lattice_interactionTwinTwin exists +!-------------------------------------------------------------------------------------------------- +function lattice_interactionTwinTwin2(Ntwin,interactionValues,structure) + use IO, only: & + IO_error + + implicit none + integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family + real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin + character(len=*), intent(in) :: structure !< lattice structure + real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: lattice_interactionTwinTwin2 + + integer(pInt), dimension(:), allocatable :: NtwinMax + integer(pInt), dimension(:,:), allocatable :: interactionTwinTwin + + select case(structure) + case('fcc') + interactionTwinTwin = lattice_fcc_interactionTwinTwin + NtwinMax = lattice_fcc_Ntwin + case('bcc') + interactionTwinTwin = lattice_bcc_interactionTwinTwin + NtwinMax = lattice_bcc_Ntwin + case('hex','hexagonal') !ToDo: "No alias policy": long or short? + interactionTwinTwin = lattice_hex_interactionTwinTwin + NtwinMax = lattice_hex_Ntwin + case default + write(6,*) 'mist' + end select + + !if (size(Ntwin) > count(Ntwin > 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + !if (any(NtwinMax(1:size(Ntwin)) -Ntwin < 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + + lattice_interactionTwinTwin2 = & + buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTwinTwin) + +end function lattice_interactionTwinTwin2 + + +!-------------------------------------------------------------------------------------------------- +!> @brief Populates reduced trans-trans interaction matrix +!> ToDo: prefix "2" needed as long as deprecated array lattice_interactionTransTrans exists +!-------------------------------------------------------------------------------------------------- +function lattice_interactionTransTrans2(Ntrans,interactionValues,structure,targetStructure) + use IO, only: & + IO_error + + implicit none + integer(pInt), dimension(:), intent(in) :: Ntrans !< number of active twin systems per family + real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin + character(len=*), intent(in) :: & + structure, & !< lattice structure of parent crystal + targetStructure !< lattice structure of transformed crystal + real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: lattice_interactionTransTrans2 + + integer(pInt), dimension(:), allocatable :: NtransMax + integer(pInt), dimension(:,:), allocatable :: interactionTransTrans + + if (trim(structure) == 'fcc' .and. trim(targetStructure) == 'hex') then + interactionTransTrans = lattice_fccToHex_interactionTransTrans + NtransMax = lattice_fcc_Ntrans + else + write(6,*) 'mist' + end if + + !if (size(Ntwin) > count(Ntwin > 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + !if (any(NtwinMax(1:size(Ntwin)) -Ntwin < 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + + lattice_interactionTransTrans2 = & + buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTransTrans) + +end function lattice_interactionTransTrans2 + + +!-------------------------------------------------------------------------------------------------- +!> @brief Populates reduced interaction matrix +!-------------------------------------------------------------------------------------------------- +pure function buildInteraction(activeA,activeB,maxA,maxB,values,matrix) + use IO, only: & + IO_error + + implicit none + integer(pInt), dimension(:), intent(in) :: & + activeA, & + activeB, & + maxA, & + maxB + real(pReal), dimension(:), intent(in) :: values !< interaction values + integer(pInt), dimension(:,:), intent(in) :: matrix + real(pReal), dimension(sum(activeA),sum(activeB)) :: buildInteraction + + integer(pInt) :: & + index_myFamily, index_otherFamily, & + mf, ms, of, os + + myFamilies: do mf = 1_pInt,size(activeA,1) + index_myFamily = sum(activeA(1:mf-1_pInt)) + mySystems: do ms = 1_pInt,activeA(mf) + otherFamilies: do of = 1_pInt,size(activeB,1) + index_otherFamily = sum(activeB(1:of-1_pInt)) + otherSystems: do os = 1_pInt,activeB(of) + buildInteraction(index_myFamily+ms,index_otherFamily+os) = & + values(matrix(sum(maxA(1:mf-1))+ms, sum(maxB(1:of-1))+os)) + enddo otherSystems; enddo otherFamilies; + enddo mySystems;enddo myFamilies + +end function buildInteraction + end module lattice From 1c76e156d976369014a1b7b87b240205f6b6b9db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 17:19:48 +0200 Subject: [PATCH 02/45] new symmetry handling will be handled in a different branch --- src/lattice.f90 | 66 ------------------------------------------------- 1 file changed, 66 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index d175e2b2d..cf4749520 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1183,72 +1183,6 @@ real(pReal), dimension(4,36), parameter, private :: & 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal & ],[4,36]) !< Symmetry operations as quaternions 24 for cubic, 12 for hexagonal = 36 - ! use this later on to substitute the matrix above - ! if self.lattice == 'cubic': - ! symQuats = [ - ! [ 1.0,0.0,0.0,0.0 ], - ! [ 0.0,1.0,0.0,0.0 ], - ! [ 0.0,0.0,1.0,0.0 ], - ! [ 0.0,0.0,0.0,1.0 ], - ! [ 0.0, 0.0, 0.5*math.sqrt(2), 0.5*math.sqrt(2) ], - ! [ 0.0, 0.0, 0.5*math.sqrt(2),-0.5*math.sqrt(2) ], - ! [ 0.0, 0.5*math.sqrt(2), 0.0, 0.5*math.sqrt(2) ], - ! [ 0.0, 0.5*math.sqrt(2), 0.0,-0.5*math.sqrt(2) ], - ! [ 0.0, 0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0 ], - ! [ 0.0,-0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0 ], - ! [ 0.5, 0.5, 0.5, 0.5 ], - ! [-0.5, 0.5, 0.5, 0.5 ], - ! [-0.5, 0.5, 0.5,-0.5 ], - ! [-0.5, 0.5,-0.5, 0.5 ], - ! [-0.5,-0.5, 0.5, 0.5 ], - ! [-0.5,-0.5, 0.5,-0.5 ], - ! [-0.5,-0.5,-0.5, 0.5 ], - ! [-0.5, 0.5,-0.5,-0.5 ], - ! [-0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], - ! [ 0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], - ! [-0.5*math.sqrt(2), 0.0, 0.5*math.sqrt(2), 0.0 ], - ! [-0.5*math.sqrt(2), 0.0,-0.5*math.sqrt(2), 0.0 ], - ! [-0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0, 0.0 ], - ! [-0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0, 0.0 ], - ! ] - ! elif self.lattice == 'hexagonal': - ! symQuats = [ - ! [ 1.0,0.0,0.0,0.0 ], - ! [ 0.0,1.0,0.0,0.0 ], - ! [ 0.0,0.0,1.0,0.0 ], - ! [ 0.0,0.0,0.0,1.0 ], - ! [-0.5*math.sqrt(3), 0.0, 0.0, 0.5 ], - ! [-0.5*math.sqrt(3), 0.0, 0.0,-0.5 ], - ! [ 0.0, 0.5*math.sqrt(3), 0.5, 0.0 ], - ! [ 0.0,-0.5*math.sqrt(3), 0.5, 0.0 ], - ! [ 0.0, 0.5,-0.5*math.sqrt(3), 0.0 ], - ! [ 0.0,-0.5,-0.5*math.sqrt(3), 0.0 ], - ! [ 0.5, 0.0, 0.0, 0.5*math.sqrt(3) ], - ! [-0.5, 0.0, 0.0, 0.5*math.sqrt(3) ], - ! ] - ! elif self.lattice == 'tetragonal': - ! symQuats = [ - ! [ 1.0,0.0,0.0,0.0 ], - ! [ 0.0,1.0,0.0,0.0 ], - ! [ 0.0,0.0,1.0,0.0 ], - ! [ 0.0,0.0,0.0,1.0 ], - ! [ 0.0, 0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0 ], - ! [ 0.0,-0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0 ], - ! [ 0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], - ! [-0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], - ! ] - ! elif self.lattice == 'orthorhombic': - ! symQuats = [ - ! [ 1.0,0.0,0.0,0.0 ], - ! [ 0.0,1.0,0.0,0.0 ], - ! [ 0.0,0.0,1.0,0.0 ], - ! [ 0.0,0.0,0.0,1.0 ], - ! ] - ! else: - ! symQuats = [ - ! [ 1.0,0.0,0.0,0.0 ], - ! ] - public :: & lattice_init, & lattice_qDisorientation, & From d3db0d816ce8260693b518e07d8f45834f550174 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 17:29:20 +0200 Subject: [PATCH 03/45] ordered by dimension and set private where possible --- src/lattice.f90 | 49 +++++++++++++++++++++---------------------------- 1 file changed, 21 insertions(+), 28 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index cf4749520..ed93ee569 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -40,32 +40,20 @@ module lattice real(pReal), allocatable, dimension(:,:,:,:), protected, public :: & lattice_Sslip_v, & !< Mandel notation of lattice_Sslip - lattice_Scleavage_v !< Mandel notation of lattice_Scleavege - - real(pReal), allocatable, dimension(:,:,:), protected, public :: & - lattice_sn, & !< normal direction of slip system - lattice_sd, & !< slip direction of slip system - lattice_st !< sd x sn - -! rotation and Schmid matrices, normal, shear direction and d x n of twin systems - real(pReal), allocatable, dimension(:,:,:,:), protected, public :: & + lattice_Scleavage_v, & !< Mandel notation of lattice_Scleavege + lattice_Qtrans, & !< Total rotation: Q = R*B + lattice_Strans, & !< Eigendeformation tensor for phase transformation lattice_Stwin, & lattice_Qtwin real(pReal), allocatable, dimension(:,:,:), protected, public :: & + lattice_sn, & !< normal direction of slip system + lattice_st, & !< sd x sn + lattice_sd, & !< slip direction of slip system lattice_Stwin_v, & - lattice_tn, & - lattice_td, & - lattice_tt - - real(pReal), allocatable, dimension(:,:,:), protected, public :: & lattice_Strans_v, & !< Eigendeformation tensor in vector form lattice_projectionTrans !< Matrix for projection of slip to fault-band (twin) systems for strain-induced martensite nucleation - real(pReal), allocatable, dimension(:,:,:,:), protected, public :: & - lattice_Qtrans, & !< Total rotation: Q = R*B - lattice_Strans !< Eigendeformation tensor for phase transformation - real(pReal), allocatable, dimension(:,:), protected, public :: & lattice_shearTwin, & !< characteristic twin shear lattice_shearTrans !< characteristic transformation shear @@ -73,6 +61,11 @@ module lattice integer(pInt), allocatable, dimension(:), protected, public :: & lattice_NnonSchmid !< total # of non-Schmid contributions for each structure + real(pReal), allocatable, dimension(:,:,:), private :: & + lattice_tn, & + lattice_td, & + lattice_tt + !-------------------------------------------------------------------------------------------------- ! face centered cubic integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & @@ -1137,9 +1130,7 @@ module lattice integer(kind(LATTICE_undefined_ID)), dimension(:), allocatable, public, protected :: & lattice_structure, trans_lattice_structure - - -integer(pInt), dimension(2), parameter, private :: & + integer(pInt), dimension(2), parameter, private :: & lattice_NsymOperations = [24_pInt,12_pInt] real(pReal), dimension(4,36), parameter, private :: & @@ -1264,16 +1255,10 @@ subroutine lattice_init allocate(lattice_Sslip_v(6,1+2*lattice_maxNnonSchmid,lattice_maxNslip,Nphases),source=0.0_pReal) allocate(lattice_Scleavage(3,3,3,lattice_maxNslip,Nphases),source=0.0_pReal) allocate(lattice_Scleavage_v(6,3,lattice_maxNslip,Nphases),source=0.0_pReal) - allocate(lattice_sd(3,lattice_maxNslip,Nphases),source=0.0_pReal) - allocate(lattice_st(3,lattice_maxNslip,Nphases),source=0.0_pReal) - allocate(lattice_sn(3,lattice_maxNslip,Nphases),source=0.0_pReal) - allocate(lattice_Qtwin(3,3,lattice_maxNtwin,Nphases),source=0.0_pReal) + allocate(lattice_Qtwin(3,3,lattice_maxNtwin,Nphases),source=0.0_pReal) allocate(lattice_Stwin(3,3,lattice_maxNtwin,Nphases),source=0.0_pReal) allocate(lattice_Stwin_v(6,lattice_maxNtwin,Nphases),source=0.0_pReal) - allocate(lattice_td(3,lattice_maxNtwin,Nphases),source=0.0_pReal) - allocate(lattice_tt(3,lattice_maxNtwin,Nphases),source=0.0_pReal) - allocate(lattice_tn(3,lattice_maxNtwin,Nphases),source=0.0_pReal) allocate(lattice_shearTwin(lattice_maxNtwin,Nphases),source=0.0_pReal) allocate(lattice_shearTrans(lattice_maxNtrans,Nphases),source=0.0_pReal) @@ -1301,6 +1286,14 @@ subroutine lattice_init allocate(a_fcc(Nphases),source=0.0_pReal) allocate(a_bcc(Nphases),source=0.0_pReal) + allocate(lattice_td(3,lattice_maxNtwin,Nphases),source=0.0_pReal) + allocate(lattice_tt(3,lattice_maxNtwin,Nphases),source=0.0_pReal) + allocate(lattice_tn(3,lattice_maxNtwin,Nphases),source=0.0_pReal) + allocate(lattice_sd(3,lattice_maxNslip,Nphases),source=0.0_pReal) + allocate(lattice_st(3,lattice_maxNslip,Nphases),source=0.0_pReal) + allocate(lattice_sn(3,lattice_maxNslip,Nphases),source=0.0_pReal) + + do p = 1, size(config_phase) tag = config_phase(p)%getString('lattice_structure') select case(trim(tag)) From 5b3bed3e48d3e420dd5719c77c9eeebdd796a31a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 18:32:55 +0200 Subject: [PATCH 04/45] functionality to calculate SchmidMatrices --- src/lattice.f90 | 105 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 105 insertions(+) diff --git a/src/lattice.f90 b/src/lattice.f90 index ed93ee569..a1fc397ae 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2163,6 +2163,46 @@ function lattice_interactionTransTrans2(Ntrans,interactionValues,structure,targe end function lattice_interactionTransTrans2 +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates Schmid matrix for active slip systems +!-------------------------------------------------------------------------------------------------- +function lattice_SchmidSlip(Nslip,structure,cOverA) + use IO, only: & + IO_error + use math, only: & + math_tensorproduct33 + + implicit none + integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family + character(len=*), intent(in) :: structure !< lattice structure + real(pReal), dimension(3,3,sum(Nslip)) :: lattice_SchmidSlip + real(pReal), intent(in), optional :: & + cOverA + + real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem + integer(pInt) :: i + + select case(structure) + case('fcc') + coordinateSystem = buildCoordinateSystem(Nslip,int(lattice_fcc_systemSlip,pInt),structure) + case('bcc') + coordinateSystem = buildCoordinateSystem(Nslip,int(lattice_bcc_systemSlip,pInt),structure) + case('hex','hexagonal') !ToDo: "No alias policy": long or short? + coordinateSystem = buildCoordinateSystem(Nslip,int(lattice_hex_systemSlip,pInt),structure,cOverA) + case('bct') + coordinateSystem = buildCoordinateSystem(Nslip,int(lattice_bct_systemslip,pInt),structure,cOverA) + case default + write(6,*) 'mist' + end select + + do i = 1, sum(Nslip) + lattice_SchmidSlip(1:3,1:3,i) = & + math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) + enddo + +end function lattice_SchmidSlip + + !-------------------------------------------------------------------------------------------------- !> @brief Populates reduced interaction matrix !-------------------------------------------------------------------------------------------------- @@ -2197,4 +2237,69 @@ pure function buildInteraction(activeA,activeB,maxA,maxB,values,matrix) end function buildInteraction + +!-------------------------------------------------------------------------------------------------- +!> @brief build a local coordinate system in a slip, twin, trans, cleavage system +!-------------------------------------------------------------------------------------------------- +pure function buildCoordinateSystem(active,system,structure,cOverA) + use math, only: & + math_crossproduct + + implicit none + integer(pInt), dimension(:), intent(in) :: & + active + integer(pInt), dimension(:,:), intent(in) :: & + system + character(len=*), intent(in) :: & + structure !< lattice structure + real(pReal), intent(in), optional :: & + cOverA + real(pReal), dimension(3,3,sum(active)) :: & + buildCoordinateSystem + + real(pReal), dimension(3) :: & + direction, normal + integer(pInt) :: ir, ig, mf, ms + + ir = 0_pInt + myFamilies: do mf = 1_pInt,size(active,1) + mySystems: do ms = 1_pInt,active(mf) + ir = ir + 1_pInt + ig = sum(system(:,1:mf-1))+ms + + select case(trim(structure)) + + case ('fcc','bcc') + direction = real(system(1:3,ig),pReal) + normal = real(system(4:6,ig),pReal) + + case ('hex') + ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)]) + direction(1) = real(system(1,ig),pReal)*1.5_pReal + direction(2) = (real(system(1,ig),pReal)+2.0_pReal*real(system(2,ig),pReal))*0.5_pReal*sqrt(3.0_pReal) + direction(3) = real(system(4,ig),pReal)*CoverA + + ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) + normal(1) = real(system(5,ig),pReal) + normal(2) = (real(system(5,ig),pReal)+2.0_pReal*real(system(6,ig),pReal))/ sqrt(3.0_pReal) + normal(3) = real(system(8,ig),pReal)/CoverA + + case ('bct') + direction(1:2) = real(system(1:2,ig),pReal) + direction(3) = real(system(3,ig),pReal)*CoverA + + normal(1:2) = real(system(4:5,ig),pReal) + normal(3) = real(system(6,ig),pReal)/CoverA + + end select + + buildCoordinateSystem(1:3,1,ir) = direction/norm2(direction) ! make unit vector + buildCoordinateSystem(1:3,2,ir) = normal/norm2(normal) ! make unit vector + buildCoordinateSystem(1:3,3,ir) = math_crossproduct(direction,normal) + + enddo mySystems + enddo myFamilies + +end function buildCoordinateSystem + end module lattice From 3d49c70dbc609c7133f21715330481a1bdf5bd97 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 8 Sep 2018 19:32:26 +0200 Subject: [PATCH 05/45] commenting --- src/lattice.f90 | 88 +++++++++++++++++++++++++++++++++++-------------- 1 file changed, 64 insertions(+), 24 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index eff7f48e5..eee35bb83 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -19,6 +19,7 @@ module lattice LATTICE_maxNtransFamily = 1_pInt, & !< max # of transformation system families over lattice structures LATTICE_maxNcleavageFamily = 3_pInt !< max # of transformation system families over lattice structures +! BEGIN DEPRECATED integer(pInt), allocatable, dimension(:,:), protected, public :: & lattice_NslipSystem, & !< total # of slip systems in each family lattice_NtwinSystem, & !< total # of twin systems in each family @@ -65,6 +66,8 @@ module lattice lattice_tn, & lattice_td, & lattice_tt +! END DEPRECATED + !-------------------------------------------------------------------------------------------------- ! face centered cubic @@ -2226,15 +2229,15 @@ function lattice_SchmidSlip(Nslip,structure,cOverA) select case(structure) case('fcc') - coordinateSystem = buildCoordinateSystem(Nslip,int(lattice_fcc_systemSlip,pInt),structure) + coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_FCC_SYSTEMSLIP,pInt),structure) case('bcc') - coordinateSystem = buildCoordinateSystem(Nslip,int(lattice_bcc_systemSlip,pInt),structure) + coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMSLIP,pInt),structure) case('hex','hexagonal') !ToDo: "No alias policy": long or short? - coordinateSystem = buildCoordinateSystem(Nslip,int(lattice_hex_systemSlip,pInt),structure,cOverA) + coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_HEX_SYSTEMSLIP,pInt),'hex',cOverA) case('bct') - coordinateSystem = buildCoordinateSystem(Nslip,int(lattice_bct_systemslip,pInt),structure,cOverA) + coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCT_SYSTEMSLIP,pInt),structure,cOverA) case default - write(6,*) 'mist' + call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidSlip)') end select do i = 1, sum(Nslip) @@ -2245,21 +2248,57 @@ function lattice_SchmidSlip(Nslip,structure,cOverA) end function lattice_SchmidSlip +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates Schmid matrix for active twin systems +!-------------------------------------------------------------------------------------------------- +function lattice_SchmidTwin(Ntwin,structure,cOverA) + use IO, only: & + IO_error + use math, only: & + math_tensorproduct33 + + implicit none + integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family + character(len=*), intent(in) :: structure !< lattice structure + real(pReal), dimension(3,3,sum(Ntwin)) :: lattice_SchmidTwin + real(pReal), intent(in), optional :: & + cOverA + + real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem + integer(pInt) :: i + + select case(structure) + case('fcc') + coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) + case('bcc') + coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) + case('hex','hexagonal') !ToDo: "No alias policy": long or short? + coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) + case default + call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidTwin)') + end select + + do i = 1, sum(Nslip) + lattice_SchmidTwin(1:3,1:3,i) = & + math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) + enddo + +end function lattice_SchmidTwin + + !-------------------------------------------------------------------------------------------------- !> @brief Populates reduced interaction matrix !-------------------------------------------------------------------------------------------------- pure function buildInteraction(activeA,activeB,maxA,maxB,values,matrix) - use IO, only: & - IO_error implicit none integer(pInt), dimension(:), intent(in) :: & - activeA, & - activeB, & - maxA, & - maxB + activeA, & !< number of active systems as specified in material.config + activeB, & !< number of active systems as specified in material.config + maxA, & !< number of maximum available systems + maxB !< number of maximum available systems real(pReal), dimension(:), intent(in) :: values !< interaction values - integer(pInt), dimension(:,:), intent(in) :: matrix + integer(pInt), dimension(:,:), intent(in) :: matrix !< full interaction matrix real(pReal), dimension(sum(activeA),sum(activeB)) :: buildInteraction integer(pInt) :: & @@ -2301,7 +2340,11 @@ pure function buildCoordinateSystem(active,system,structure,cOverA) real(pReal), dimension(3) :: & direction, normal - integer(pInt) :: ir, ig, mf, ms + integer(pInt) :: & + ir, & !< index in reduced matrix + ig, & !< index in full matrix + mf, & !< index of my family + ms !< index of my system in current family ir = 0_pInt myFamilies: do mf = 1_pInt,size(active,1) @@ -2317,21 +2360,18 @@ pure function buildCoordinateSystem(active,system,structure,cOverA) case ('hex') ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)]) - direction(1) = real(system(1,ig),pReal)*1.5_pReal - direction(2) = (real(system(1,ig),pReal)+2.0_pReal*real(system(2,ig),pReal))*0.5_pReal*sqrt(3.0_pReal) - direction(3) = real(system(4,ig),pReal)*CoverA + direction = [ real(system(1,ig),pReal)*1.5_pReal, & + (real(system(1,ig),pReal)+2.0_pReal*real(system(2,ig),pReal))*sqrt(0.75_pReal), & + real(system(4,ig),pReal)*CoverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) - normal(1) = real(system(5,ig),pReal) - normal(2) = (real(system(5,ig),pReal)+2.0_pReal*real(system(6,ig),pReal))/ sqrt(3.0_pReal) - normal(3) = real(system(8,ig),pReal)/CoverA + normal = [ real(system(5,ig),pReal), & + (real(system(5,ig),pReal)+2.0_pReal*real(system(6,ig),pReal))/ sqrt(3.0_pReal), & + real(system(8,ig),pReal)/CoverA ] case ('bct') - direction(1:2) = real(system(1:2,ig),pReal) - direction(3) = real(system(3,ig),pReal)*CoverA - - normal(1:2) = real(system(4:5,ig),pReal) - normal(3) = real(system(6,ig),pReal)/CoverA + direction = [real(system(1:2,ig),pReal),real(system(3,ig),pReal)*CoverA] + normal = [real(system(4:5,ig),pReal),real(system(6,ig),pReal)/CoverA] end select From f1f8922ab81199c0ddba6372aeb0b1cb607cf24c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 12 Sep 2018 14:03:45 +0200 Subject: [PATCH 06/45] more meaningful error messages --- src/IO.f90 | 4 +++ src/lattice.f90 | 73 +++++++++++++++++++++++++++++++++---------------- 2 files changed, 54 insertions(+), 23 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 8e1b9e80f..2daae2fa6 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -1458,6 +1458,10 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'empty list' case (143_pInt) msg = 'no value found for key' + case (144_pInt) + msg = 'negative number of systems' + case (145_pInt) + msg = 'to many systems requested' !-------------------------------------------------------------------------------------------------- ! material error messages and related messages in mesh diff --git a/src/lattice.f90 b/src/lattice.f90 index eee35bb83..11abef2a5 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1957,6 +1957,7 @@ pure function lattice_symmetrizeC66(struct,C66) end function lattice_symmetrizeC66 + !-------------------------------------------------------------------------------------------------- !> @brief Symmetrizes 2nd order tensor according to lattice type !-------------------------------------------------------------------------------------------------- @@ -2090,6 +2091,32 @@ pure function lattice_qDisorientation(Q1, Q2, struct) end function lattice_qDisorientation +!function lattice_nonSchmidMatrix +! coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMSLIP,pInt),structure) +! +! do i = 1_pInt,myNslip ! assign slip system vectors +! direction = coordinateSystem(1:3,1,i) +! normal = coordinateSystem(1:3,1,i) +! ! "np" and "nn" according to Gröger_etal2008, Acta Materialia 56 (2008) 5412–5425, table 1 (corresponds to their "n1" for positive and negative slip direction respectively) +! np = math_mul33x3(math_axisAngleToR(sdU,+60.0_pReal*INRAD), normal) +! ! Schmid matrices with non-Schmid contributions according to Koester_etal2012, Acta Materialia 60 (2012) 3894–3901, eq. (17) ("n1" is replaced by either "np" or "nn" according to either positive or negative slip direction) +! nonSchmidMatrix_pos(1:3,1:3,1,i) = math_tensorproduct33(direction, np) +! nonSchmidMatrix_pos(1:3,1:3,2,i) = math_tensorproduct33(math_crossproduct(normal, direction), normal) +! nonSchmidMatrix_pos(1:3,1:3,3,i) = math_tensorproduct33(math_crossproduct(np, direction), np) +! nonSchmidMatrix_pos(1:3,1:3,4,i) = math_tensorproduct33(normal, normal) +! nonSchmidMatrix_pos(1:3,1:3,5,i) = math_tensorproduct33(math_crossproduct(normal, direction), math_crossproduct(normal, direction)) +! nonSchmidMatrix_pos(1:3,1:3,6,i) = math_tensorproduct33(direction, direction) +! +! nn = math_mul33x3(math_axisAngleToR(-sdU,60.0_pReal*INRAD), normal) +! sns(1:3,1:3,2,1,i) = math_tensorproduct33(-sdU, nn) +! sns(1:3,1:3,2,2,i) = math_tensorproduct33(math_crossproduct(snU, -sdU), snU) +! sns(1:3,1:3,2,3,i) = math_tensorproduct33(math_crossproduct(nn, -sdU), nn) +! sns(1:3,1:3,2,4,i) = math_tensorproduct33(snU, snU) +! sns(1:3,1:3,2,5,i) = math_tensorproduct33(math_crossproduct(snU, -sdU), math_crossproduct(snU, -sdU)) +! sns(1:3,1:3,2,6,i) = math_tensorproduct33(-sdU, -sdU) +! enddo +!end function lattice_nonSchmidMatrix + !-------------------------------------------------------------------------------------------------- !> @brief Populates reduced slip-slip interaction matrix !> ToDo: prefix "2" needed as long as deprecated array lattice_interactionSlipSlip exists @@ -2121,11 +2148,11 @@ function lattice_interactionSlipSlip2(Nslip,interactionValues,structure) interactionSlipSlip = lattice_bct_interactionSlipSlip NslipMax = lattice_bct_Nslip case default - write(6,*) 'mist' + call IO_error(132_pInt,ext_msg=trim(structure)//' (slip slip interaction)') end select - !if (size(Ntwin) > count(Ntwin > 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') - !if (any(NtwinMax(1:size(Ntwin)) -Ntwin < 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + if (any(Nslip(1:size(Nslip)) - Nslip < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) lattice_interactionSlipSlip2 = & buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionSlipSlip) @@ -2161,11 +2188,11 @@ function lattice_interactionTwinTwin2(Ntwin,interactionValues,structure) interactionTwinTwin = lattice_hex_interactionTwinTwin NtwinMax = lattice_hex_Ntwin case default - write(6,*) 'mist' + call IO_error(132_pInt,ext_msg=trim(structure)//' (twin twin interaction)') end select - !if (size(Ntwin) > count(Ntwin > 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') - !if (any(NtwinMax(1:size(Ntwin)) -Ntwin < 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + if (any(Ntwin(1:size(Ntwin)) - Ntwin < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) lattice_interactionTwinTwin2 = & buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTwinTwin) @@ -2196,11 +2223,11 @@ function lattice_interactionTransTrans2(Ntrans,interactionValues,structure,targe interactionTransTrans = lattice_fccToHex_interactionTransTrans NtransMax = lattice_fcc_Ntrans else - write(6,*) 'mist' + call IO_error(132_pInt,ext_msg=trim(structure)//' => '//trim(targetStructure)) end if - !if (size(Ntwin) > count(Ntwin > 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') - !if (any(NtwinMax(1:size(Ntwin)) -Ntwin < 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + if (any(Ntrans(1:size(Ntrans)) - Ntrans < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) lattice_interactionTransTrans2 = & buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTransTrans) @@ -2211,7 +2238,7 @@ end function lattice_interactionTransTrans2 !-------------------------------------------------------------------------------------------------- !> @brief Calculates Schmid matrix for active slip systems !-------------------------------------------------------------------------------------------------- -function lattice_SchmidSlip(Nslip,structure,cOverA) +function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) use IO, only: & IO_error use math, only: & @@ -2220,7 +2247,7 @@ function lattice_SchmidSlip(Nslip,structure,cOverA) implicit none integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(3,3,sum(Nslip)) :: lattice_SchmidSlip + real(pReal), dimension(3,3,sum(Nslip)) :: lattice_SchmidMatrix_slip real(pReal), intent(in), optional :: & cOverA @@ -2237,21 +2264,21 @@ function lattice_SchmidSlip(Nslip,structure,cOverA) case('bct') coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCT_SYSTEMSLIP,pInt),structure,cOverA) case default - call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidSlip)') + call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_slip)') end select do i = 1, sum(Nslip) - lattice_SchmidSlip(1:3,1:3,i) = & + lattice_SchmidMatrix_slip(1:3,1:3,i) = & math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) enddo -end function lattice_SchmidSlip +end function lattice_SchmidMatrix_slip !-------------------------------------------------------------------------------------------------- !> @brief Calculates Schmid matrix for active twin systems !-------------------------------------------------------------------------------------------------- -function lattice_SchmidTwin(Ntwin,structure,cOverA) +function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) use IO, only: & IO_error use math, only: & @@ -2260,7 +2287,7 @@ function lattice_SchmidTwin(Ntwin,structure,cOverA) implicit none integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(3,3,sum(Ntwin)) :: lattice_SchmidTwin + real(pReal), dimension(3,3,sum(Ntwin)) :: lattice_SchmidMatrix_Twin real(pReal), intent(in), optional :: & cOverA @@ -2269,21 +2296,21 @@ function lattice_SchmidTwin(Ntwin,structure,cOverA) select case(structure) case('fcc') - coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) + coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) case('bcc') - coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) + coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) case('hex','hexagonal') !ToDo: "No alias policy": long or short? - coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) + coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) case default - call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidTwin)') + call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_twin)') end select - do i = 1, sum(Nslip) - lattice_SchmidTwin(1:3,1:3,i) = & + do i = 1, sum(Ntwin) + lattice_SchmidMatrix_twin(1:3,1:3,i) = & math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) enddo -end function lattice_SchmidTwin +end function lattice_SchmidMatrix_twin !-------------------------------------------------------------------------------------------------- From 0a621ece8b90f285f1cc3cd2effc2571d15e96d6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 12 Sep 2018 14:29:22 +0200 Subject: [PATCH 07/45] WIP: polishing and implementing remaining functions --- src/lattice.f90 | 29 +++++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 11abef2a5..7fcf735c7 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2090,6 +2090,26 @@ pure function lattice_qDisorientation(Q1, Q2, struct) end function lattice_qDisorientation +!function lattice_C66_twin +! +! select case(structure) +! case('fcc') +! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) +! case('bcc') +! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) +! case('hex','hexagonal') !ToDo: "No alias policy": long or short? +! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) +! case default +! call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_twin)') +! end select +! +! do i = 1, sum(Ntwin) +! R = math_axisAngleToR(coordinateSystem(1:3,2,i), 180.0_pReal * INRAD) ! ToDo: Why always 180 deg? +! math_rotate_forward3333(C,R) +! C_twin66(1:6,1:6,i) = math_Mandel3333to66(C_twin) +! enddo + +end function !function lattice_nonSchmidMatrix ! coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMSLIP,pInt),structure) @@ -2238,7 +2258,7 @@ end function lattice_interactionTransTrans2 !-------------------------------------------------------------------------------------------------- !> @brief Calculates Schmid matrix for active slip systems !-------------------------------------------------------------------------------------------------- -function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) +function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) use IO, only: & IO_error use math, only: & @@ -2247,7 +2267,7 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) implicit none integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(3,3,sum(Nslip)) :: lattice_SchmidMatrix_slip + real(pReal), dimension(3,3,sum(Nslip)) :: SchmidMatrix real(pReal), intent(in), optional :: & cOverA @@ -2268,8 +2288,9 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) end select do i = 1, sum(Nslip) - lattice_SchmidMatrix_slip(1:3,1:3,i) = & - math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) + SchmidMatrix(1:3,1:3,i) = math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) + if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & + call IO_error(0_pInt,i,ext_msg = 'dilatational Schmid matrix for slip') enddo end function lattice_SchmidMatrix_slip From 192bb6453d7b4e30d4d9d7246a12ffb024c9d671 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 12 Sep 2018 15:57:54 +0200 Subject: [PATCH 08/45] cleaning --- src/lattice.f90 | 75 ++++++++++++++++++++++++++----------------------- 1 file changed, 40 insertions(+), 35 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 7fcf735c7..e5bd1e65f 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1135,10 +1135,7 @@ module lattice real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & lattice_C3333, lattice_trans_C3333 real(pReal), dimension(:), allocatable, public, protected :: & - lattice_mu, & - lattice_nu, & - lattice_trans_mu, & - lattice_trans_nu + lattice_mu, lattice_nu real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & ! with higher-order parameters (e.g. temperature-dependent) lattice_thermalExpansion33 real(pReal), dimension(:,:,:), allocatable, public, protected :: & @@ -1292,8 +1289,6 @@ subroutine lattice_init allocate(lattice_mu(Nphases), source=0.0_pReal) allocate(lattice_nu(Nphases), source=0.0_pReal) - allocate(lattice_trans_mu(Nphases), source=0.0_pReal) - allocate(lattice_trans_nu(Nphases), source=0.0_pReal) allocate(lattice_NnonSchmid(Nphases), source=0_pInt) allocate(lattice_Sslip(3,3,1+2*lattice_maxNnonSchmid,lattice_maxNslip,Nphases),source=0.0_pReal) @@ -1528,8 +1523,6 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) select case(trans_lattice_structure(myPhase)) case (LATTICE_bcc_ID) lattice_trans_C66(1:6,1:6,myPhase) = lattice_C66(1:6,1:6,myPhase) - lattice_trans_mu(myPhase) = lattice_mu(myPhase) - lattice_trans_nu(myPhase) = lattice_nu(myPhase) lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = lattice_C3333(1:3,1:3,1:3,1:3,myPhase) lattice_trans_C66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) do i = 1_pInt, 6_pInt @@ -1554,15 +1547,6 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) lattice_trans_C66(1:6,1:6,myPhase) = lattice_symmetrizeC66(trans_lattice_structure(myPhase),& lattice_trans_C66(1:6,1:6,myPhase)) - lattice_trans_mu(myPhase) = 0.2_pReal *( lattice_trans_C66(1,1,myPhase) & - - lattice_trans_C66(1,2,myPhase) & - + 3.0_pReal*lattice_trans_C66(4,4,myPhase)) - lattice_trans_nu(myPhase) = ( lattice_trans_C66(1,1,myPhase) & - + 4.0_pReal*lattice_trans_C66(1,2,myPhase) & - - 2.0_pReal*lattice_trans_C66(4,4,myPhase)) & - /( 4.0_pReal*lattice_trans_C66(1,1,myPhase) & - + 6.0_pReal*lattice_trans_C66(1,2,myPhase) & - + 2.0_pReal*lattice_trans_C66(4,4,myPhase)) lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = math_Voigt66to3333(lattice_trans_C66(1:6,1:6,myPhase)) lattice_trans_C66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) do i = 1_pInt, 6_pInt @@ -2090,24 +2074,42 @@ pure function lattice_qDisorientation(Q1, Q2, struct) end function lattice_qDisorientation -!function lattice_C66_twin -! -! select case(structure) -! case('fcc') -! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) -! case('bcc') -! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) -! case('hex','hexagonal') !ToDo: "No alias policy": long or short? -! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) -! case default -! call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_twin)') -! end select -! -! do i = 1, sum(Ntwin) -! R = math_axisAngleToR(coordinateSystem(1:3,2,i), 180.0_pReal * INRAD) ! ToDo: Why always 180 deg? -! math_rotate_forward3333(C,R) -! C_twin66(1:6,1:6,i) = math_Mandel3333to66(C_twin) -! enddo +function lattice_C66_twin(Ntwin,C66,structure,CoverA) + use IO, only: & + IO_error + use math, only: & + INRAD, & + math_axisAngleToR, & + math_Mandel3333to66, & + math_Mandel66to3333, & + math_rotate_forward3333 + + implicit none + integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family + character(len=*), intent(in) :: structure !< lattice structure + real(pReal), dimension(6,6), intent(in) :: C66 + real(pReal), intent(in) :: cOverA + real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin + + real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem + + real(pReal), dimension(3,3) :: R + integer(pInt) :: i + + select case(structure) + case('fcc') + coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) + case('bcc') + coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) + case('hex','hexagonal') !ToDo: "No alias policy": long or short? + coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) + case default + call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_C66_twin)') + end select + do i = 1, sum(Ntwin) + R = math_axisAngleToR(coordinateSystem(1:3,2,i), 180.0_pReal * INRAD) ! ToDo: Why always 180 deg? + lattice_C66_twin(1:6,1:6,i) = math_Mandel3333to66(math_rotate_forward3333(math_Mandel66to3333(C66),R)) + enddo end function @@ -2259,9 +2261,12 @@ end function lattice_interactionTransTrans2 !> @brief Calculates Schmid matrix for active slip systems !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) + use prec, only: & + tol_math_check use IO, only: & IO_error use math, only: & + math_trace33, & math_tensorproduct33 implicit none From b95174a8b7b231c0abbb97700c7d70002cdbbf35 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 12 Sep 2018 16:01:57 +0200 Subject: [PATCH 09/45] reducing global variables --- src/lattice.f90 | 123 ++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 103 insertions(+), 20 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index e5bd1e65f..28bd49f34 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1130,6 +1130,8 @@ module lattice !maxval(lattice_bct_interactionTwinTwin))) ) !< max # of interaction types (in hardening matrix part) #endif + real(pReal), dimension(:,:,:), allocatable, private :: & + temp66 real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_C66, lattice_trans_C66 real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & @@ -1370,15 +1372,15 @@ subroutine lattice_init lattice_C66(5,5,p) = config_phase(p)%getFloat('c55',defaultVal=0.0_pReal) lattice_C66(6,6,p) = config_phase(p)%getFloat('c66',defaultVal=0.0_pReal) - lattice_trans_C66(1,1,p) = config_phase(p)%getFloat('c11_trans',defaultVal=0.0_pReal) - lattice_trans_C66(1,2,p) = config_phase(p)%getFloat('c12_trans',defaultVal=0.0_pReal) - lattice_trans_C66(1,3,p) = config_phase(p)%getFloat('c13_trans',defaultVal=0.0_pReal) - lattice_trans_C66(2,2,p) = config_phase(p)%getFloat('c22_trans',defaultVal=0.0_pReal) - lattice_trans_C66(2,3,p) = config_phase(p)%getFloat('c23_trans',defaultVal=0.0_pReal) - lattice_trans_C66(3,3,p) = config_phase(p)%getFloat('c33_trans',defaultVal=0.0_pReal) - lattice_trans_C66(4,4,p) = config_phase(p)%getFloat('c44_trans',defaultVal=0.0_pReal) - lattice_trans_C66(5,5,p) = config_phase(p)%getFloat('c55_trans',defaultVal=0.0_pReal) - lattice_trans_C66(6,6,p) = config_phase(p)%getFloat('c66_trans',defaultVal=0.0_pReal) + temp66(1,1,p) = config_phase(p)%getFloat('c11_trans',defaultVal=0.0_pReal) + temp66(1,2,p) = config_phase(p)%getFloat('c12_trans',defaultVal=0.0_pReal) + temp66(1,3,p) = config_phase(p)%getFloat('c13_trans',defaultVal=0.0_pReal) + temp66(2,2,p) = config_phase(p)%getFloat('c22_trans',defaultVal=0.0_pReal) + temp66(2,3,p) = config_phase(p)%getFloat('c23_trans',defaultVal=0.0_pReal) + temp66(3,3,p) = config_phase(p)%getFloat('c33_trans',defaultVal=0.0_pReal) + temp66(4,4,p) = config_phase(p)%getFloat('c44_trans',defaultVal=0.0_pReal) + temp66(5,5,p) = config_phase(p)%getFloat('c55_trans',defaultVal=0.0_pReal) + temp66(6,6,p) = config_phase(p)%getFloat('c66_trans',defaultVal=0.0_pReal) CoverA(p) = config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal) CoverA_trans(p) = config_phase(p)%getFloat('c/a_trans',defaultVal=0.0_pReal) @@ -1522,9 +1524,9 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) case (LATTICE_fcc_ID) select case(trans_lattice_structure(myPhase)) case (LATTICE_bcc_ID) - lattice_trans_C66(1:6,1:6,myPhase) = lattice_C66(1:6,1:6,myPhase) + temp66(1:6,1:6,myPhase) = lattice_C66(1:6,1:6,myPhase) lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = lattice_C3333(1:3,1:3,1:3,1:3,myPhase) - lattice_trans_C66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) + temp66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) do i = 1_pInt, 6_pInt if (abs(lattice_trans_C66(i,i,myPhase))bcc transformation') @@ -1539,18 +1541,18 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) /(3.0_pReal*sqrt(2.0_pReal)) A = c14bar**(2.0_pReal)/c44bar B = c14bar**(2.0_pReal)/(0.5_pReal*(c11bar - c12bar)) - lattice_trans_C66(1,1,myPhase) = c11bar - A - lattice_trans_C66(1,2,myPhase) = c12bar + A - lattice_trans_C66(1,3,myPhase) = c13bar - lattice_trans_C66(3,3,myPhase) = c33bar - lattice_trans_C66(4,4,myPhase) = c44bar - B + temp66(1,1,myPhase) = c11bar - A + temp66(1,2,myPhase) = c12bar + A + temp66(1,3,myPhase) = c13bar + temp66(3,3,myPhase) = c33bar + temp66(4,4,myPhase) = c44bar - B - lattice_trans_C66(1:6,1:6,myPhase) = lattice_symmetrizeC66(trans_lattice_structure(myPhase),& + temp66(1:6,1:6,myPhase) = lattice_symmetrizeC66(trans_lattice_structure(myPhase),& lattice_trans_C66(1:6,1:6,myPhase)) - lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = math_Voigt66to3333(lattice_trans_C66(1:6,1:6,myPhase)) - lattice_trans_C66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) + lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = math_Voigt66to3333(temp66(1:6,1:6,myPhase)) + temp66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) do i = 1_pInt, 6_pInt - if (abs(lattice_trans_C66(i,i,myPhase))hex transformation') enddo end select @@ -2113,6 +2115,86 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) end function + +!function lattice_C66_trans(Ntrans,C66_parent,C66_targetstructure_parent,cOverA_parent,structure_target,cOverA_target) +! use IO, only: & +! IO_error +! use math, only: & +! INRAD, & +! math_axisAngleToR, & +! math_Mandel3333to66, & +! math_Mandel66to3333, & +! math_rotate_forward3333 +! +! implicit none +! integer(pInt), dimension(:), intent(in) :: Ntrans !< number of active twin systems per family +! character(len=*), intent(in) :: & +! structure_target, & !< lattice structure +! structure_parent !< lattice structure +! real(pReal), dimension(6,6), intent(in) :: C66_target +! real(pReal), intent(in) :: cOverA_parent, cOverA_target +! real(pReal), dimension(6,6,sum(Ntarget)) :: lattice_C66_trans +! +! real(pReal), dimension(3,3,sum(Ntarget)) :: coordinateSystem +! +! real(pReal), dimension(3,3) :: R +! integer(pInt) :: i +! +! if (trim(structure) == 'fcc' .and. trim(targetStructure) == 'hex') then +! c11bar = (lattice_C66(1,1,myPhase) + lattice_C66(1,2,myPhase) + 2.0_pReal*lattice_C66(4,4,myPhase))/2.0_pReal +! c12bar = (lattice_C66(1,1,myPhase) + 5.0_pReal*lattice_C66(1,2,myPhase) - 2.0_pReal*lattice_C66(4,4,myPhase))/6.0_pReal +! c33bar = (lattice_C66(1,1,myPhase) + 2.0_pReal*lattice_C66(1,2,myPhase) + 4.0_pReal*lattice_C66(4,4,myPhase))/3.0_pReal +! c13bar = (lattice_C66(1,1,myPhase) + 2.0_pReal*lattice_C66(1,2,myPhase) - 2.0_pReal*lattice_C66(4,4,myPhase))/3.0_pReal +! c44bar = (lattice_C66(1,1,myPhase) - lattice_C66(1,2,myPhase) + lattice_C66(4,4,myPhase))/3.0_pReal +! c14bar = (lattice_C66(1,1,myPhase) - lattice_C66(1,2,myPhase) - 2.0_pReal*lattice_C66(4,4,myPhase)) & +! /(3.0_pReal*sqrt(2.0_pReal)) +! A = c14bar**(2.0_pReal)/c44bar +! B = c14bar**(2.0_pReal)/(0.5_pReal*(c11bar - c12bar)) +! temp66(1,1,myPhase) = c11bar - A +! temp66(1,2,myPhase) = c12bar + A +! temp66(1,3,myPhase) = c13bar +! temp66(3,3,myPhase) = c33bar +! temp66(4,4,myPhase) = c44bar - B +! +! temp66(1:6,1:6,myPhase) = lattice_symmetrizeC66(trans_lattice_structure(myPhase),& +! lattice_trans_C66(1:6,1:6,myPhase)) +! lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = math_Voigt66to3333(temp66(1:6,1:6,myPhase)) +! temp66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) +! do i = 1_pInt, 6_pInt +! if (abs(temp66(i,i,myPhase))hex transformation') +! enddo +! +!! Elasticity matrices for transformed phase +! select case(lattice_structure(myPhase)) +! case (LATTICE_fcc_ID) +! select case(trans_lattice_structure(myPhase)) +! case (LATTICE_bcc_ID) +! temp66(1:6,1:6,myPhase) = lattice_C66(1:6,1:6,myPhase) +! lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = lattice_C3333(1:3,1:3,1:3,1:3,myPhase) +! temp66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) +! do i = 1_pInt, 6_pInt +! if (abs(lattice_trans_C66(i,i,myPhase))bcc transformation') +! enddo +! case (LATTICE_hex_ID) +!select case(structure) +! case('fcc') +! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) +! case('bcc') +! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) +! case('hex','hexagonal') !ToDo: "No alias policy": long or short? +! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) +! case default +! call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_C66_twin)') +! end select +! do i = 1, sum(Ntwin) +! R = math_axisAngleToR(coordinateSystem(1:3,2,i), 180.0_pReal * INRAD) ! ToDo: Why always 180 deg? +! lattice_C66_twin(1:6,1:6,i) = math_Mandel3333to66(math_rotate_forward3333(math_Mandel66to3333(C66),R)) +! enddo +! +!end function + !function lattice_nonSchmidMatrix ! coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMSLIP,pInt),structure) ! @@ -2374,6 +2456,7 @@ end function buildInteraction !-------------------------------------------------------------------------------------------------- !> @brief build a local coordinate system in a slip, twin, trans, cleavage system +!> @details: Order: Direction, plane (normal), and common perpendicular !-------------------------------------------------------------------------------------------------- pure function buildCoordinateSystem(active,system,structure,cOverA) use math, only: & From 4e68d049a87bcfba1cac9212ebf6ba0417616f60 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 12 Sep 2018 20:37:55 +0200 Subject: [PATCH 10/45] WIP: calculating rotated stiffness matrices for transformation --- src/lattice.f90 | 170 +++++++++++++++++++++++++++--------------------- 1 file changed, 95 insertions(+), 75 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 28bd49f34..a8eacb2f3 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1524,7 +1524,6 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) case (LATTICE_fcc_ID) select case(trans_lattice_structure(myPhase)) case (LATTICE_bcc_ID) - temp66(1:6,1:6,myPhase) = lattice_C66(1:6,1:6,myPhase) lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = lattice_C3333(1:3,1:3,1:3,1:3,myPhase) temp66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) do i = 1_pInt, 6_pInt @@ -2116,84 +2115,105 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) end function -!function lattice_C66_trans(Ntrans,C66_parent,C66_targetstructure_parent,cOverA_parent,structure_target,cOverA_target) -! use IO, only: & -! IO_error -! use math, only: & -! INRAD, & -! math_axisAngleToR, & -! math_Mandel3333to66, & -! math_Mandel66to3333, & -! math_rotate_forward3333 +function lattice_C66_trans(Ntrans,C66_parent,structure_parent,cOverA_parent, & + C66_target,structure_target,cOverA_target) + use IO, only: & + IO_error + use math, only: & + INRAD, & + math_axisAngleToR, & + math_Mandel3333to66, & + math_Mandel66to3333, & + math_rotate_forward3333 + + implicit none + integer(pInt), dimension(:), intent(in) :: Ntrans !< number of active twin systems per family + character(len=*), intent(in) :: & + structure_target, & !< lattice structure + structure_parent !< lattice structure + real(pReal), dimension(6,6), intent(in) :: C66_target, C66_parent + real(pReal), intent(in) :: cOverA_parent, cOverA_target + real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans, C_bar66, C_target_unrotated66 + + real(pReal), dimension(3,3,sum(Ntrans)) :: coordinateSystem + + real(pReal), dimension(3,3) :: R + real(pReal) :: A, B + integer(pInt) :: i + +! if (trim(structure_parent) == 'fcc' .and. trim(structure_target) == 'hex') then +! C_bar66(1,1) = (C66_parent(1,1) + C66_parent(1,2) + 2.0_pReal*C66_parent(4,4))/2.0_pReal +! C_bar66(1,2) = (C66_parent(1,1) + 5.0_pReal*C66_parent(1,2) - 2.0_pReal*C66_parent(4,4))/6.0_pReal +! C_bar66(3,3) = (C66_parent(1,1) + 2.0_pReal*C66_parent(1,2) + 4.0_pReal*C66_parent(4,4))/3.0_pReal +! C_bar66(1,3) = (C66_parent(1,1) + 2.0_pReal*C66_parent(1,2) - 2.0_pReal*C66_parent(4,4))/3.0_pReal +! C_bar66(4,4) = (C66_parent(1,1) - C66_parent(1,2) + C66_parent(4,4))/3.0_pReal +! C_bar66(1,4) = (C66_parent(1,1) - C66_parent(1,2) - 2.0_pReal*C66_parent(4,4)) /(3.0_pReal*sqrt(2.0_pReal)) +! A = C_bar(1,4)**2.0_pReal/C_bar66(4,4) +! B = C_bar(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2))) +! C_target_unrotated66(1,1) = C_bar66(1,1) - A +! C_target_unrotated66(1,2) = C_bar66(1,2) + A +! C_target_unrotated66(1,3) = C_bar66(1,3) +! C_target_unrotated66(3,3) = C_bar66(3,3) +! C_target_unrotated66(4,4) = C_bar66(4,4) - B +! C_target_unrotated66 = lattice_symmetrizeC66(LATTICE_HEX_ID,C_target_unrotated) +! elseif (trim(structure_parent) == 'fcc' .and. trim(structure_target) == 'bcc') then +! C_target_unrotated66 = C_parent66 +! endif +! do i = 1_pInt, 6_pInt +! if (abs(C_target_unrotated66(i,i))hex transformation') -! enddo -! -!! Elasticity matrices for transformed phase -! select case(lattice_structure(myPhase)) -! case (LATTICE_fcc_ID) -! select case(trans_lattice_structure(myPhase)) -! case (LATTICE_bcc_ID) -! temp66(1:6,1:6,myPhase) = lattice_C66(1:6,1:6,myPhase) -! lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = lattice_C3333(1:3,1:3,1:3,1:3,myPhase) -! temp66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) -! do i = 1_pInt, 6_pInt -! if (abs(lattice_trans_C66(i,i,myPhase))bcc transformation') -! enddo -! case (LATTICE_hex_ID) -!select case(structure) -! case('fcc') -! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) -! case('bcc') -! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) -! case('hex','hexagonal') !ToDo: "No alias policy": long or short? -! coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) -! case default -! call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_C66_twin)') -! end select ! do i = 1, sum(Ntwin) ! R = math_axisAngleToR(coordinateSystem(1:3,2,i), 180.0_pReal * INRAD) ! ToDo: Why always 180 deg? -! lattice_C66_twin(1:6,1:6,i) = math_Mandel3333to66(math_rotate_forward3333(math_Mandel66to3333(C66),R)) +! lattice_C66_trans(1:6,1:6,i) = math_Mandel3333to66(math_rotate_forward3333(math_Mandel66to3333(C66),R)) ! enddo -! -!end function +! ! Phase transformation +! select case(trans_lattice_structure(myPhase)) +! case (LATTICE_bcc_ID) ! fcc to bcc transformation +! do i = 1_pInt,myNtrans +! Rtr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), & ! Pitsch rotation +! lattice_fccTobcc_systemTrans(4,i)*INRAD) +! Btr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_bainRot(1:3,i), & ! Rotation of fcc to Bain coordinate system +! lattice_fccTobcc_bainRot(4,i)*INRAD) +! xtr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(1:3,i),pReal) +! ytr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(4:6,i),pReal) +! ztr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(7:9,i),pReal) +! Utr(1:3,1:3,i) = 0.0_pReal ! Bain deformation +! if ((a_fcc > 0.0_pReal) .and. (a_bcc > 0.0_pReal)) then +! Utr(1:3,1:3,i) = (a_bcc/a_fcc)*math_tensorproduct33(xtr(1:3,i), xtr(1:3,i)) + & +! sqrt(2.0_pReal)*(a_bcc/a_fcc)*math_tensorproduct33(ytr(1:3,i), ytr(1:3,i)) + & +! sqrt(2.0_pReal)*(a_bcc/a_fcc)*math_tensorproduct33(ztr(1:3,i), ztr(1:3,i)) +! endif +! Qtr(1:3,1:3,i) = math_mul33x33(Rtr(1:3,1:3,i), Btr(1:3,1:3,i)) +! Str(1:3,1:3,i) = math_mul33x33(Rtr(1:3,1:3,i), Utr(1:3,1:3,i)) - MATH_I3 +! enddo +! case (LATTICE_hex_ID) +! sstr(1:3,1:3) = MATH_I3 +! sstr(1,3) = sqrt(2.0_pReal)/4.0_pReal +! sdtr(1:3,1:3) = MATH_I3 +! if (CoverA_trans > 1.0_pReal .and. CoverA_trans < 2.0_pReal) then +! sdtr(3,3) = CoverA_trans/sqrt(8.0_pReal/3.0_pReal) +! endif +! sttr = math_mul33x33(sdtr, sstr) +! do i = 1_pInt,myNtrans +! xtr(1:3,i) = lattice_fccTohex_systemTrans(1:3,i)/norm2(lattice_fccTohex_systemTrans(1:3,i)) +! ztr(1:3,i) = lattice_fccTohex_systemTrans(4:6,i)/norm2(lattice_fccTohex_systemTrans(4:6,i)) +! ytr(1:3,i) = -math_crossproduct(xtr(1:3,i), ztr(1:3,i)) +! Rtr(1:3,1,i) = xtr(1:3,i) +! Rtr(1:3,2,i) = ytr(1:3,i) +! Rtr(1:3,3,i) = ztr(1:3,i) +! Qtr(1:3,1:3,i) = Rtr(1:3,1:3,i) +! Str(1:3,1:3,i) = math_mul33x33(Rtr(1:3,1:3,i), math_mul33x33(sttr, math_transpose33(Rtr(1:3,1:3,i)))) +! Str(1:3,1:3,i) = Str(1:3,1:3,i) - MATH_I3 +! trs(i) = lattice_fccTohex_shearTrans(i) +! enddo +! case default +! Qtr = 0.0_pReal +! Str = 0.0_pReal +! end select + + end function !function lattice_nonSchmidMatrix ! coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMSLIP,pInt),structure) From cc20b044a6920c54dd2c3bf25b833e7d3e5294b0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 13 Sep 2018 09:27:54 +0200 Subject: [PATCH 11/45] forgot to allocate array now used internally only --- src/lattice.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/lattice.f90 b/src/lattice.f90 index a8eacb2f3..cbd4b2e06 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1264,6 +1264,7 @@ subroutine lattice_init allocate(lattice_structure(Nphases),source = LATTICE_undefined_ID) allocate(trans_lattice_structure(Nphases),source = LATTICE_undefined_ID) allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) + allocate(temp66(6,6,Nphases), source=0.0_pReal) allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal) allocate(lattice_trans_C66(6,6,Nphases), source=0.0_pReal) allocate(lattice_trans_C3333(3,3,3,3,Nphases), source=0.0_pReal) From 55511782e8ca57c4921045f74120d29a2fa0e34b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 14 Sep 2018 08:19:39 +0200 Subject: [PATCH 12/45] incomplete transition/renaming caused segmentation fault --- src/lattice.f90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index cbd4b2e06..ad3b39e39 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1133,7 +1133,7 @@ module lattice real(pReal), dimension(:,:,:), allocatable, private :: & temp66 real(pReal), dimension(:,:,:), allocatable, public, protected :: & - lattice_C66, lattice_trans_C66 + lattice_C66 real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & lattice_C3333, lattice_trans_C3333 real(pReal), dimension(:), allocatable, public, protected :: & @@ -1266,7 +1266,6 @@ subroutine lattice_init allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) allocate(temp66(6,6,Nphases), source=0.0_pReal) allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal) - allocate(lattice_trans_C66(6,6,Nphases), source=0.0_pReal) allocate(lattice_trans_C3333(3,3,3,3,Nphases), source=0.0_pReal) allocate(lattice_thermalExpansion33 (3,3,3,Nphases), source=0.0_pReal) ! constant, linear, quadratic coefficients allocate(lattice_thermalConductivity33 (3,3,Nphases), source=0.0_pReal) @@ -1528,7 +1527,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = lattice_C3333(1:3,1:3,1:3,1:3,myPhase) temp66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) do i = 1_pInt, 6_pInt - if (abs(lattice_trans_C66(i,i,myPhase))bcc transformation') enddo case (LATTICE_hex_ID) @@ -1548,7 +1547,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) temp66(4,4,myPhase) = c44bar - B temp66(1:6,1:6,myPhase) = lattice_symmetrizeC66(trans_lattice_structure(myPhase),& - lattice_trans_C66(1:6,1:6,myPhase)) + temp66(1:6,1:6,myPhase)) lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase) = math_Voigt66to3333(temp66(1:6,1:6,myPhase)) temp66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_trans_C3333(1:3,1:3,1:3,1:3,myPhase)) do i = 1_pInt, 6_pInt From ff2614b757de377c06bdacaed175415d560e7014 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 1 Oct 2018 22:48:14 +0200 Subject: [PATCH 13/45] using 'service functions' from lattice get rid of many global array on the long run --- src/lattice.f90 | 4 +++- src/plastic_phenopowerlaw.f90 | 16 ++++++++++------ 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index ad3b39e39..4e94a790f 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1224,7 +1224,9 @@ real(pReal), dimension(4,36), parameter, private :: & LATTICE_fcc_ID, & LATTICE_bcc_ID, & LATTICE_bct_ID, & - LATTICE_hex_ID + LATTICE_hex_ID, & + lattice_SchmidMatrix_slip, & + lattice_SchmidMatrix_twin contains diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index f7c723521..1dc37bd29 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -168,7 +168,8 @@ subroutine plastic_phenopowerlaw_init outputID !< ID of each post result output character(len=512) :: & - extmsg = '' + extmsg = '', & + structure = '' character(len=65536), dimension(:), allocatable :: outputs write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' @@ -193,6 +194,8 @@ subroutine plastic_phenopowerlaw_init associate(prm => param(instance),stt => state(instance),dot => dotState(instance)) extmsg = '' + structure = config_phase(p)%getString('lattice_structure') + prm%Nslip = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray) prm%totalNslip = sum(prm%Nslip) if (size(prm%Nslip) > count(lattice_NslipSystem(:,p) > 0_pInt)) & @@ -201,6 +204,9 @@ subroutine plastic_phenopowerlaw_init call IO_error(150_pInt,ext_msg='Nslip') slipActive: if (prm%totalNslip > 0_pInt) then + + prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& + config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) ! reading in slip related parameters prm%xi_slip_0 = config_phase(p)%getFloats('tau0_slip', requiredShape=shape(prm%Nslip)) prm%xi_slip_sat = config_phase(p)%getFloats('tausat_slip', requiredShape=shape(prm%Nslip)) @@ -242,6 +248,8 @@ subroutine plastic_phenopowerlaw_init call IO_error(150_pInt,ext_msg='Ntwin') twinActive: if (prm%totalNtwin > 0_pInt) then + prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,structure(1:3),& + config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) ! reading in twin related parameters prm%xi_twin_0 = config_phase(p)%getFloats('tau0_twin',requiredShape=shape(prm%Ntwin)) prm%interaction_TwinTwin = spread(config_phase(p)%getFloats('interaction_twintwin', & @@ -371,7 +379,6 @@ subroutine plastic_phenopowerlaw_init ! calculate hardening matrices allocate(temp1(prm%totalNslip,prm%totalNslip),source = 0.0_pReal) allocate(temp2(prm%totalNslip,prm%totalNtwin),source = 0.0_pReal) - allocate(prm%Schmid_slip(3,3,prm%totalNslip),source = 0.0_pReal) allocate(prm%nonSchmid_pos(3,3,size(prm%nonSchmidCoeff),prm%totalNslip),source = 0.0_pReal) allocate(prm%nonSchmid_neg(3,3,size(prm%nonSchmidCoeff),prm%totalNslip),source = 0.0_pReal) i = 0_pInt @@ -380,7 +387,6 @@ subroutine plastic_phenopowerlaw_init mySlipSystems: do j = 1_pInt,prm%Nslip(f) i = i + 1_pInt - prm%Schmid_slip(1:3,1:3,i) = lattice_Sslip(1:3,1:3,1,sum(lattice_Nslipsystem(1:f-1,p))+j,p) do k = 1,size(prm%nonSchmidCoeff) prm%nonSchmid_pos(1:3,1:3,k,i) = lattice_Sslip(1:3,1:3,2*k, index_myFamily+j,p) & * prm%nonSchmidCoeff(k) @@ -414,14 +420,12 @@ subroutine plastic_phenopowerlaw_init allocate(temp1(prm%totalNtwin,prm%totalNslip),source = 0.0_pReal) allocate(temp2(prm%totalNtwin,prm%totalNtwin),source = 0.0_pReal) - allocate(prm%Schmid_twin(3,3,prm%totalNtwin),source = 0.0_pReal) - allocate(prm%gamma_twin_char(prm%totalNtwin),source = 0.0_pReal) + allocate(prm%gamma_twin_char(prm%totalNtwin),source = 0.0_pReal) i = 0_pInt myTwinFamilies: do f = 1_pInt,size(prm%Ntwin,1) ! >>> interaction twin -- X index_myFamily = sum(prm%Ntwin(1:f-1_pInt)) myTwinSystems: do j = 1_pInt,prm%Ntwin(f) i = i + 1_pInt - prm%Schmid_twin(1:3,1:3,i) = lattice_Stwin(1:3,1:3,sum(lattice_NTwinsystem(1:f-1,p))+j,p) prm%gamma_twin_char(i) = lattice_shearTwin(sum(lattice_Ntwinsystem(1:f-1,p))+j,p) slipFamilies: do o = 1_pInt,size(prm%Nslip,1) index_otherFamily = sum(prm%Nslip(1:o-1_pInt)) From 384a7858050bb00ca3587ebc91017a6c713fe428 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 1 Oct 2018 23:02:31 +0200 Subject: [PATCH 14/45] more sevice functions use still trying --- src/lattice.f90 | 4 +++- src/plastic_phenopowerlaw.f90 | 38 +++++++---------------------------- 2 files changed, 10 insertions(+), 32 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 4e94a790f..be628b299 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1226,7 +1226,9 @@ real(pReal), dimension(4,36), parameter, private :: & LATTICE_bct_ID, & LATTICE_hex_ID, & lattice_SchmidMatrix_slip, & - lattice_SchmidMatrix_twin + lattice_SchmidMatrix_twin, & + lattice_interactionSlipSlip2, & + lattice_interactionTwinTwin2 contains diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 1dc37bd29..920f2e028 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -152,7 +152,7 @@ subroutine plastic_phenopowerlaw_init sizeState,sizeDotState, & startIndex, endIndex - real(pReal), dimension(:,:), allocatable :: temp1, temp2 + real(pReal), dimension(:,:), allocatable :: temp1 integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] @@ -210,8 +210,8 @@ subroutine plastic_phenopowerlaw_init ! reading in slip related parameters prm%xi_slip_0 = config_phase(p)%getFloats('tau0_slip', requiredShape=shape(prm%Nslip)) prm%xi_slip_sat = config_phase(p)%getFloats('tausat_slip', requiredShape=shape(prm%Nslip)) - prm%interaction_SlipSlip = spread(config_phase(p)%getFloats('interaction_slipslip', & - requiredShape=shape(prm%Nslip)),2,1) + prm%interaction_SlipSlip = lattice_interactionSlipSlip2(prm%Nslip,config_phase(p)%getFloats('interaction_slipslip'), & + structure(1:3)) prm%H_int = config_phase(p)%getFloats('h_int', requiredShape=shape(prm%Nslip), & defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) prm%nonSchmidCoeff = config_phase(p)%getFloats('nonschmid_coefficients',& @@ -252,8 +252,8 @@ subroutine plastic_phenopowerlaw_init config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) ! reading in twin related parameters prm%xi_twin_0 = config_phase(p)%getFloats('tau0_twin',requiredShape=shape(prm%Ntwin)) - prm%interaction_TwinTwin = spread(config_phase(p)%getFloats('interaction_twintwin', & - requiredShape=shape(prm%Ntwin)),2,1) + prm%interaction_TwinTwin = lattice_interactionTwinTwin2(prm%Ntwin,config_phase(p)%getFloats('interaction_twintwin'), & + structure(1:3)) prm%gdot0_twin = config_phase(p)%getFloat('gdot0_twin') prm%n_twin = config_phase(p)%getFloat('n_twin') @@ -378,7 +378,6 @@ subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! calculate hardening matrices allocate(temp1(prm%totalNslip,prm%totalNslip),source = 0.0_pReal) - allocate(temp2(prm%totalNslip,prm%totalNtwin),source = 0.0_pReal) allocate(prm%nonSchmid_pos(3,3,size(prm%nonSchmidCoeff),prm%totalNslip),source = 0.0_pReal) allocate(prm%nonSchmid_neg(3,3,size(prm%nonSchmidCoeff),prm%totalNslip),source = 0.0_pReal) i = 0_pInt @@ -393,20 +392,10 @@ subroutine plastic_phenopowerlaw_init prm%nonSchmid_neg(1:3,1:3,k,i) = lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+j,p) & * prm%nonSchmidCoeff(k) enddo - otherSlipFamilies: do o = 1_pInt,size(prm%Nslip,1) - index_otherFamily = sum(prm%Nslip(1:o-1_pInt)) - otherSlipSystems: do k = 1_pInt,prm%Nslip(o) - temp1(index_myFamily+j,index_otherFamily+k) = & - prm%interaction_SlipSlip(lattice_interactionSlipSlip( & - sum(lattice_NslipSystem(1:f-1,p))+j, & - sum(lattice_NslipSystem(1:o-1,p))+k, & - p),1) - enddo otherSlipSystems; enddo otherSlipFamilies - twinFamilies: do o = 1_pInt,size(prm%Ntwin,1) index_otherFamily = sum(prm%Ntwin(1:o-1_pInt)) twinSystems: do k = 1_pInt,prm%Ntwin(o) - temp2(index_myFamily+j,index_otherFamily+k) = & + temp1(index_myFamily+j,index_otherFamily+k) = & prm%interaction_SlipTwin(lattice_interactionSlipTwin( & sum(lattice_NslipSystem(1:f-1_pInt,p))+j, & sum(lattice_NtwinSystem(1:o-1_pInt,p))+k, & @@ -414,12 +403,10 @@ subroutine plastic_phenopowerlaw_init enddo twinSystems; enddo twinFamilies enddo mySlipSystems enddo mySlipFamilies - prm%interaction_SlipSlip = temp1; deallocate(temp1) - prm%interaction_SlipTwin = temp2; deallocate(temp2) + prm%interaction_SlipTwin = temp1; deallocate(temp1) allocate(temp1(prm%totalNtwin,prm%totalNslip),source = 0.0_pReal) - allocate(temp2(prm%totalNtwin,prm%totalNtwin),source = 0.0_pReal) allocate(prm%gamma_twin_char(prm%totalNtwin),source = 0.0_pReal) i = 0_pInt myTwinFamilies: do f = 1_pInt,size(prm%Ntwin,1) ! >>> interaction twin -- X @@ -436,20 +423,9 @@ subroutine plastic_phenopowerlaw_init sum(lattice_NslipSystem(1:o-1_pInt,p))+k, & p),1) enddo slipSystems; enddo slipFamilies - - otherTwinFamilies: do o = 1_pInt,size(prm%Ntwin,1) - index_otherFamily = sum(prm%Ntwin(1:o-1_pInt)) - otherTwinSystems: do k = 1_pInt,prm%Ntwin(o) - temp2(index_myFamily+j,index_otherFamily+k) = & - prm%interaction_TwinTwin(lattice_interactionTwinTwin( & - sum(lattice_NtwinSystem(1:f-1_pInt,p))+j, & - sum(lattice_NtwinSystem(1:o-1_pInt,p))+k, & - p),1) - enddo otherTwinSystems; enddo otherTwinFamilies enddo myTwinSystems enddo myTwinFamilies prm%interaction_TwinSlip = temp1; deallocate(temp1) - prm%interaction_TwinTwin = temp2; deallocate(temp2) !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState From 67c55965e9cd7b656c7d64aad989eaea8810aaa6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Oct 2018 06:45:42 +0200 Subject: [PATCH 15/45] segmentation fault, tried LHS with scalar Nxx per system is stored in LATTICE_yyy_NxxxSYSTEM Using CAPITALS for parameters improves readbility --- src/lattice.f90 | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index be628b299..8e4e79d93 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2264,17 +2264,17 @@ function lattice_interactionSlipSlip2(Nslip,interactionValues,structure) select case(structure) case('fcc') - interactionSlipSlip = lattice_fcc_interactionSlipSlip - NslipMax = lattice_fcc_Nslip + interactionSlipSlip = LATTICE_FCC_INTERACTIONSLIPSLIP + NslipMax = LATTICE_FCC_NSLIPSYSTEM case('bcc') - interactionSlipSlip = lattice_bcc_interactionSlipSlip - NslipMax = lattice_bcc_Nslip + interactionSlipSlip = LATTICE_BCC_INTERACTIONSLIPSLIP + NslipMax = LATTICE_BCC_NSLIPSYSTEM case('hex','hexagonal') !ToDo: "No alias policy": long or short? - interactionSlipSlip = lattice_hex_interactionTwinTwin - NslipMax = lattice_hex_Nslip + interactionSlipSlip = LATTICE_HEX_INTERACTIONSLIPSLIP + NslipMax = LATTICE_HEX_NSLIPSYSTEM case('bct') - interactionSlipSlip = lattice_bct_interactionSlipSlip - NslipMax = lattice_bct_Nslip + interactionSlipSlip = LATTICE_BCT_INTERACTIONSLIPSLIP + NslipMax = LATTICE_BCT_NSLIPSYSTEM case default call IO_error(132_pInt,ext_msg=trim(structure)//' (slip slip interaction)') end select @@ -2307,14 +2307,14 @@ function lattice_interactionTwinTwin2(Ntwin,interactionValues,structure) select case(structure) case('fcc') - interactionTwinTwin = lattice_fcc_interactionTwinTwin - NtwinMax = lattice_fcc_Ntwin + interactionTwinTwin = LATTICE_FCC_INTERACTIONTWINTWIN + NtwinMax = LATTICE_FCC_NTWINSYSTEM case('bcc') - interactionTwinTwin = lattice_bcc_interactionTwinTwin - NtwinMax = lattice_bcc_Ntwin + interactionTwinTwin = LATTICE_BCC_INTERACTIONTWINTWIN + NtwinMax = LATTICE_BCC_NTWINSYSTEM case('hex','hexagonal') !ToDo: "No alias policy": long or short? - interactionTwinTwin = lattice_hex_interactionTwinTwin - NtwinMax = lattice_hex_Ntwin + interactionTwinTwin = LATTICE_HEX_INTERACTIONTWINTWIN + NtwinMax = LATTICE_HEX_NTWINSYSTEM case default call IO_error(132_pInt,ext_msg=trim(structure)//' (twin twin interaction)') end select From d718521890c04754490bc25c6bce2f0302c79b7c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Oct 2018 07:13:30 +0200 Subject: [PATCH 16/45] polishing --- src/plastic_phenopowerlaw.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 920f2e028..b4fecda7b 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -349,7 +349,7 @@ subroutine plastic_phenopowerlaw_init NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase sizeState = size(['tau_slip ','gamma_slip']) * prm%TotalNslip & + size(['tau_twin ','gamma_twin']) * prm%TotalNtwin & - + size(['sum(gamma)', 'sum(f) ']) + + size(['sum(gamma)','sum(f) ']) sizeDotState = sizeState plasticState(p)%sizeState = sizeState @@ -551,7 +551,7 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) of integer(pInt) :: & - i,k + i real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & xi_slip_sat_offset @@ -643,7 +643,6 @@ subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, & real(pReal), dimension(prm%totalNslip) :: & tau_slip_pos, & tau_slip_neg - integer(pInt) :: i, j do i = 1_pInt, prm%totalNslip From ba383336f5b8382584e05e277f74f6724dae6983 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Oct 2018 07:23:18 +0200 Subject: [PATCH 17/45] only one projection matrix required for all (non)schmid projections following the example from plastic_nonlocal --- src/plastic_phenopowerlaw.f90 | 35 ++++++++++++++++------------------- 1 file changed, 16 insertions(+), 19 deletions(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index b4fecda7b..65fb7c513 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -70,8 +70,7 @@ module plastic_phenopowerlaw interaction_TwinTwin !< twin resistance from twin activity real(pReal), dimension(:,:,:), allocatable :: & Schmid_slip, & - Schmid_twin - real(pReal), dimension(:,:,:,:), allocatable :: & + Schmid_twin, & nonSchmid_pos, & nonSchmid_neg integer(kind(undefined_ID)), dimension(:), allocatable :: & @@ -378,19 +377,23 @@ subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! calculate hardening matrices allocate(temp1(prm%totalNslip,prm%totalNslip),source = 0.0_pReal) - allocate(prm%nonSchmid_pos(3,3,size(prm%nonSchmidCoeff),prm%totalNslip),source = 0.0_pReal) - allocate(prm%nonSchmid_neg(3,3,size(prm%nonSchmidCoeff),prm%totalNslip),source = 0.0_pReal) + allocate(prm%nonSchmid_pos(3,3,prm%totalNslip),source = 0.0_pReal) + allocate(prm%nonSchmid_neg(3,3,prm%totalNslip),source = 0.0_pReal) i = 0_pInt mySlipFamilies: do f = 1_pInt,size(prm%Nslip,1) ! >>> interaction slip -- X index_myFamily = sum(prm%Nslip(1:f-1_pInt)) mySlipSystems: do j = 1_pInt,prm%Nslip(f) i = i + 1_pInt + prm%nonSchmid_pos(1:3,1:3,i) = lattice_Sslip(1:3,1:3,1, index_myFamily+j,p) + prm%nonSchmid_neg(1:3,1:3,i) = lattice_Sslip(1:3,1:3,1, index_myFamily+j,p) do k = 1,size(prm%nonSchmidCoeff) - prm%nonSchmid_pos(1:3,1:3,k,i) = lattice_Sslip(1:3,1:3,2*k, index_myFamily+j,p) & - * prm%nonSchmidCoeff(k) - prm%nonSchmid_neg(1:3,1:3,k,i) = lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+j,p) & - * prm%nonSchmidCoeff(k) + prm%nonSchmid_pos(1:3,1:3,i) = prm%nonSchmid_pos(1:3,1:3,i) & + + lattice_Sslip(1:3,1:3,2*k, index_myFamily+j,p) & + * prm%nonSchmidCoeff(k) + prm%nonSchmid_neg(1:3,1:3,i) = prm%nonSchmid_neg(1:3,1:3,i) & + + lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+j,p) & + * prm%nonSchmidCoeff(k) enddo twinFamilies: do o = 1_pInt,size(prm%Ntwin,1) index_otherFamily = sum(prm%Ntwin(1:o-1_pInt)) @@ -517,12 +520,10 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) Lp = Lp + (1.0_pReal-stt%sumF(of))*(gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) & - *(prm%Schmid_slip(m,n,i) + sum(prm%nonSchmid_pos(m,n,:,i))) + + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) & - *(prm%Schmid_slip(m,n,i) + sum(prm%nonSchmid_neg(m,n,:,i))) + + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems call kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtautwin) @@ -643,15 +644,11 @@ subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, & real(pReal), dimension(prm%totalNslip) :: & tau_slip_pos, & tau_slip_neg - integer(pInt) :: i, j + integer(pInt) :: i do i = 1_pInt, prm%totalNslip - tau_slip_pos(i) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i)) - tau_slip_neg(i) = tau_slip_pos(i) - do j = 1,size(prm%nonSchmidCoeff) - tau_slip_pos(i) = tau_slip_pos(i) + math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j,i)) - tau_slip_neg(i) = tau_slip_neg(i) + math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j,i)) - enddo + tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) + tau_slip_neg(i) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) enddo gdot_slip_pos = 0.5_pReal*prm%gdot0_slip & From 0cac6ce0d2a2953b5425b149a876b9028ec913ac Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Oct 2018 07:36:57 +0200 Subject: [PATCH 18/45] allocating to wrong size --- src/plastic_phenopowerlaw.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 65fb7c513..e836cb3fa 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -376,7 +376,7 @@ subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! calculate hardening matrices - allocate(temp1(prm%totalNslip,prm%totalNslip),source = 0.0_pReal) + allocate(temp1(prm%totalNslip,prm%totalNtwin),source = 0.0_pReal) allocate(prm%nonSchmid_pos(3,3,prm%totalNslip),source = 0.0_pReal) allocate(prm%nonSchmid_neg(3,3,prm%totalNslip),source = 0.0_pReal) i = 0_pInt From 89e03be69942ee582358422129bf46e16c35fed4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Oct 2018 07:51:11 +0200 Subject: [PATCH 19/45] also building interactions for slip-twin and twin-slip --- src/lattice.f90 | 90 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) diff --git a/src/lattice.f90 b/src/lattice.f90 index 8e4e79d93..09476ae4d 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2328,6 +2328,96 @@ function lattice_interactionTwinTwin2(Ntwin,interactionValues,structure) end function lattice_interactionTwinTwin2 +!-------------------------------------------------------------------------------------------------- +!> @brief Populates reduced slip-twin interaction matrix +!> ToDo: prefix "2" needed as long as deprecated array lattice_interactionTwinTwin exists +!-------------------------------------------------------------------------------------------------- +function lattice_interactionSlipTwin2(Nslip,Ntwin,interactionValues,structure) + use IO, only: & + IO_error + + implicit none + integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family + integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family + real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin + character(len=*), intent(in) :: structure !< lattice structure + real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: lattice_interactionSlipTwin2 + + integer(pInt), dimension(:), allocatable :: NslipMax + integer(pInt), dimension(:), allocatable :: NtwinMax + integer(pInt), dimension(:,:), allocatable :: interactionSlipTwin + + select case(structure) + case('fcc') + interactionSlipTwin = LATTICE_FCC_INTERACTIONSLIPTWIN + NslipMax = LATTICE_FCC_NSLIPSYSTEM + NtwinMax = LATTICE_FCC_NTWINSYSTEM + case('bcc') + interactionSlipTwin = LATTICE_BCC_INTERACTIONSLIPTWIN + NslipMax = LATTICE_BCC_NSLIPSYSTEM + NtwinMax = LATTICE_BCC_NTWINSYSTEM + case('hex','hexagonal') !ToDo: "No alias policy": long or short? + interactionSlipTwin = LATTICE_HEX_INTERACTIONSLIPTWIN + NslipMax = LATTICE_HEX_NSLIPSYSTEM + NtwinMax = LATTICE_HEX_NTWINSYSTEM + case default + call IO_error(132_pInt,ext_msg=trim(structure)//' (slip twin interaction)') + end select + + if (any(Ntwin(1:size(Ntwin)) - Ntwin < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) + + lattice_interactionSlipTwin2 = & + buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionSlipTwin) + +end function lattice_interactionSlipTwin2 + + +!-------------------------------------------------------------------------------------------------- +!> @brief Populates reduced slip-twin interaction matrix +!> ToDo: prefix "2" needed as long as deprecated array lattice_interactionTwinTwin exists +!-------------------------------------------------------------------------------------------------- +function lattice_interactionTwinSlip2(Ntwin,Nslip,interactionValues,structure) + use IO, only: & + IO_error + + implicit none + integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family + integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family + real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin + character(len=*), intent(in) :: structure !< lattice structure + real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: lattice_interactionTwinSlip2 + + integer(pInt), dimension(:), allocatable :: NslipMax + integer(pInt), dimension(:), allocatable :: NtwinMax + integer(pInt), dimension(:,:), allocatable :: interactionTwinSlip + + select case(structure) + case('fcc') + interactionTwinSlip = LATTICE_FCC_INTERACTIONTWINSLIP + NtwinMax = LATTICE_FCC_NTWINSYSTEM + NslipMax = LATTICE_FCC_NSLIPSYSTEM + case('bcc') + interactionTwinSlip = LATTICE_BCC_INTERACTIONTWINSLIP + NtwinMax = LATTICE_BCC_NTWINSYSTEM + NslipMax = LATTICE_BCC_NSLIPSYSTEM + case('hex','hexagonal') !ToDo: "No alias policy": long or short? + interactionTwinSlip = LATTICE_HEX_INTERACTIONTWINSLIP + NtwinMax = LATTICE_HEX_NTWINSYSTEM + NslipMax = LATTICE_HEX_NSLIPSYSTEM + case default + call IO_error(132_pInt,ext_msg=trim(structure)//' (slip twin interaction)') + end select + + if (any(Ntwin(1:size(Ntwin)) - Ntwin < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) + + lattice_interactionTwinSlip2 = & + buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTwinSlip) + +end function lattice_interactionTwinSlip2 + + !-------------------------------------------------------------------------------------------------- !> @brief Populates reduced trans-trans interaction matrix !> ToDo: prefix "2" needed as long as deprecated array lattice_interactionTransTrans exists From 103b46ecbb6e579d2e479b2f93e4f78958edacf3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Oct 2018 07:58:02 +0200 Subject: [PATCH 20/45] don't repeate code over different plastic_xxx modules --- src/lattice.f90 | 4 +++- src/plastic_phenopowerlaw.f90 | 35 ++++++++--------------------------- 2 files changed, 11 insertions(+), 28 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 09476ae4d..c970f7e91 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1228,7 +1228,9 @@ real(pReal), dimension(4,36), parameter, private :: & lattice_SchmidMatrix_slip, & lattice_SchmidMatrix_twin, & lattice_interactionSlipSlip2, & - lattice_interactionTwinTwin2 + lattice_interactionTwinTwin2, & + lattice_interactionSlipTwin2, & + lattice_interactionTwinSlip2 contains diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index e836cb3fa..519d2ac33 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -151,8 +151,6 @@ subroutine plastic_phenopowerlaw_init sizeState,sizeDotState, & startIndex, endIndex - real(pReal), dimension(:,:), allocatable :: temp1 - integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] @@ -272,8 +270,12 @@ subroutine plastic_phenopowerlaw_init endif twinActive slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then - prm%interaction_SlipTwin = spread(config_phase(p)%getFloats('interaction_sliptwin'),2,1) - prm%interaction_TwinSlip = spread(config_phase(p)%getFloats('interaction_twinslip'),2,1) + prm%interaction_SlipTwin = lattice_interactionSlipTwin2(prm%Nslip,prm%Ntwin,& + config_phase(p)%getFloats('interaction_sliptwin'), & + structure(1:3)) + prm%interaction_TwinSlip = lattice_interactionTwinSlip2(prm%Ntwin,prm%Nslip,& + config_phase(p)%getFloats('interaction_twinslip'), & + structure(1:3)) prm%h0_TwinSlip = config_phase(p)%getFloat('h0_twinslip') else slipAndTwinActive prm%h0_TwinSlip = 0.0_pReal @@ -350,6 +352,8 @@ subroutine plastic_phenopowerlaw_init + size(['tau_twin ','gamma_twin']) * prm%TotalNtwin & + size(['sum(gamma)','sum(f) ']) +!-------------------------------------------------------------------------------------------------- +! ToDo: This could be done by a function (in constitutive?) sizeDotState = sizeState plasticState(p)%sizeState = sizeState plasticState(p)%sizeDotState = sizeDotState @@ -376,7 +380,6 @@ subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! calculate hardening matrices - allocate(temp1(prm%totalNslip,prm%totalNtwin),source = 0.0_pReal) allocate(prm%nonSchmid_pos(3,3,prm%totalNslip),source = 0.0_pReal) allocate(prm%nonSchmid_neg(3,3,prm%totalNslip),source = 0.0_pReal) i = 0_pInt @@ -395,21 +398,9 @@ subroutine plastic_phenopowerlaw_init + lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+j,p) & * prm%nonSchmidCoeff(k) enddo - twinFamilies: do o = 1_pInt,size(prm%Ntwin,1) - index_otherFamily = sum(prm%Ntwin(1:o-1_pInt)) - twinSystems: do k = 1_pInt,prm%Ntwin(o) - temp1(index_myFamily+j,index_otherFamily+k) = & - prm%interaction_SlipTwin(lattice_interactionSlipTwin( & - sum(lattice_NslipSystem(1:f-1_pInt,p))+j, & - sum(lattice_NtwinSystem(1:o-1_pInt,p))+k, & - p),1) - enddo twinSystems; enddo twinFamilies enddo mySlipSystems enddo mySlipFamilies - prm%interaction_SlipTwin = temp1; deallocate(temp1) - - allocate(temp1(prm%totalNtwin,prm%totalNslip),source = 0.0_pReal) allocate(prm%gamma_twin_char(prm%totalNtwin),source = 0.0_pReal) i = 0_pInt myTwinFamilies: do f = 1_pInt,size(prm%Ntwin,1) ! >>> interaction twin -- X @@ -417,18 +408,8 @@ subroutine plastic_phenopowerlaw_init myTwinSystems: do j = 1_pInt,prm%Ntwin(f) i = i + 1_pInt prm%gamma_twin_char(i) = lattice_shearTwin(sum(lattice_Ntwinsystem(1:f-1,p))+j,p) - slipFamilies: do o = 1_pInt,size(prm%Nslip,1) - index_otherFamily = sum(prm%Nslip(1:o-1_pInt)) - slipSystems: do k = 1_pInt,prm%Nslip(o) - temp1(index_myFamily+j,index_otherFamily+k) = & - prm%interaction_TwinSlip(lattice_interactionTwinSlip( & - sum(lattice_NtwinSystem(1:f-1_pInt,p))+j, & - sum(lattice_NslipSystem(1:o-1_pInt,p))+k, & - p),1) - enddo slipSystems; enddo slipFamilies enddo myTwinSystems enddo myTwinFamilies - prm%interaction_TwinSlip = temp1; deallocate(temp1) !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState From db645efda7c8a9c4ee6d31acf1612514ff88eb17 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Oct 2018 08:49:23 +0200 Subject: [PATCH 21/45] function to calculate nonSchmid matrix --- src/lattice.f90 | 70 +++++++++++++++++++++++++++++++------------------ 1 file changed, 44 insertions(+), 26 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index c970f7e91..7d70396d5 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1304,7 +1304,7 @@ subroutine lattice_init allocate(lattice_Scleavage(3,3,3,lattice_maxNslip,Nphases),source=0.0_pReal) allocate(lattice_Scleavage_v(6,3,lattice_maxNslip,Nphases),source=0.0_pReal) - allocate(lattice_Qtwin(3,3,lattice_maxNtwin,Nphases),source=0.0_pReal) + allocate(lattice_Qtwin(3,3,lattice_maxNtwin,Nphases),source=0.0_pReal) allocate(lattice_Stwin(3,3,lattice_maxNtwin,Nphases),source=0.0_pReal) allocate(lattice_Stwin_v(6,lattice_maxNtwin,Nphases),source=0.0_pReal) @@ -2081,6 +2081,7 @@ pure function lattice_qDisorientation(Q1, Q2, struct) end function lattice_qDisorientation + function lattice_C66_twin(Ntwin,C66,structure,CoverA) use IO, only: & IO_error @@ -2221,31 +2222,48 @@ function lattice_C66_trans(Ntrans,C66_parent,structure_parent,cOverA_parent, & end function -!function lattice_nonSchmidMatrix -! coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMSLIP,pInt),structure) -! -! do i = 1_pInt,myNslip ! assign slip system vectors -! direction = coordinateSystem(1:3,1,i) -! normal = coordinateSystem(1:3,1,i) -! ! "np" and "nn" according to Gröger_etal2008, Acta Materialia 56 (2008) 5412–5425, table 1 (corresponds to their "n1" for positive and negative slip direction respectively) -! np = math_mul33x3(math_axisAngleToR(sdU,+60.0_pReal*INRAD), normal) ! ! Schmid matrices with non-Schmid contributions according to Koester_etal2012, Acta Materialia 60 (2012) 3894–3901, eq. (17) ("n1" is replaced by either "np" or "nn" according to either positive or negative slip direction) -! nonSchmidMatrix_pos(1:3,1:3,1,i) = math_tensorproduct33(direction, np) -! nonSchmidMatrix_pos(1:3,1:3,2,i) = math_tensorproduct33(math_crossproduct(normal, direction), normal) -! nonSchmidMatrix_pos(1:3,1:3,3,i) = math_tensorproduct33(math_crossproduct(np, direction), np) -! nonSchmidMatrix_pos(1:3,1:3,4,i) = math_tensorproduct33(normal, normal) -! nonSchmidMatrix_pos(1:3,1:3,5,i) = math_tensorproduct33(math_crossproduct(normal, direction), math_crossproduct(normal, direction)) -! nonSchmidMatrix_pos(1:3,1:3,6,i) = math_tensorproduct33(direction, direction) -! -! nn = math_mul33x3(math_axisAngleToR(-sdU,60.0_pReal*INRAD), normal) -! sns(1:3,1:3,2,1,i) = math_tensorproduct33(-sdU, nn) -! sns(1:3,1:3,2,2,i) = math_tensorproduct33(math_crossproduct(snU, -sdU), snU) -! sns(1:3,1:3,2,3,i) = math_tensorproduct33(math_crossproduct(nn, -sdU), nn) -! sns(1:3,1:3,2,4,i) = math_tensorproduct33(snU, snU) -! sns(1:3,1:3,2,5,i) = math_tensorproduct33(math_crossproduct(snU, -sdU), math_crossproduct(snU, -sdU)) -! sns(1:3,1:3,2,6,i) = math_tensorproduct33(-sdU, -sdU) -! enddo -!end function lattice_nonSchmidMatrix +! ! "np" and "nn" according to Gröger_etal2008, Acta Materialia 56 (2008) 5412–5425, table 1 (corresponds to their "n1" for positive and negative slip direction respectively) +function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix) + use math, only: & + INRAD, & + math_tensorproduct33, & + math_crossproduct, & + math_mul33x3, & + math_axisAngleToR + implicit none + integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family + real(pReal), dimension(6), intent(in) :: nonSchmidCoefficients + integer(pInt), intent(in) :: sense !< sense (-1,+1) + + real(pReal), dimension(1:3,1:3,sum(Nslip)) :: nonSchmidMatrix + + real(pReal), dimension(1:3,1:3,sum(Nslip)) :: coordinateSystem + real(pReal), dimension(:), allocatable :: direction + real(pReal), dimension(:), allocatable :: normal,np + integer(pInt) :: i + + if (abs(sense) /= 1_pInt) write(6,*) 'mist' + coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMSLIP,pInt),'bcc') + coordinateSystem(1:3,1,sum(Nslip)) = coordinateSystem(1:3,1,sum(Nslip)) *real(sense,pReal) + nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'bcc') + + do i = 1_pInt,sum(Nslip) + direction = coordinateSystem(1:3,1,i) + normal = coordinateSystem(1:3,2,i) + np = math_mul33x3(math_axisAngleToR(+direction,60.0_pReal*INRAD), normal) + nonSchmidMatrix(1:3,1:3,i) & + = nonSchmidMatrix(1:3,1:3,i) & + + nonSchmidCoefficients(1) * math_tensorproduct33(direction, np) & + + nonSchmidCoefficients(2) * math_tensorproduct33(math_crossproduct(normal, direction), normal) & + + nonSchmidCoefficients(3) * math_tensorproduct33(math_crossproduct(np, direction), np) & + + nonSchmidCoefficients(4) * math_tensorproduct33(normal, normal) & + + nonSchmidCoefficients(5) * math_tensorproduct33(math_crossproduct(normal, direction), & + math_crossproduct(normal, normal)) & + + nonSchmidCoefficients(6) * math_tensorproduct33(direction, direction) + enddo +end function lattice_nonSchmidMatrix + !-------------------------------------------------------------------------------------------------- !> @brief Populates reduced slip-slip interaction matrix @@ -2376,7 +2394,7 @@ end function lattice_interactionSlipTwin2 !-------------------------------------------------------------------------------------------------- -!> @brief Populates reduced slip-twin interaction matrix +!> @brief Populates reduced twin-slip interaction matrix !> ToDo: prefix "2" needed as long as deprecated array lattice_interactionTwinTwin exists !-------------------------------------------------------------------------------------------------- function lattice_interactionTwinSlip2(Ntwin,Nslip,interactionValues,structure) From ec77d4450b91d6541157bc2e79731c71c0097b3b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Oct 2018 09:17:06 +0200 Subject: [PATCH 22/45] checking for required things --- src/lattice.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 7d70396d5..090dad385 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -13,13 +13,13 @@ module lattice implicit none private +! BEGIN DEPRECATED integer(pInt), parameter, public :: & LATTICE_maxNslipFamily = 13_pInt, & !< max # of slip system families over lattice structures LATTICE_maxNtwinFamily = 4_pInt, & !< max # of twin system families over lattice structures LATTICE_maxNtransFamily = 1_pInt, & !< max # of transformation system families over lattice structures LATTICE_maxNcleavageFamily = 3_pInt !< max # of transformation system families over lattice structures -! BEGIN DEPRECATED integer(pInt), allocatable, dimension(:,:), protected, public :: & lattice_NslipSystem, & !< total # of slip systems in each family lattice_NtwinSystem, & !< total # of twin systems in each family @@ -1091,6 +1091,7 @@ module lattice 1, 0, 0, 0, 0, 1 & ],pReal),[ 3_pInt + 3_pInt,LATTICE_ortho_Ncleavage]) +! BEGIN DEPRECATED integer(pInt), parameter, public :: & LATTICE_maxNslip = max(LATTICE_fcc_Nslip,LATTICE_bcc_Nslip,LATTICE_hex_Nslip, & LATTICE_bct_Nslip,LATTICE_iso_Nslip,LATTICE_ortho_Nslip), & !< max # of slip systems over lattice structures @@ -1130,6 +1131,7 @@ module lattice !maxval(lattice_bct_interactionTwinTwin))) ) !< max # of interaction types (in hardening matrix part) #endif +!END DEPRECATED real(pReal), dimension(:,:,:), allocatable, private :: & temp66 real(pReal), dimension(:,:,:), allocatable, public, protected :: & From f494803c76eb570ac79fd96fe1e36af6a67aa99f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Oct 2018 16:40:40 +0200 Subject: [PATCH 23/45] fixed tangent calculation --- src/plastic_phenopowerlaw.f90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 519d2ac33..46902eb60 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -501,10 +501,8 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) Lp = Lp + (1.0_pReal-stt%sumF(of))*(gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) + + dgdot_dtauslip_pos(i) * (prm%Schmid_slip(k,l,i) + prm%nonSchmid_pos(m,n,i)) & + + dgdot_dtauslip_neg(i) * (prm%Schmid_slip(k,l,i) + prm%nonSchmid_neg(m,n,i)) enddo slipSystems call kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtautwin) From cd7f285686dc1a4a20931fedb123ba74ba561016 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 4 Oct 2018 07:17:00 +0200 Subject: [PATCH 24/45] prevent use of unallocated array interaction matrices are used even if no slip/twin systems defined --- src/plastic_phenopowerlaw.f90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 46902eb60..e9663ad4e 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -234,6 +234,7 @@ subroutine plastic_phenopowerlaw_init prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip) prm%H_int = math_expand(prm%H_int,prm%Nslip) else slipActive + allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%xi_slip_0(0)) endif slipActive @@ -266,6 +267,7 @@ subroutine plastic_phenopowerlaw_init ! expand slip related parameters from system => family prm%xi_twin_0 = math_expand(prm%xi_twin_0,prm%Ntwin) else twinActive + allocate(prm%interaction_TwinTwin(0,0)) allocate(prm%xi_twin_0(0)) endif twinActive @@ -276,8 +278,9 @@ subroutine plastic_phenopowerlaw_init prm%interaction_TwinSlip = lattice_interactionTwinSlip2(prm%Ntwin,prm%Nslip,& config_phase(p)%getFloats('interaction_twinslip'), & structure(1:3)) - prm%h0_TwinSlip = config_phase(p)%getFloat('h0_twinslip') else slipAndTwinActive + allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension 0 + allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension 0 prm%h0_TwinSlip = 0.0_pReal endif slipAndTwinActive From 35956ce0c83fa0d7bd293fa135663c1a045ac071 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 4 Oct 2018 21:04:19 +0200 Subject: [PATCH 25/45] needs to be multiplikative --- src/plastic_phenopowerlaw.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index e9663ad4e..792889c13 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -62,7 +62,7 @@ module plastic_phenopowerlaw xi_slip_sat, & !< maximum critical shear stress for slip nonSchmidCoeff, & H_int, & !< per family hardening activity (optional) !ToDo: Better name! - gamma_twin_char !< characteristic shear for twins + gamma_twin_char !< characteristic shear for twins real(pReal), dimension(:,:), allocatable :: & interaction_SlipSlip, & !< slip resistance from slip activity interaction_SlipTwin, & !< slip resistance from twin activity @@ -504,8 +504,8 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) Lp = Lp + (1.0_pReal-stt%sumF(of))*(gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + dgdot_dtauslip_pos(i) * (prm%Schmid_slip(k,l,i) + prm%nonSchmid_pos(m,n,i)) & - + dgdot_dtauslip_neg(i) * (prm%Schmid_slip(k,l,i) + prm%nonSchmid_neg(m,n,i)) + + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & + + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems call kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtautwin) From 79b0efc678348d31a291358e9d473beab1ce613b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 4 Oct 2018 21:18:13 +0200 Subject: [PATCH 26/45] also using nonSchmid matrix from lattice --- src/lattice.f90 | 14 ++++++++------ src/plastic_phenopowerlaw.f90 | 28 +++------------------------- 2 files changed, 11 insertions(+), 31 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 090dad385..46f243561 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1229,6 +1229,7 @@ real(pReal), dimension(4,36), parameter, private :: & LATTICE_hex_ID, & lattice_SchmidMatrix_slip, & lattice_SchmidMatrix_twin, & + lattice_nonSchmidMatrix, & lattice_interactionSlipSlip2, & lattice_interactionTwinTwin2, & lattice_interactionSlipTwin2, & @@ -2226,13 +2227,14 @@ function lattice_C66_trans(Ntrans,C66_parent,structure_parent,cOverA_parent, & ! ! Schmid matrices with non-Schmid contributions according to Koester_etal2012, Acta Materialia 60 (2012) 3894–3901, eq. (17) ("n1" is replaced by either "np" or "nn" according to either positive or negative slip direction) ! ! "np" and "nn" according to Gröger_etal2008, Acta Materialia 56 (2008) 5412–5425, table 1 (corresponds to their "n1" for positive and negative slip direction respectively) + function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix) use math, only: & - INRAD, & - math_tensorproduct33, & - math_crossproduct, & - math_mul33x3, & - math_axisAngleToR + INRAD, & + math_tensorproduct33, & + math_crossproduct, & + math_mul33x3, & + math_axisAngleToR implicit none integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family real(pReal), dimension(6), intent(in) :: nonSchmidCoefficients @@ -2253,7 +2255,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc do i = 1_pInt,sum(Nslip) direction = coordinateSystem(1:3,1,i) normal = coordinateSystem(1:3,2,i) - np = math_mul33x3(math_axisAngleToR(+direction,60.0_pReal*INRAD), normal) + np = math_mul33x3(math_axisAngleToR(direction,60.0_pReal*INRAD), normal) nonSchmidMatrix(1:3,1:3,i) & = nonSchmidMatrix(1:3,1:3,i) & + nonSchmidCoefficients(1) * math_tensorproduct33(direction, np) & diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 792889c13..a53f0103f 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -213,7 +213,8 @@ subroutine plastic_phenopowerlaw_init defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) prm%nonSchmidCoeff = config_phase(p)%getFloats('nonschmid_coefficients',& defaultVal = emptyRealArray ) - + prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) prm%gdot0_slip = config_phase(p)%getFloat('gdot0_slip') prm%n_slip = config_phase(p)%getFloat('n_slip') prm%a_slip = config_phase(p)%getFloat('a_slip') @@ -381,32 +382,9 @@ subroutine plastic_phenopowerlaw_init allocate(plasticState(p)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal) -!-------------------------------------------------------------------------------------------------- -! calculate hardening matrices - allocate(prm%nonSchmid_pos(3,3,prm%totalNslip),source = 0.0_pReal) - allocate(prm%nonSchmid_neg(3,3,prm%totalNslip),source = 0.0_pReal) - i = 0_pInt - mySlipFamilies: do f = 1_pInt,size(prm%Nslip,1) ! >>> interaction slip -- X - index_myFamily = sum(prm%Nslip(1:f-1_pInt)) - - mySlipSystems: do j = 1_pInt,prm%Nslip(f) - i = i + 1_pInt - prm%nonSchmid_pos(1:3,1:3,i) = lattice_Sslip(1:3,1:3,1, index_myFamily+j,p) - prm%nonSchmid_neg(1:3,1:3,i) = lattice_Sslip(1:3,1:3,1, index_myFamily+j,p) - do k = 1,size(prm%nonSchmidCoeff) - prm%nonSchmid_pos(1:3,1:3,i) = prm%nonSchmid_pos(1:3,1:3,i) & - + lattice_Sslip(1:3,1:3,2*k, index_myFamily+j,p) & - * prm%nonSchmidCoeff(k) - prm%nonSchmid_neg(1:3,1:3,i) = prm%nonSchmid_neg(1:3,1:3,i) & - + lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+j,p) & - * prm%nonSchmidCoeff(k) - enddo - enddo mySlipSystems - enddo mySlipFamilies - allocate(prm%gamma_twin_char(prm%totalNtwin),source = 0.0_pReal) i = 0_pInt - myTwinFamilies: do f = 1_pInt,size(prm%Ntwin,1) ! >>> interaction twin -- X + myTwinFamilies: do f = 1_pInt,size(prm%Ntwin,1) index_myFamily = sum(prm%Ntwin(1:f-1_pInt)) myTwinSystems: do j = 1_pInt,prm%Ntwin(f) i = i + 1_pInt From fda267636545da389cf4473e7d687e7dda3ad016 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 5 Oct 2018 04:54:47 +0200 Subject: [PATCH 27/45] calculate only if required --- src/lattice.f90 | 46 +++++++++++++++++++++++++++++++++++ src/plastic_phenopowerlaw.f90 | 2 +- 2 files changed, 47 insertions(+), 1 deletion(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 46f243561..573e67a7f 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2084,6 +2084,52 @@ pure function lattice_qDisorientation(Q1, Q2, struct) end function lattice_qDisorientation +function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(characteristicShear) + use IO, only: & + IO_error + + implicit none + integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family + character(len=3), intent(in) :: structure + real(pReal), intent(in) :: cOverA + real(pReal), dimension(sum(Ntwin)) :: characteristicShear + integer(pInt) :: & + ir, & !< index in reduced list + ig, & !< index in full list + mf, & !< index of my family + ms !< index of my system in current family + + select case(structure) + case('fcc') + characteristicShear = LATTICE_FCC_SHEARTWIN + case('bcc') + characteristicShear = LATTICE_BCC_SHEARTWIN + case('hex') + ir = 0_pInt + myFamilies: do mf = 1_pInt,size(Ntwin,1) + mySystems: do ms = 1_pInt,Ntwin(mf) + ir = ir + 1_pInt + ig = sum(LATTICE_HEX_NTWINSYSTEM(1:mf-1))+ms + characteristicShear = LATTICE_BCC_SHEARTWIN + select case(LATTICE_HEX_SHEARTWIN(ig)) ! from Christian & Mahajan 1995 p.29 + case (1_pInt) ! <-10.1>{10.2} + characteristicShear(ir) = (3.0_pReal-cOverA*cOverA)/sqrt(3.0_pReal)/CoverA + case (2_pInt) ! <11.6>{-1-1.1} + characteristicShear(ir) = 1.0_pReal/cOverA + case (3_pInt) ! <10.-2>{10.1} + characteristicShear(ir) = (4.0_pReal*cOverA*cOverA-9.0_pReal)/4.0_pReal & + / sqrt(3.0_pReal)/cOverA + !characteristicShear(ir) = (4.0_pReal*cOverA*cOverA-9.0_pReal)/sqrt(48.0_pReal)/cOverA + case (4_pInt) ! <11.-3>{11.2} + characteristicShear(ir) = 2.0_pReal*(cOverA*cOverA-2.0_pReal)/3.0_pReal/cOverA + end select + enddo mySystems + enddo myFamilies + end select +end function lattice_characteristicShear_Twin + + + function lattice_C66_twin(Ntwin,C66,structure,CoverA) use IO, only: & diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index a53f0103f..a209a1da6 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -162,7 +162,7 @@ subroutine plastic_phenopowerlaw_init dot integer(kind(undefined_ID)) :: & - outputID !< ID of each post result output + outputID !< ID of each post result output character(len=512) :: & extmsg = '', & From 3a00fcde16bbac37f86ca4686805adc654b2130d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 5 Oct 2018 05:09:31 +0200 Subject: [PATCH 28/45] WIP: characteristic shear calculated in lattice --- src/lattice.f90 | 3 ++- src/plastic_phenopowerlaw.f90 | 12 +++--------- 2 files changed, 5 insertions(+), 10 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 573e67a7f..15fd71b8a 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1233,7 +1233,8 @@ real(pReal), dimension(4,36), parameter, private :: & lattice_interactionSlipSlip2, & lattice_interactionTwinTwin2, & lattice_interactionSlipTwin2, & - lattice_interactionTwinSlip2 + lattice_interactionTwinSlip2, & + lattice_characteristicShear_Twin contains diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index a209a1da6..1532b9d52 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -272,6 +272,9 @@ subroutine plastic_phenopowerlaw_init allocate(prm%xi_twin_0(0)) endif twinActive + prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),& + config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) + slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then prm%interaction_SlipTwin = lattice_interactionSlipTwin2(prm%Nslip,prm%Ntwin,& config_phase(p)%getFloats('interaction_sliptwin'), & @@ -382,15 +385,6 @@ subroutine plastic_phenopowerlaw_init allocate(plasticState(p)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal) - allocate(prm%gamma_twin_char(prm%totalNtwin),source = 0.0_pReal) - i = 0_pInt - myTwinFamilies: do f = 1_pInt,size(prm%Ntwin,1) - index_myFamily = sum(prm%Ntwin(1:f-1_pInt)) - myTwinSystems: do j = 1_pInt,prm%Ntwin(f) - i = i + 1_pInt - prm%gamma_twin_char(i) = lattice_shearTwin(sum(lattice_Ntwinsystem(1:f-1,p))+j,p) - enddo myTwinSystems - enddo myTwinFamilies !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState From 6ccbc5e4e66999c7d1dcf14ad8bc121764604359 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 5 Oct 2018 07:20:51 +0200 Subject: [PATCH 29/45] characteristing shear is defined per system --- src/lattice.f90 | 32 ++++++++++++++++---------------- src/plastic_phenopowerlaw.f90 | 4 +--- 2 files changed, 17 insertions(+), 19 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 15fd71b8a..3ccd65fea 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2099,19 +2099,18 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact ig, & !< index in full list mf, & !< index of my family ms !< index of my system in current family - - select case(structure) - case('fcc') - characteristicShear = LATTICE_FCC_SHEARTWIN - case('bcc') - characteristicShear = LATTICE_BCC_SHEARTWIN - case('hex') - ir = 0_pInt - myFamilies: do mf = 1_pInt,size(Ntwin,1) - mySystems: do ms = 1_pInt,Ntwin(mf) - ir = ir + 1_pInt - ig = sum(LATTICE_HEX_NTWINSYSTEM(1:mf-1))+ms - characteristicShear = LATTICE_BCC_SHEARTWIN + return + ir = 0_pInt + myFamilies: do mf = 1_pInt,size(Ntwin,1) + mySystems: do ms = 1_pInt,Ntwin(mf) + ir = ir + 1_pInt + ig = sum(LATTICE_HEX_NTWINSYSTEM(1:mf-1))+ms + select case(structure) + case('fcc') + characteristicShear(ir) = LATTICE_FCC_SHEARTWIN(ig) + case('bcc') + characteristicShear(ir) = LATTICE_BCC_SHEARTWIN(ig) + case('hex') select case(LATTICE_HEX_SHEARTWIN(ig)) ! from Christian & Mahajan 1995 p.29 case (1_pInt) ! <-10.1>{10.2} characteristicShear(ir) = (3.0_pReal-cOverA*cOverA)/sqrt(3.0_pReal)/CoverA @@ -2124,9 +2123,10 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact case (4_pInt) ! <11.-3>{11.2} characteristicShear(ir) = 2.0_pReal*(cOverA*cOverA-2.0_pReal)/3.0_pReal/cOverA end select - enddo mySystems - enddo myFamilies - end select + end select + enddo mySystems + enddo myFamilies + end function lattice_characteristicShear_Twin diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 1532b9d52..b957e76f0 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -145,9 +145,8 @@ subroutine plastic_phenopowerlaw_init integer(pInt) :: & maxNinstance, & - instance,p,j,k, f,o, i,& + instance,p,j,k, o, i,& NipcMyPhase, outputSize, & - index_myFamily, index_otherFamily, & sizeState,sizeDotState, & startIndex, endIndex @@ -385,7 +384,6 @@ subroutine plastic_phenopowerlaw_init allocate(plasticState(p)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal) - !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1_pInt From aa34b83978994f0e9825706e3449a552653cfa93 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 5 Oct 2018 10:36:44 +0200 Subject: [PATCH 30/45] avoid array out of bounds allows to handle un- or partly defined non schmid coefficients --- src/lattice.f90 | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 46f243561..a963d947d 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2237,7 +2237,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc math_axisAngleToR implicit none integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family - real(pReal), dimension(6), intent(in) :: nonSchmidCoefficients + real(pReal), dimension(:), intent(in) :: nonSchmidCoefficients integer(pInt), intent(in) :: sense !< sense (-1,+1) real(pReal), dimension(1:3,1:3,sum(Nslip)) :: nonSchmidMatrix @@ -2256,15 +2256,19 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc direction = coordinateSystem(1:3,1,i) normal = coordinateSystem(1:3,2,i) np = math_mul33x3(math_axisAngleToR(direction,60.0_pReal*INRAD), normal) - nonSchmidMatrix(1:3,1:3,i) & - = nonSchmidMatrix(1:3,1:3,i) & - + nonSchmidCoefficients(1) * math_tensorproduct33(direction, np) & - + nonSchmidCoefficients(2) * math_tensorproduct33(math_crossproduct(normal, direction), normal) & - + nonSchmidCoefficients(3) * math_tensorproduct33(math_crossproduct(np, direction), np) & - + nonSchmidCoefficients(4) * math_tensorproduct33(normal, normal) & - + nonSchmidCoefficients(5) * math_tensorproduct33(math_crossproduct(normal, direction), & - math_crossproduct(normal, normal)) & - + nonSchmidCoefficients(6) * math_tensorproduct33(direction, direction) + if (size(nonSchmidCoefficients)>1) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + + nonSchmidCoefficients(1) * math_tensorproduct33(direction, np) + if (size(nonSchmidCoefficients)>2) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + + nonSchmidCoefficients(2) * math_tensorproduct33(math_crossproduct(normal, direction), normal) + if (size(nonSchmidCoefficients)>3) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + + nonSchmidCoefficients(3) * math_tensorproduct33(math_crossproduct(np, direction), np) + if (size(nonSchmidCoefficients)>4) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + + nonSchmidCoefficients(4) * math_tensorproduct33(normal, normal) + if (size(nonSchmidCoefficients)>5) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + + nonSchmidCoefficients(5) * math_tensorproduct33(math_crossproduct(normal, direction), & + math_crossproduct(normal, normal)) + if (size(nonSchmidCoefficients)>6) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + + nonSchmidCoefficients(6) * math_tensorproduct33(direction, direction) enddo end function lattice_nonSchmidMatrix From 81a77b508a991e360494da4bfc36e50fd33adc5f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 5 Oct 2018 22:46:53 +0200 Subject: [PATCH 31/45] check for array size was off by 1 --- src/lattice.f90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 25d34f5a7..9b0d39db6 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2099,7 +2099,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact ig, & !< index in full list mf, & !< index of my family ms !< index of my system in current family - return + ir = 0_pInt myFamilies: do mf = 1_pInt,size(Ntwin,1) mySystems: do ms = 1_pInt,Ntwin(mf) @@ -2303,18 +2303,18 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc direction = coordinateSystem(1:3,1,i) normal = coordinateSystem(1:3,2,i) np = math_mul33x3(math_axisAngleToR(direction,60.0_pReal*INRAD), normal) - if (size(nonSchmidCoefficients)>1) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + if (size(nonSchmidCoefficients)>0) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + nonSchmidCoefficients(1) * math_tensorproduct33(direction, np) - if (size(nonSchmidCoefficients)>2) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + if (size(nonSchmidCoefficients)>1) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + nonSchmidCoefficients(2) * math_tensorproduct33(math_crossproduct(normal, direction), normal) - if (size(nonSchmidCoefficients)>3) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + if (size(nonSchmidCoefficients)>2) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + nonSchmidCoefficients(3) * math_tensorproduct33(math_crossproduct(np, direction), np) - if (size(nonSchmidCoefficients)>4) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + if (size(nonSchmidCoefficients)>3) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + nonSchmidCoefficients(4) * math_tensorproduct33(normal, normal) - if (size(nonSchmidCoefficients)>5) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + if (size(nonSchmidCoefficients)>4) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + nonSchmidCoefficients(5) * math_tensorproduct33(math_crossproduct(normal, direction), & math_crossproduct(normal, normal)) - if (size(nonSchmidCoefficients)>6) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + if (size(nonSchmidCoefficients)>5) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + nonSchmidCoefficients(6) * math_tensorproduct33(direction, direction) enddo end function lattice_nonSchmidMatrix From 3a2f86df1c5198677f257332e8af5790eb96a96a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 5 Oct 2018 23:04:06 +0200 Subject: [PATCH 32/45] non-schmid is always bcc --- src/plastic_phenopowerlaw.f90 | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index f695fdf67..1ca5cd66c 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -212,8 +212,13 @@ subroutine plastic_phenopowerlaw_init defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) prm%nonSchmidCoeff = config_phase(p)%getFloats('nonschmid_coefficients',& defaultVal = emptyRealArray ) - prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) + if(structure=='bcc') then + prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) + else + prm%nonSchmid_pos = prm%Schmid_slip + prm%nonSchmid_neg = prm%Schmid_slip + endif prm%gdot0_slip = config_phase(p)%getFloat('gdot0_slip') prm%n_slip = config_phase(p)%getFloat('n_slip') prm%a_slip = config_phase(p)%getFloat('a_slip') From 032c35a4998162f310061f74715d779188f36613 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 6 Oct 2018 10:42:25 +0200 Subject: [PATCH 33/45] WIP: Calculation of elasticity matrices for twin and trans --- src/lattice.f90 | 197 +++++++++++++++++++++++++++--------------------- 1 file changed, 111 insertions(+), 86 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 9b0d39db6..901cbda2d 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2085,6 +2085,10 @@ pure function lattice_qDisorientation(Q1, Q2, struct) end function lattice_qDisorientation + +!-------------------------------------------------------------------------------------------------- +!> @brief Provides characteristtic shear for twinning +!-------------------------------------------------------------------------------------------------- function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(characteristicShear) use IO, only: & IO_error @@ -2107,10 +2111,13 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact ig = sum(LATTICE_HEX_NTWINSYSTEM(1:mf-1))+ms select case(structure) case('fcc') + ig = sum(LATTICE_FCC_NTWINSYSTEM(1:mf-1))+ms characteristicShear(ir) = LATTICE_FCC_SHEARTWIN(ig) case('bcc') + ig = sum(LATTICE_BCC_NTWINSYSTEM(1:mf-1))+ms characteristicShear(ir) = LATTICE_BCC_SHEARTWIN(ig) case('hex') + ig = sum(LATTICE_HEX_NTWINSYSTEM(1:mf-1))+ms select case(LATTICE_HEX_SHEARTWIN(ig)) ! from Christian & Mahajan 1995 p.29 case (1_pInt) ! <-10.1>{10.2} characteristicShear(ir) = (3.0_pReal-cOverA*cOverA)/sqrt(3.0_pReal)/CoverA @@ -2130,8 +2137,9 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact end function lattice_characteristicShear_Twin - - +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates rotated elasticity matrices for twinning +!-------------------------------------------------------------------------------------------------- function lattice_C66_twin(Ntwin,C66,structure,CoverA) use IO, only: & IO_error @@ -2172,109 +2180,126 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) end function -function lattice_C66_trans(Ntrans,C66_parent,structure_parent,cOverA_parent, & - C66_target,structure_target,cOverA_target) +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates rotated elasticity matrices for transformation +!-------------------------------------------------------------------------------------------------- +function lattice_C66_trans(Ntrans,C_parent66,structure_parent, & + C_target66,structure_target) + use prec, only: & + tol_math_check use IO, only: & IO_error use math, only: & INRAD, & + MATH_I3, & math_axisAngleToR, & math_Mandel3333to66, & math_Mandel66to3333, & - math_rotate_forward3333 + math_rotate_forward3333, & + math_mul33x33, & + math_tensorproduct33, & + math_crossproduct implicit none integer(pInt), dimension(:), intent(in) :: Ntrans !< number of active twin systems per family character(len=*), intent(in) :: & structure_target, & !< lattice structure structure_parent !< lattice structure - real(pReal), dimension(6,6), intent(in) :: C66_target, C66_parent - real(pReal), intent(in) :: cOverA_parent, cOverA_target - real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans, C_bar66, C_target_unrotated66 + real(pReal), dimension(6,6), intent(in) :: C_parent66, C_target66 + real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66 + real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans - real(pReal), dimension(3,3,sum(Ntrans)) :: coordinateSystem - - real(pReal), dimension(3,3) :: R - real(pReal) :: A, B + real(pReal), dimension(3,3) :: R,B,U,Q,S,ss,sd,st + real(pReal), dimension(3) :: x,y,z + real(pReal) :: a_bcc, a_fcc, CoverA_trans integer(pInt) :: i -! if (trim(structure_parent) == 'fcc' .and. trim(structure_target) == 'hex') then -! C_bar66(1,1) = (C66_parent(1,1) + C66_parent(1,2) + 2.0_pReal*C66_parent(4,4))/2.0_pReal -! C_bar66(1,2) = (C66_parent(1,1) + 5.0_pReal*C66_parent(1,2) - 2.0_pReal*C66_parent(4,4))/6.0_pReal -! C_bar66(3,3) = (C66_parent(1,1) + 2.0_pReal*C66_parent(1,2) + 4.0_pReal*C66_parent(4,4))/3.0_pReal -! C_bar66(1,3) = (C66_parent(1,1) + 2.0_pReal*C66_parent(1,2) - 2.0_pReal*C66_parent(4,4))/3.0_pReal -! C_bar66(4,4) = (C66_parent(1,1) - C66_parent(1,2) + C66_parent(4,4))/3.0_pReal -! C_bar66(1,4) = (C66_parent(1,1) - C66_parent(1,2) - 2.0_pReal*C66_parent(4,4)) /(3.0_pReal*sqrt(2.0_pReal)) -! A = C_bar(1,4)**2.0_pReal/C_bar66(4,4) -! B = C_bar(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2))) -! C_target_unrotated66(1,1) = C_bar66(1,1) - A -! C_target_unrotated66(1,2) = C_bar66(1,2) + A -! C_target_unrotated66(1,3) = C_bar66(1,3) -! C_target_unrotated66(3,3) = C_bar66(3,3) -! C_target_unrotated66(4,4) = C_bar66(4,4) - B -! C_target_unrotated66 = lattice_symmetrizeC66(LATTICE_HEX_ID,C_target_unrotated) -! elseif (trim(structure_parent) == 'fcc' .and. trim(structure_target) == 'bcc') then -! C_target_unrotated66 = C_parent66 -! endif -! do i = 1_pInt, 6_pInt -! if (abs(C_target_unrotated66(i,i)) 0.0_pReal) .and. (a_bcc > 0.0_pReal)) then + U = (a_bcc/a_fcc)*math_tensorproduct33(x,x) + & + (a_bcc/a_fcc)*math_tensorproduct33(y,y) * sqrt(2.0_pReal) + & + (a_bcc/a_fcc)*math_tensorproduct33(z,z) * sqrt(2.0_pReal) + else BainDeformation + U = 0.0_pReal + endif BainDeformation + Q = math_mul33x33(R,B) + S = math_mul33x33(R,U) - MATH_I3 + enddo + elseif (trim(structure_target) == 'bcc') then + ss = MATH_I3 + ss(1,3) = sqrt(0.125_pReal) + sd = MATH_I3 + if (CoverA_trans > 1.0_pReal .and. CoverA_trans < 2.0_pReal) then + sd(3,3) = CoverA_trans/sqrt(8.0_pReal/3.0_pReal) + endif + st = math_mul33x33(sd,ss) + do i = 1_pInt,sum(Ntrans)!!!!!!!!!!!!!! NEED TO BE FIXED + R(1:3,1) = lattice_fccTohex_systemTrans(1:3,i)/norm2(lattice_fccTohex_systemTrans(1:3,i)) + R(1:3,3) = lattice_fccTohex_systemTrans(4:6,i)/norm2(lattice_fccTohex_systemTrans(4:6,i)) + R(1:3,2) = -math_crossproduct(R(1:3,1),R(1:3,3)) + Q = R + S = math_mul33x33(R, math_mul33x33(st, transpose(R))) - MATH_I3 + ! trs(i) = lattice_fccTohex_shearTrans(i) + enddo + else + write(6,*) "Mist" + endif + + + do i = 1, sum(Ntrans) ! R = math_axisAngleToR(coordinateSystem(1:3,2,i), 180.0_pReal * INRAD) ! ToDo: Why always 180 deg? ! lattice_C66_trans(1:6,1:6,i) = math_Mandel3333to66(math_rotate_forward3333(math_Mandel66to3333(C66),R)) -! enddo -! ! Phase transformation -! select case(trans_lattice_structure(myPhase)) -! case (LATTICE_bcc_ID) ! fcc to bcc transformation -! do i = 1_pInt,myNtrans -! Rtr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), & ! Pitsch rotation -! lattice_fccTobcc_systemTrans(4,i)*INRAD) -! Btr(1:3,1:3,i) = math_axisAngleToR(lattice_fccTobcc_bainRot(1:3,i), & ! Rotation of fcc to Bain coordinate system -! lattice_fccTobcc_bainRot(4,i)*INRAD) -! xtr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(1:3,i),pReal) -! ytr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(4:6,i),pReal) -! ztr(1:3,i) = real(LATTICE_fccTobcc_bainVariant(7:9,i),pReal) -! Utr(1:3,1:3,i) = 0.0_pReal ! Bain deformation -! if ((a_fcc > 0.0_pReal) .and. (a_bcc > 0.0_pReal)) then -! Utr(1:3,1:3,i) = (a_bcc/a_fcc)*math_tensorproduct33(xtr(1:3,i), xtr(1:3,i)) + & -! sqrt(2.0_pReal)*(a_bcc/a_fcc)*math_tensorproduct33(ytr(1:3,i), ytr(1:3,i)) + & -! sqrt(2.0_pReal)*(a_bcc/a_fcc)*math_tensorproduct33(ztr(1:3,i), ztr(1:3,i)) -! endif -! Qtr(1:3,1:3,i) = math_mul33x33(Rtr(1:3,1:3,i), Btr(1:3,1:3,i)) -! Str(1:3,1:3,i) = math_mul33x33(Rtr(1:3,1:3,i), Utr(1:3,1:3,i)) - MATH_I3 -! enddo -! case (LATTICE_hex_ID) -! sstr(1:3,1:3) = MATH_I3 -! sstr(1,3) = sqrt(2.0_pReal)/4.0_pReal -! sdtr(1:3,1:3) = MATH_I3 -! if (CoverA_trans > 1.0_pReal .and. CoverA_trans < 2.0_pReal) then -! sdtr(3,3) = CoverA_trans/sqrt(8.0_pReal/3.0_pReal) -! endif -! sttr = math_mul33x33(sdtr, sstr) -! do i = 1_pInt,myNtrans -! xtr(1:3,i) = lattice_fccTohex_systemTrans(1:3,i)/norm2(lattice_fccTohex_systemTrans(1:3,i)) -! ztr(1:3,i) = lattice_fccTohex_systemTrans(4:6,i)/norm2(lattice_fccTohex_systemTrans(4:6,i)) -! ytr(1:3,i) = -math_crossproduct(xtr(1:3,i), ztr(1:3,i)) -! Rtr(1:3,1,i) = xtr(1:3,i) -! Rtr(1:3,2,i) = ytr(1:3,i) -! Rtr(1:3,3,i) = ztr(1:3,i) -! Qtr(1:3,1:3,i) = Rtr(1:3,1:3,i) -! Str(1:3,1:3,i) = math_mul33x33(Rtr(1:3,1:3,i), math_mul33x33(sttr, math_transpose33(Rtr(1:3,1:3,i)))) -! Str(1:3,1:3,i) = Str(1:3,1:3,i) - MATH_I3 -! trs(i) = lattice_fccTohex_shearTrans(i) -! enddo -! case default -! Qtr = 0.0_pReal -! Str = 0.0_pReal -! end select + enddo +end function - end function -! ! Schmid matrices with non-Schmid contributions according to Koester_etal2012, Acta Materialia 60 (2012) 3894–3901, eq. (17) ("n1" is replaced by either "np" or "nn" according to either positive or negative slip direction) -! ! "np" and "nn" according to Gröger_etal2008, Acta Materialia 56 (2008) 5412–5425, table 1 (corresponds to their "n1" for positive and negative slip direction respectively) +!-------------------------------------------------------------------------------------------------- +!> @brief Non-schmid tensor +! Schmid matrices with non-Schmid contributions according to Koester_etal2012, Acta Materialia 60 (2012) +! 3894–3901, eq. (17) ("n1" is replaced by either "np" or "nn" according to either positive or negative slip direction) +! "np" and "nn" according to Gröger_etal2008, Acta Materialia 56 (2008) 5412–5425, table 1 +! (corresponds to their "n1" for positive and negative slip direction respectively) +!-------------------------------------------------------------------------------------------------- function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix) use math, only: & INRAD, & @@ -2700,8 +2725,8 @@ pure function buildCoordinateSystem(active,system,structure,cOverA) end select - buildCoordinateSystem(1:3,1,ir) = direction/norm2(direction) ! make unit vector - buildCoordinateSystem(1:3,2,ir) = normal/norm2(normal) ! make unit vector + buildCoordinateSystem(1:3,1,ir) = direction/norm2(direction) + buildCoordinateSystem(1:3,2,ir) = normal/norm2(normal) buildCoordinateSystem(1:3,3,ir) = math_crossproduct(direction,normal) enddo mySystems From a53488d666e614fecdc0af5259187f2b80ceae1e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 7 Oct 2018 15:44:13 +0200 Subject: [PATCH 34/45] calculation of slip/twin/trans/damage-coordinate system was wrong --- src/lattice.f90 | 95 +++++++++++++++++++++++++------------------------ 1 file changed, 49 insertions(+), 46 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 901cbda2d..724857b6d 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2164,11 +2164,11 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) select case(structure) case('fcc') - coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) case('bcc') - coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) case('hex','hexagonal') !ToDo: "No alias policy": long or short? - coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NSLIPSYSTEM,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) case default call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_C66_twin)') end select @@ -2320,8 +2320,8 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc integer(pInt) :: i if (abs(sense) /= 1_pInt) write(6,*) 'mist' - coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMSLIP,pInt),'bcc') - coordinateSystem(1:3,1,sum(Nslip)) = coordinateSystem(1:3,1,sum(Nslip)) *real(sense,pReal) + coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,int(LATTICE_BCC_SYSTEMSLIP,pInt),'bcc') + coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip)) *real(sense,pReal) nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'bcc') do i = 1_pInt,sum(Nslip) @@ -2577,13 +2577,13 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) select case(structure) case('fcc') - coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_FCC_SYSTEMSLIP,pInt),structure) + coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,int(LATTICE_FCC_SYSTEMSLIP,pInt),structure) case('bcc') - coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCC_SYSTEMSLIP,pInt),structure) + coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_FCC_NSLIPSYSTEM,int(LATTICE_BCC_SYSTEMSLIP,pInt),structure) case('hex','hexagonal') !ToDo: "No alias policy": long or short? - coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_HEX_SYSTEMSLIP,pInt),'hex',cOverA) + coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_HEX_NSLIPSYSTEM,int(LATTICE_HEX_SYSTEMSLIP,pInt),'hex',cOverA) case('bct') - coordinateSystem = buildCoordinateSystem(Nslip,int(LATTICE_BCT_SYSTEMSLIP,pInt),structure,cOverA) + coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCT_NSLIPSYSTEM,int(LATTICE_BCT_SYSTEMSLIP,pInt),structure,cOverA) case default call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_slip)') end select @@ -2618,11 +2618,11 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) select case(structure) case('fcc') - coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NTWINSYSTEM,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) case('bcc') - coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NTWINSYSTEM,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) case('hex','hexagonal') !ToDo: "No alias policy": long or short? - coordinateSystem = buildCoordinateSystem(Ntwin,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NTWINSYSTEM,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) case default call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_twin)') end select @@ -2655,14 +2655,16 @@ pure function buildInteraction(activeA,activeB,maxA,maxB,values,matrix) mf, ms, of, os myFamilies: do mf = 1_pInt,size(activeA,1) - index_myFamily = sum(activeA(1:mf-1_pInt)) - mySystems: do ms = 1_pInt,activeA(mf) - otherFamilies: do of = 1_pInt,size(activeB,1) - index_otherFamily = sum(activeB(1:of-1_pInt)) - otherSystems: do os = 1_pInt,activeB(of) - buildInteraction(index_myFamily+ms,index_otherFamily+os) = & - values(matrix(sum(maxA(1:mf-1))+ms, sum(maxB(1:of-1))+os)) - enddo otherSystems; enddo otherFamilies; + index_myFamily = sum(activeA(1:mf-1_pInt)) + mySystems: do ms = 1_pInt,activeA(mf) + + otherFamilies: do of = 1_pInt,size(activeB,1) + index_otherFamily = sum(activeB(1:of-1_pInt)) + otherSystems: do os = 1_pInt,activeB(of) + buildInteraction(index_myFamily+ms,index_otherFamily+os) = & + values(matrix(sum(maxA(1:mf-1))+ms, sum(maxB(1:of-1))+os)) + enddo otherSystems; enddo otherFamilies; + enddo mySystems;enddo myFamilies end function buildInteraction @@ -2672,13 +2674,14 @@ end function buildInteraction !> @brief build a local coordinate system in a slip, twin, trans, cleavage system !> @details: Order: Direction, plane (normal), and common perpendicular !-------------------------------------------------------------------------------------------------- -pure function buildCoordinateSystem(active,system,structure,cOverA) +function buildCoordinateSystem(active,maximum,system,structure,cOverA) use math, only: & math_crossproduct implicit none integer(pInt), dimension(:), intent(in) :: & - active + active, & + maximum integer(pInt), dimension(:,:), intent(in) :: & system character(len=*), intent(in) :: & @@ -2691,46 +2694,46 @@ pure function buildCoordinateSystem(active,system,structure,cOverA) real(pReal), dimension(3) :: & direction, normal integer(pInt) :: & - ir, & !< index in reduced matrix - ig, & !< index in full matrix - mf, & !< index of my family - ms !< index of my system in current family + i, & !< index in reduced matrix + j, & !< index in full matrix + f, & !< index of my family + s !< index of my system in current family - ir = 0_pInt - myFamilies: do mf = 1_pInt,size(active,1) - mySystems: do ms = 1_pInt,active(mf) - ir = ir + 1_pInt - ig = sum(system(:,1:mf-1))+ms + i = 0_pInt + activeFamilies: do f = 1_pInt,size(active,1) + activeSystems: do s = 1_pInt,active(f) + i = i + 1_pInt + j = sum(maximum(1:f-1))+s select case(trim(structure)) case ('fcc','bcc') - direction = real(system(1:3,ig),pReal) - normal = real(system(4:6,ig),pReal) + direction = real(system(1:3,j),pReal) + normal = real(system(4:6,j),pReal) case ('hex') ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)]) - direction = [ real(system(1,ig),pReal)*1.5_pReal, & - (real(system(1,ig),pReal)+2.0_pReal*real(system(2,ig),pReal))*sqrt(0.75_pReal), & - real(system(4,ig),pReal)*CoverA ] + direction = [ real(system(1,j),pReal)*1.5_pReal, & + (real(system(1,j),pReal)+2.0_pReal*real(system(2,j),pReal))*sqrt(0.75_pReal), & + real(system(4,j),pReal)*CoverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) - normal = [ real(system(5,ig),pReal), & - (real(system(5,ig),pReal)+2.0_pReal*real(system(6,ig),pReal))/ sqrt(3.0_pReal), & - real(system(8,ig),pReal)/CoverA ] + normal = [ real(system(5,j),pReal), & + (real(system(5,j),pReal)+2.0_pReal*real(system(6,j),pReal))/ sqrt(3.0_pReal), & + real(system(8,j),pReal)/CoverA ] case ('bct') - direction = [real(system(1:2,ig),pReal),real(system(3,ig),pReal)*CoverA] - normal = [real(system(4:5,ig),pReal),real(system(6,ig),pReal)/CoverA] + direction = [real(system(1:2,j),pReal),real(system(3,i),pReal)*CoverA] + normal = [real(system(4:5,j),pReal),real(system(6,i),pReal)/CoverA] end select - buildCoordinateSystem(1:3,1,ir) = direction/norm2(direction) - buildCoordinateSystem(1:3,2,ir) = normal/norm2(normal) - buildCoordinateSystem(1:3,3,ir) = math_crossproduct(direction,normal) + buildCoordinateSystem(1:3,1,i) = direction/norm2(direction) + buildCoordinateSystem(1:3,2,i) = normal/norm2(normal) + buildCoordinateSystem(1:3,3,i) = math_crossproduct(direction,normal) - enddo mySystems - enddo myFamilies + enddo activeSystems + enddo activeFamilies end function buildCoordinateSystem From 34e0aca564f39e2b5203a9f02581533983181b79 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 7 Oct 2018 18:18:24 +0200 Subject: [PATCH 35/45] using full precision --- src/lattice.f90 | 72 ++++++++++++++++++++++++------------------------- 1 file changed, 36 insertions(+), 36 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 724857b6d..a3a790574 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1181,43 +1181,43 @@ module lattice real(pReal), dimension(4,36), parameter, private :: & lattice_symOperations = reshape([& - 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! cubic symmetry operations - 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal, & ! 2-fold symmetry - 0.0_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, & - 0.0_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, -0.7071067811865476_pReal, & - 0.0_pReal, -0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, & - 0.0_pReal, 0.7071067811865476_pReal, -0.7071067811865476_pReal, 0.0_pReal, & - 0.5_pReal, 0.5_pReal, 0.5_pReal, 0.5_pReal, & ! 3-fold symmetry - -0.5_pReal, 0.5_pReal, 0.5_pReal, 0.5_pReal, & - 0.5_pReal, -0.5_pReal, 0.5_pReal, 0.5_pReal, & - -0.5_pReal, -0.5_pReal, 0.5_pReal, 0.5_pReal, & - 0.5_pReal, 0.5_pReal, -0.5_pReal, 0.5_pReal, & - -0.5_pReal, 0.5_pReal, -0.5_pReal, 0.5_pReal, & - 0.5_pReal, 0.5_pReal, 0.5_pReal, -0.5_pReal, & - -0.5_pReal, 0.5_pReal, 0.5_pReal, -0.5_pReal, & - 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, & ! 4-fold symmetry - 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, & - -0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, & - 0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, & - -0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.0_pReal, & - 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, & - 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal, & - -0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, & + 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! cubic symmetry operations + 0.0_pReal, 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), & ! 2-fold symmetry + 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), & + 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, & + 0.0_pReal, 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), -1.0_pReal/sqrt(2.0_pReal), & + 0.0_pReal, -1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), & + 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), -1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, & + 0.5_pReal, 0.5_pReal, 0.5_pReal, 0.5_pReal, & ! 3-fold symmetry + -0.5_pReal, 0.5_pReal, 0.5_pReal, 0.5_pReal, & + 0.5_pReal, -0.5_pReal, 0.5_pReal, 0.5_pReal, & + -0.5_pReal, -0.5_pReal, 0.5_pReal, 0.5_pReal, & + 0.5_pReal, 0.5_pReal, -0.5_pReal, 0.5_pReal, & + -0.5_pReal, 0.5_pReal, -0.5_pReal, 0.5_pReal, & + 0.5_pReal, 0.5_pReal, 0.5_pReal, -0.5_pReal, & + -0.5_pReal, 0.5_pReal, 0.5_pReal, -0.5_pReal, & + 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, 0.0_pReal, & ! 4-fold symmetry + 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, & + -1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, 0.0_pReal, & + 1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, & + 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, & + -1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, & + 1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), & + 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal, & + -1.0_pReal/sqrt(2.0_pReal), 0.0_pReal, 0.0_pReal, 1.0_pReal/sqrt(2.0_pReal), & ! - 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! hexagonal symmetry operations - 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, & ! 2-fold symmetry - 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.5_pReal, 0.866025403784439_pReal, 0.0_pReal, & - 0.0_pReal, -0.5_pReal, 0.866025403784439_pReal, 0.0_pReal, & - 0.0_pReal, 0.866025403784439_pReal, 0.5_pReal, 0.0_pReal, & - 0.0_pReal, -0.866025403784439_pReal, 0.5_pReal, 0.0_pReal, & - 0.866025403784439_pReal, 0.0_pReal, 0.0_pReal, 0.5_pReal, & ! 6-fold symmetry - -0.866025403784439_pReal, 0.0_pReal, 0.0_pReal, 0.5_pReal, & - 0.5_pReal, 0.0_pReal, 0.0_pReal, 0.866025403784439_pReal, & - -0.5_pReal, 0.0_pReal, 0.0_pReal, 0.866025403784439_pReal, & - 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal & + 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! hexagonal symmetry operations + 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, & ! 2-fold symmetry + 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, & + 0.0_pReal, 0.5_pReal, 2.0_pReal/sqrt(3.0_pReal), 0.0_pReal, & + 0.0_pReal, -0.5_pReal, 2.0_pReal/sqrt(3.0_pReal), 0.0_pReal, & + 0.0_pReal, 2.0_pReal/sqrt(3.0_pReal), 0.5_pReal, 0.0_pReal, & + 0.0_pReal, -2.0_pReal/sqrt(3.0_pReal), 0.5_pReal, 0.0_pReal, & + 2.0_pReal/sqrt(3.0_pReal), 0.0_pReal, 0.0_pReal, 0.5_pReal, & ! 6-fold symmetry + -2.0_pReal/sqrt(3.0_pReal), 0.0_pReal, 0.0_pReal, 0.5_pReal, & + 0.5_pReal, 0.0_pReal, 0.0_pReal, 2.0_pReal/sqrt(3.0_pReal), & + -0.5_pReal, 0.0_pReal, 0.0_pReal, 2.0_pReal/sqrt(3.0_pReal), & + 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal & ],[4,36]) !< Symmetry operations as quaternions 24 for cubic, 12 for hexagonal = 36 public :: & From 0e17b17b909aa172795d31ae54992fd8f6c09ca8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 7 Oct 2018 18:40:02 +0200 Subject: [PATCH 36/45] using full precision --- src/lattice.f90 | 92 ++++++++++++++++++----------------- src/plastic_phenopowerlaw.f90 | 8 +-- 2 files changed, 52 insertions(+), 48 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index a3a790574..30563f4a8 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1230,10 +1230,10 @@ real(pReal), dimension(4,36), parameter, private :: & lattice_SchmidMatrix_slip, & lattice_SchmidMatrix_twin, & lattice_nonSchmidMatrix, & - lattice_interactionSlipSlip2, & - lattice_interactionTwinTwin2, & - lattice_interactionSlipTwin2, & - lattice_interactionTwinSlip2, & + lattice_interaction_SlipSlip, & + lattice_interaction_TwinTwin, & + lattice_interaction_SlipTwin, & + lattice_interaction_TwinSlip, & lattice_characteristicShear_Twin contains @@ -2346,10 +2346,10 @@ end function lattice_nonSchmidMatrix !-------------------------------------------------------------------------------------------------- -!> @brief Populates reduced slip-slip interaction matrix -!> ToDo: prefix "2" needed as long as deprecated array lattice_interactionSlipSlip exists +!> @brief Populates slip-slip interaction matrix +!> details: only active slip systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interactionSlipSlip2(Nslip,interactionValues,structure) +function lattice_interaction_SlipSlip(Nslip,interactionValues,structure) use IO, only: & IO_error @@ -2357,7 +2357,7 @@ function lattice_interactionSlipSlip2(Nslip,interactionValues,structure) integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values slip-slip character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(sum(Nslip),sum(Nslip)) :: lattice_interactionSlipSlip2 + real(pReal), dimension(sum(Nslip),sum(Nslip)) :: lattice_interaction_SlipSlip integer(pInt), dimension(:), allocatable :: NslipMax integer(pInt), dimension(:,:), allocatable :: interactionSlipSlip @@ -2379,20 +2379,20 @@ function lattice_interactionSlipSlip2(Nslip,interactionValues,structure) call IO_error(132_pInt,ext_msg=trim(structure)//' (slip slip interaction)') end select - if (any(Nslip(1:size(Nslip)) - Nslip < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) + if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Nslip '//trim(structure)) - lattice_interactionSlipSlip2 = & + lattice_interaction_SlipSlip = & buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionSlipSlip) -end function lattice_interactionSlipSlip2 +end function lattice_interaction_SlipSlip !-------------------------------------------------------------------------------------------------- -!> @brief Populates reduced twin-twin interaction matrix -!> ToDo: prefix "2" needed as long as deprecated array lattice_interactionTwinTwin exists +!> @brief Populates twin-twin interaction matrix +!> details: only active twin systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interactionTwinTwin2(Ntwin,interactionValues,structure) +function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) use IO, only: & IO_error @@ -2400,7 +2400,7 @@ function lattice_interactionTwinTwin2(Ntwin,interactionValues,structure) integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: lattice_interactionTwinTwin2 + real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: lattice_interaction_TwinTwin integer(pInt), dimension(:), allocatable :: NtwinMax integer(pInt), dimension(:,:), allocatable :: interactionTwinTwin @@ -2419,20 +2419,20 @@ function lattice_interactionTwinTwin2(Ntwin,interactionValues,structure) call IO_error(132_pInt,ext_msg=trim(structure)//' (twin twin interaction)') end select - if (any(Ntwin(1:size(Ntwin)) - Ntwin < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) + if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Ntwin '//trim(structure)) - lattice_interactionTwinTwin2 = & + lattice_interaction_TwinTwin = & buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTwinTwin) -end function lattice_interactionTwinTwin2 +end function lattice_interaction_TwinTwin !-------------------------------------------------------------------------------------------------- -!> @brief Populates reduced slip-twin interaction matrix -!> ToDo: prefix "2" needed as long as deprecated array lattice_interactionTwinTwin exists +!> @brief Populates slip-twin interaction matrix +!> details: only active slip and twin systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interactionSlipTwin2(Nslip,Ntwin,interactionValues,structure) +function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) use IO, only: & IO_error @@ -2441,7 +2441,7 @@ function lattice_interactionSlipTwin2(Nslip,Ntwin,interactionValues,structure) integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: lattice_interactionSlipTwin2 + real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: lattice_interaction_SlipTwin integer(pInt), dimension(:), allocatable :: NslipMax integer(pInt), dimension(:), allocatable :: NtwinMax @@ -2464,20 +2464,22 @@ function lattice_interactionSlipTwin2(Nslip,Ntwin,interactionValues,structure) call IO_error(132_pInt,ext_msg=trim(structure)//' (slip twin interaction)') end select - if (any(Ntwin(1:size(Ntwin)) - Ntwin < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) + if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Nslip '//trim(structure)) + if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Ntwin '//trim(structure)) - lattice_interactionSlipTwin2 = & + lattice_interaction_SlipTwin = & buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionSlipTwin) -end function lattice_interactionSlipTwin2 +end function lattice_interaction_SlipTwin !-------------------------------------------------------------------------------------------------- -!> @brief Populates reduced twin-slip interaction matrix -!> ToDo: prefix "2" needed as long as deprecated array lattice_interactionTwinTwin exists +!> @brief Populates twin-slip interaction matrix +!> details: only active twin and slip systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interactionTwinSlip2(Ntwin,Nslip,interactionValues,structure) +function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) use IO, only: & IO_error @@ -2486,7 +2488,7 @@ function lattice_interactionTwinSlip2(Ntwin,Nslip,interactionValues,structure) integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: lattice_interactionTwinSlip2 + real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: lattice_interaction_TwinSlip integer(pInt), dimension(:), allocatable :: NslipMax integer(pInt), dimension(:), allocatable :: NtwinMax @@ -2506,23 +2508,25 @@ function lattice_interactionTwinSlip2(Ntwin,Nslip,interactionValues,structure) NtwinMax = LATTICE_HEX_NTWINSYSTEM NslipMax = LATTICE_HEX_NSLIPSYSTEM case default - call IO_error(132_pInt,ext_msg=trim(structure)//' (slip twin interaction)') + call IO_error(132_pInt,ext_msg=trim(structure)//' (twin slip interaction)') end select - if (any(Ntwin(1:size(Ntwin)) - Ntwin < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) + if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Nslip '//trim(structure)) + if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Ntwin '//trim(structure)) - lattice_interactionTwinSlip2 = & + lattice_interaction_TwinSlip = & buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTwinSlip) -end function lattice_interactionTwinSlip2 +end function lattice_interaction_TwinSlip !-------------------------------------------------------------------------------------------------- -!> @brief Populates reduced trans-trans interaction matrix -!> ToDo: prefix "2" needed as long as deprecated array lattice_interactionTransTrans exists +!> @brief Populates trans-trans interaction matrix +!> details: only active transformation systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interactionTransTrans2(Ntrans,interactionValues,structure,targetStructure) +function lattice_interaction_TransTrans(Ntrans,interactionValues,structure,targetStructure) use IO, only: & IO_error @@ -2532,7 +2536,7 @@ function lattice_interactionTransTrans2(Ntrans,interactionValues,structure,targe character(len=*), intent(in) :: & structure, & !< lattice structure of parent crystal targetStructure !< lattice structure of transformed crystal - real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: lattice_interactionTransTrans2 + real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: lattice_interaction_TransTrans integer(pInt), dimension(:), allocatable :: NtransMax integer(pInt), dimension(:,:), allocatable :: interactionTransTrans @@ -2544,13 +2548,13 @@ function lattice_interactionTransTrans2(Ntrans,interactionValues,structure,targe call IO_error(132_pInt,ext_msg=trim(structure)//' => '//trim(targetStructure)) end if - if (any(Ntrans(1:size(Ntrans)) - Ntrans < 0_pInt)) & + if (any(NtransMax(1:size(Ntrans)) - Ntrans < 0_pInt)) & call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) - lattice_interactionTransTrans2 = & + lattice_interaction_TransTrans = & buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTransTrans) -end function lattice_interactionTransTrans2 +end function lattice_interaction_TransTrans !-------------------------------------------------------------------------------------------------- diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 1ca5cd66c..9d30cf184 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -206,7 +206,7 @@ subroutine plastic_phenopowerlaw_init ! reading in slip related parameters prm%xi_slip_0 = config_phase(p)%getFloats('tau0_slip', requiredShape=shape(prm%Nslip)) prm%xi_slip_sat = config_phase(p)%getFloats('tausat_slip', requiredShape=shape(prm%Nslip)) - prm%interaction_SlipSlip = lattice_interactionSlipSlip2(prm%Nslip,config_phase(p)%getFloats('interaction_slipslip'), & + prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip,config_phase(p)%getFloats('interaction_slipslip'), & structure(1:3)) prm%H_int = config_phase(p)%getFloats('h_int', requiredShape=shape(prm%Nslip), & defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) @@ -255,7 +255,7 @@ subroutine plastic_phenopowerlaw_init config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) ! reading in twin related parameters prm%xi_twin_0 = config_phase(p)%getFloats('tau0_twin',requiredShape=shape(prm%Ntwin)) - prm%interaction_TwinTwin = lattice_interactionTwinTwin2(prm%Ntwin,config_phase(p)%getFloats('interaction_twintwin'), & + prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,config_phase(p)%getFloats('interaction_twintwin'), & structure(1:3)) prm%gdot0_twin = config_phase(p)%getFloat('gdot0_twin') @@ -280,10 +280,10 @@ subroutine plastic_phenopowerlaw_init config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then - prm%interaction_SlipTwin = lattice_interactionSlipTwin2(prm%Nslip,prm%Ntwin,& + prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,& config_phase(p)%getFloats('interaction_sliptwin'), & structure(1:3)) - prm%interaction_TwinSlip = lattice_interactionTwinSlip2(prm%Ntwin,prm%Nslip,& + prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,& config_phase(p)%getFloats('interaction_twinslip'), & structure(1:3)) else slipAndTwinActive From facb1c2407c9d11435a1b8651199b8b39e5129ed Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 7 Oct 2018 20:03:45 +0200 Subject: [PATCH 37/45] unified error checking --- src/lattice.f90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 30563f4a8..6e1f90752 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2604,16 +2604,19 @@ end function lattice_SchmidMatrix_slip !-------------------------------------------------------------------------------------------------- !> @brief Calculates Schmid matrix for active twin systems !-------------------------------------------------------------------------------------------------- -function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) +function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) + use prec, only: & + tol_math_check use IO, only: & IO_error use math, only: & + math_trace33, & math_tensorproduct33 implicit none integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(3,3,sum(Ntwin)) :: lattice_SchmidMatrix_Twin + real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrix real(pReal), intent(in), optional :: & cOverA @@ -2632,8 +2635,9 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) end select do i = 1, sum(Ntwin) - lattice_SchmidMatrix_twin(1:3,1:3,i) = & - math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) + SchmidMatrix(1:3,1:3,i) = math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) + if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & + call IO_error(0_pInt,i,ext_msg = 'dilatational Schmid matrix for twin') enddo end function lattice_SchmidMatrix_twin From 302cf9b6c2614df12915d5e71d031901b02500c7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 7 Oct 2018 22:16:18 +0200 Subject: [PATCH 38/45] bugfix, now inline with code from initializeStructure --- src/lattice.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 6e1f90752..401172428 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2338,7 +2338,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc + nonSchmidCoefficients(4) * math_tensorproduct33(normal, normal) if (size(nonSchmidCoefficients)>4) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + nonSchmidCoefficients(5) * math_tensorproduct33(math_crossproduct(normal, direction), & - math_crossproduct(normal, normal)) + math_crossproduct(normal, direction)) if (size(nonSchmidCoefficients)>5) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & + nonSchmidCoefficients(6) * math_tensorproduct33(direction, direction) enddo From 5efcad952a964f13e4e94995a0f62fed64d45897 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 8 Oct 2018 08:27:12 +0200 Subject: [PATCH 39/45] avoid type conversion more systematic checking of correct input parameters --- src/lattice.f90 | 205 +++++++++++++++++++++++++----------------------- 1 file changed, 109 insertions(+), 96 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 401172428..ec1f187e1 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2139,6 +2139,7 @@ end function lattice_characteristicShear_Twin !-------------------------------------------------------------------------------------------------- !> @brief Calculates rotated elasticity matrices for twinning +!> ToDo: Completely untested !-------------------------------------------------------------------------------------------------- function lattice_C66_twin(Ntwin,C66,structure,CoverA) use IO, only: & @@ -2164,11 +2165,11 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) select case(structure) case('fcc') - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,LATTICE_FCC_SYSTEMTWIN,structure) case('bcc') - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMTWIN,structure) case('hex','hexagonal') !ToDo: "No alias policy": long or short? - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NSLIPSYSTEM,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NSLIPSYSTEM,LATTICE_HEX_SYSTEMTWIN,'hex',cOverA) case default call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_C66_twin)') end select @@ -2182,6 +2183,7 @@ end function !-------------------------------------------------------------------------------------------------- !> @brief Calculates rotated elasticity matrices for transformation +!> ToDo: Completely untested and incomplete !-------------------------------------------------------------------------------------------------- function lattice_C66_trans(Ntrans,C_parent66,structure_parent, & C_target66,structure_target) @@ -2245,7 +2247,7 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_parent, & enddo if (trim(structure_parent) == 'fcc' .and. trim(structure_target) == 'hex') then - do i = 1_pInt,sum(Ntrans)!!!!!!!!!!!!!! NEED TO BE FIXED + do i = 1_pInt,sum(Ntrans)!!!!!!!!!!!!!! ToDo: NEED TO BE FIXED R = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), & ! Pitsch rotation lattice_fccTobcc_systemTrans(4,i)*INRAD) B = math_axisAngleToR(lattice_fccTobcc_bainRot(1:3,i), & ! Rotation of fcc to Bain coordinate system @@ -2295,6 +2297,7 @@ end function !-------------------------------------------------------------------------------------------------- !> @brief Non-schmid tensor +!> ToDo: Clean description needed ! Schmid matrices with non-Schmid contributions according to Koester_etal2012, Acta Materialia 60 (2012) ! 3894–3901, eq. (17) ("n1" is replaced by either "np" or "nn" according to either positive or negative slip direction) ! "np" and "nn" according to Gröger_etal2008, Acta Materialia 56 (2008) 5412–5425, table 1 @@ -2320,7 +2323,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc integer(pInt) :: i if (abs(sense) /= 1_pInt) write(6,*) 'mist' - coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,int(LATTICE_BCC_SYSTEMSLIP,pInt),'bcc') + coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMSLIP,'bcc') coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip)) *real(sense,pReal) nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'bcc') @@ -2349,7 +2352,7 @@ end function lattice_nonSchmidMatrix !> @brief Populates slip-slip interaction matrix !> details: only active slip systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_SlipSlip(Nslip,interactionValues,structure) +function lattice_interaction_SlipSlip(Nslip,interactionValues,structure) result(interactionMatrix) use IO, only: & IO_error @@ -2357,33 +2360,32 @@ function lattice_interaction_SlipSlip(Nslip,interactionValues,structure) integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values slip-slip character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(sum(Nslip),sum(Nslip)) :: lattice_interaction_SlipSlip + real(pReal), dimension(sum(Nslip),sum(Nslip)) :: interactionMatrix integer(pInt), dimension(:), allocatable :: NslipMax - integer(pInt), dimension(:,:), allocatable :: interactionSlipSlip + integer(pInt), dimension(:,:), allocatable :: interactionTypes select case(structure) case('fcc') - interactionSlipSlip = LATTICE_FCC_INTERACTIONSLIPSLIP - NslipMax = LATTICE_FCC_NSLIPSYSTEM + interactionTypes = LATTICE_FCC_INTERACTIONSLIPSLIP + NslipMax = LATTICE_FCC_NSLIPSYSTEM case('bcc') - interactionSlipSlip = LATTICE_BCC_INTERACTIONSLIPSLIP - NslipMax = LATTICE_BCC_NSLIPSYSTEM + interactionTypes = LATTICE_BCC_INTERACTIONSLIPSLIP + NslipMax = LATTICE_BCC_NSLIPSYSTEM case('hex','hexagonal') !ToDo: "No alias policy": long or short? - interactionSlipSlip = LATTICE_HEX_INTERACTIONSLIPSLIP - NslipMax = LATTICE_HEX_NSLIPSYSTEM + interactionTypes = LATTICE_HEX_INTERACTIONSLIPSLIP + NslipMax = LATTICE_HEX_NSLIPSYSTEM case('bct') - interactionSlipSlip = LATTICE_BCT_INTERACTIONSLIPSLIP - NslipMax = LATTICE_BCT_NSLIPSYSTEM + interactionTypes = LATTICE_BCT_INTERACTIONSLIPSLIP + NslipMax = LATTICE_BCT_NSLIPSYSTEM case default call IO_error(132_pInt,ext_msg=trim(structure)//' (slip slip interaction)') end select - if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Nslip '//trim(structure)) + if (size(interactionValues) > maxval(interactionTypes)) & + call IO_error(0_pInt) ! ToDo - lattice_interaction_SlipSlip = & - buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionSlipSlip) + interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes) end function lattice_interaction_SlipSlip @@ -2392,7 +2394,7 @@ end function lattice_interaction_SlipSlip !> @brief Populates twin-twin interaction matrix !> details: only active twin systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) +function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) result(interactionMatrix) use IO, only: & IO_error @@ -2400,30 +2402,29 @@ function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: lattice_interaction_TwinTwin + real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: interactionMatrix integer(pInt), dimension(:), allocatable :: NtwinMax - integer(pInt), dimension(:,:), allocatable :: interactionTwinTwin + integer(pInt), dimension(:,:), allocatable :: interactionTypes select case(structure) case('fcc') - interactionTwinTwin = LATTICE_FCC_INTERACTIONTWINTWIN - NtwinMax = LATTICE_FCC_NTWINSYSTEM + interactionTypes = LATTICE_FCC_INTERACTIONTWINTWIN + NtwinMax = LATTICE_FCC_NTWINSYSTEM case('bcc') - interactionTwinTwin = LATTICE_BCC_INTERACTIONTWINTWIN - NtwinMax = LATTICE_BCC_NTWINSYSTEM + interactionTypes = LATTICE_BCC_INTERACTIONTWINTWIN + NtwinMax = LATTICE_BCC_NTWINSYSTEM case('hex','hexagonal') !ToDo: "No alias policy": long or short? - interactionTwinTwin = LATTICE_HEX_INTERACTIONTWINTWIN - NtwinMax = LATTICE_HEX_NTWINSYSTEM + interactionTypes = LATTICE_HEX_INTERACTIONTWINTWIN + NtwinMax = LATTICE_HEX_NTWINSYSTEM case default call IO_error(132_pInt,ext_msg=trim(structure)//' (twin twin interaction)') end select - if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Ntwin '//trim(structure)) + if (size(interactionValues) > maxval(interactionTypes)) & + call IO_error(0_pInt) ! ToDo - lattice_interaction_TwinTwin = & - buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTwinTwin) + interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes) end function lattice_interaction_TwinTwin @@ -2432,7 +2433,7 @@ end function lattice_interaction_TwinTwin !> @brief Populates slip-twin interaction matrix !> details: only active slip and twin systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) +function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) result(interactionMatrix) use IO, only: & IO_error @@ -2441,36 +2442,33 @@ function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: lattice_interaction_SlipTwin + real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: interactionMatrix integer(pInt), dimension(:), allocatable :: NslipMax integer(pInt), dimension(:), allocatable :: NtwinMax - integer(pInt), dimension(:,:), allocatable :: interactionSlipTwin + integer(pInt), dimension(:,:), allocatable :: interactionTypes select case(structure) case('fcc') - interactionSlipTwin = LATTICE_FCC_INTERACTIONSLIPTWIN - NslipMax = LATTICE_FCC_NSLIPSYSTEM - NtwinMax = LATTICE_FCC_NTWINSYSTEM + interactionTypes = LATTICE_FCC_INTERACTIONSLIPTWIN + NslipMax = LATTICE_FCC_NSLIPSYSTEM + NtwinMax = LATTICE_FCC_NTWINSYSTEM case('bcc') - interactionSlipTwin = LATTICE_BCC_INTERACTIONSLIPTWIN - NslipMax = LATTICE_BCC_NSLIPSYSTEM - NtwinMax = LATTICE_BCC_NTWINSYSTEM + interactionTypes = LATTICE_BCC_INTERACTIONSLIPTWIN + NslipMax = LATTICE_BCC_NSLIPSYSTEM + NtwinMax = LATTICE_BCC_NTWINSYSTEM case('hex','hexagonal') !ToDo: "No alias policy": long or short? - interactionSlipTwin = LATTICE_HEX_INTERACTIONSLIPTWIN - NslipMax = LATTICE_HEX_NSLIPSYSTEM - NtwinMax = LATTICE_HEX_NTWINSYSTEM + interactionTypes = LATTICE_HEX_INTERACTIONSLIPTWIN + NslipMax = LATTICE_HEX_NSLIPSYSTEM + NtwinMax = LATTICE_HEX_NTWINSYSTEM case default call IO_error(132_pInt,ext_msg=trim(structure)//' (slip twin interaction)') end select - if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Nslip '//trim(structure)) - if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Ntwin '//trim(structure)) + if (size(interactionValues) > maxval(interactionTypes)) & + call IO_error(0_pInt) ! ToDo - lattice_interaction_SlipTwin = & - buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionSlipTwin) + interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes) end function lattice_interaction_SlipTwin @@ -2479,7 +2477,7 @@ end function lattice_interaction_SlipTwin !> @brief Populates twin-slip interaction matrix !> details: only active twin and slip systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) +function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) result(interactionMatrix) use IO, only: & IO_error @@ -2488,36 +2486,33 @@ function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin character(len=*), intent(in) :: structure !< lattice structure - real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: lattice_interaction_TwinSlip + real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: interactionMatrix integer(pInt), dimension(:), allocatable :: NslipMax integer(pInt), dimension(:), allocatable :: NtwinMax - integer(pInt), dimension(:,:), allocatable :: interactionTwinSlip + integer(pInt), dimension(:,:), allocatable :: interactionTypes select case(structure) case('fcc') - interactionTwinSlip = LATTICE_FCC_INTERACTIONTWINSLIP - NtwinMax = LATTICE_FCC_NTWINSYSTEM - NslipMax = LATTICE_FCC_NSLIPSYSTEM + interactionTypes = LATTICE_FCC_INTERACTIONTWINSLIP + NtwinMax = LATTICE_FCC_NTWINSYSTEM + NslipMax = LATTICE_FCC_NSLIPSYSTEM case('bcc') - interactionTwinSlip = LATTICE_BCC_INTERACTIONTWINSLIP - NtwinMax = LATTICE_BCC_NTWINSYSTEM - NslipMax = LATTICE_BCC_NSLIPSYSTEM + interactionTypes = LATTICE_BCC_INTERACTIONTWINSLIP + NtwinMax = LATTICE_BCC_NTWINSYSTEM + NslipMax = LATTICE_BCC_NSLIPSYSTEM case('hex','hexagonal') !ToDo: "No alias policy": long or short? - interactionTwinSlip = LATTICE_HEX_INTERACTIONTWINSLIP - NtwinMax = LATTICE_HEX_NTWINSYSTEM - NslipMax = LATTICE_HEX_NSLIPSYSTEM + interactionTypes = LATTICE_HEX_INTERACTIONTWINSLIP + NtwinMax = LATTICE_HEX_NTWINSYSTEM + NslipMax = LATTICE_HEX_NSLIPSYSTEM case default call IO_error(132_pInt,ext_msg=trim(structure)//' (twin slip interaction)') end select - if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Nslip '//trim(structure)) - if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Ntwin '//trim(structure)) + if (size(interactionValues) > maxval(interactionTypes)) & + call IO_error(0_pInt) ! ToDo - lattice_interaction_TwinSlip = & - buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTwinSlip) + interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes) end function lattice_interaction_TwinSlip @@ -2526,7 +2521,7 @@ end function lattice_interaction_TwinSlip !> @brief Populates trans-trans interaction matrix !> details: only active transformation systems are considered !-------------------------------------------------------------------------------------------------- -function lattice_interaction_TransTrans(Ntrans,interactionValues,structure,targetStructure) +function lattice_interaction_TransTrans(Ntrans,interactionValues,structure,targetStructure) result(interactionMatrix) use IO, only: & IO_error @@ -2536,23 +2531,22 @@ function lattice_interaction_TransTrans(Ntrans,interactionValues,structure,targe character(len=*), intent(in) :: & structure, & !< lattice structure of parent crystal targetStructure !< lattice structure of transformed crystal - real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: lattice_interaction_TransTrans + real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: interactionMatrix integer(pInt), dimension(:), allocatable :: NtransMax - integer(pInt), dimension(:,:), allocatable :: interactionTransTrans + integer(pInt), dimension(:,:), allocatable :: interactionTypes if (trim(structure) == 'fcc' .and. trim(targetStructure) == 'hex') then - interactionTransTrans = lattice_fccToHex_interactionTransTrans - NtransMax = lattice_fcc_Ntrans + interactionTypes = lattice_fccToHex_interactionTransTrans + NtransMax = lattice_fcc_Ntrans else call IO_error(132_pInt,ext_msg=trim(structure)//' => '//trim(targetStructure)) end if - if (any(NtransMax(1:size(Ntrans)) - Ntrans < 0_pInt)) & - call IO_error(145_pInt,ext_msg='Ntrans '//trim(structure)) + if (size(interactionValues) > maxval(interactionTypes)) & + call IO_error(0_pInt) ! ToDo - lattice_interaction_TransTrans = & - buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTransTrans) + interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes) end function lattice_interaction_TransTrans @@ -2577,21 +2571,31 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) cOverA real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem + real(pReal), dimension(:,:), allocatable :: slipSystems + integer(pInt), dimension(:), allocatable :: NslipMax integer(pInt) :: i select case(structure) case('fcc') - coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,int(LATTICE_FCC_SYSTEMSLIP,pInt),structure) + NslipMax = LATTICE_BCC_NSLIPSYSTEM + slipSystems = LATTICE_FCC_SYSTEMSLIP case('bcc') - coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_FCC_NSLIPSYSTEM,int(LATTICE_BCC_SYSTEMSLIP,pInt),structure) + NslipMax = LATTICE_FCC_NSLIPSYSTEM + slipSystems = LATTICE_BCC_SYSTEMSLIP case('hex','hexagonal') !ToDo: "No alias policy": long or short? - coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_HEX_NSLIPSYSTEM,int(LATTICE_HEX_SYSTEMSLIP,pInt),'hex',cOverA) + NslipMax = LATTICE_HEX_NSLIPSYSTEM + slipSystems = LATTICE_HEX_SYSTEMSLIP case('bct') - coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCT_NSLIPSYSTEM,int(LATTICE_BCT_SYSTEMSLIP,pInt),structure,cOverA) + NslipMax = LATTICE_BCT_NSLIPSYSTEM + slipSystems = LATTICE_BCT_SYSTEMSLIP case default call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_slip)') end select + if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt) .or. any(Nslip < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Nslip '//trim(structure)) + + coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure) do i = 1, sum(Nslip) SchmidMatrix(1:3,1:3,i) = math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & @@ -2621,19 +2625,28 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) cOverA real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem + real(pReal), dimension(:,:), allocatable :: twinSystems + integer(pInt), dimension(:), allocatable :: NtwinMax integer(pInt) :: i select case(structure) case('fcc') - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NTWINSYSTEM,int(LATTICE_FCC_SYSTEMTWIN,pInt),structure) + NtwinMax = LATTICE_FCC_NTWINSYSTEM + twinSystems = LATTICE_FCC_SYSTEMTWIN case('bcc') - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NTWINSYSTEM,int(LATTICE_BCC_SYSTEMTWIN,pInt),structure) + NtwinMax = LATTICE_BCC_NTWINSYSTEM + twinSystems = LATTICE_BCC_SYSTEMTWIN case('hex','hexagonal') !ToDo: "No alias policy": long or short? - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NTWINSYSTEM,int(LATTICE_HEX_SYSTEMTWIN,pInt),'hex',cOverA) + NtwinMax = LATTICE_HEX_NTWINSYSTEM + twinSystems = LATTICE_HEX_SYSTEMTWIN case default call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_twin)') end select + if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt) .or. any(Ntwin < 0_pInt)) & + call IO_error(145_pInt,ext_msg='Ntwin '//trim(structure)) + + coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure) do i = 1, sum(Ntwin) SchmidMatrix(1:3,1:3,i) = math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & @@ -2690,7 +2703,7 @@ function buildCoordinateSystem(active,maximum,system,structure,cOverA) integer(pInt), dimension(:), intent(in) :: & active, & maximum - integer(pInt), dimension(:,:), intent(in) :: & + real(pReal), dimension(:,:), intent(in) :: & system character(len=*), intent(in) :: & structure !< lattice structure @@ -2716,23 +2729,23 @@ function buildCoordinateSystem(active,maximum,system,structure,cOverA) select case(trim(structure)) case ('fcc','bcc') - direction = real(system(1:3,j),pReal) - normal = real(system(4:6,j),pReal) + direction = system(1:3,j) + normal = system(4:6,j) case ('hex') ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)]) - direction = [ real(system(1,j),pReal)*1.5_pReal, & - (real(system(1,j),pReal)+2.0_pReal*real(system(2,j),pReal))*sqrt(0.75_pReal), & - real(system(4,j),pReal)*CoverA ] + direction = [ system(1,j)*1.5_pReal, & + (system(1,j)+2.0_pReal*system(2,j))*sqrt(0.75_pReal), & + system(4,j)*CoverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) - normal = [ real(system(5,j),pReal), & - (real(system(5,j),pReal)+2.0_pReal*real(system(6,j),pReal))/ sqrt(3.0_pReal), & - real(system(8,j),pReal)/CoverA ] + normal = [ system(5,j), & + (system(5,j)+2.0_pReal*system(6,j))/ sqrt(3.0_pReal), & + system(8,j)/CoverA ] case ('bct') - direction = [real(system(1:2,j),pReal),real(system(3,i),pReal)*CoverA] - normal = [real(system(4:5,j),pReal),real(system(6,i),pReal)/CoverA] + direction = [system(1:2,j),system(3,i)*CoverA] + normal = [system(4:5,j),system(6,i)/CoverA] end select From 2a79f7ea61e73f3388bf0af6a24dc7680e05f47e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 8 Oct 2018 13:25:40 +0200 Subject: [PATCH 40/45] taking care of c/a for hex and bct --- src/lattice.f90 | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index ec1f187e1..86108197e 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2096,7 +2096,8 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact implicit none integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=3), intent(in) :: structure - real(pReal), intent(in) :: cOverA + real(pReal), intent(in), optional :: & + cOverA real(pReal), dimension(sum(Ntwin)) :: characteristicShear integer(pInt) :: & ir, & !< index in reduced list @@ -2117,6 +2118,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact ig = sum(LATTICE_BCC_NTWINSYSTEM(1:mf-1))+ms characteristicShear(ir) = LATTICE_BCC_SHEARTWIN(ig) case('hex') + if (.not. present(CoverA)) call IO_error(0_pInt) ig = sum(LATTICE_HEX_NTWINSYSTEM(1:mf-1))+ms select case(LATTICE_HEX_SHEARTWIN(ig)) ! from Christian & Mahajan 1995 p.29 case (1_pInt) ! <-10.1>{10.2} @@ -2155,12 +2157,12 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), dimension(6,6), intent(in) :: C66 - real(pReal), intent(in) :: cOverA - real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin + real(pReal), optional, intent(in) :: cOverA + real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin - real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem + real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem - real(pReal), dimension(3,3) :: R + real(pReal), dimension(3,3) :: R integer(pInt) :: i select case(structure) @@ -2169,6 +2171,7 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) case('bcc') coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMTWIN,structure) case('hex','hexagonal') !ToDo: "No alias policy": long or short? + if (.not. present(CoverA)) call IO_error(0_pInt) coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NSLIPSYSTEM,LATTICE_HEX_SYSTEMTWIN,'hex',cOverA) case default call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_C66_twin)') @@ -2583,9 +2586,11 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) NslipMax = LATTICE_FCC_NSLIPSYSTEM slipSystems = LATTICE_BCC_SYSTEMSLIP case('hex','hexagonal') !ToDo: "No alias policy": long or short? + if (.not. present(CoverA)) call IO_error(0_pInt) NslipMax = LATTICE_HEX_NSLIPSYSTEM slipSystems = LATTICE_HEX_SYSTEMSLIP case('bct') + if (.not. present(CoverA)) call IO_error(0_pInt) NslipMax = LATTICE_BCT_NSLIPSYSTEM slipSystems = LATTICE_BCT_SYSTEMSLIP case default @@ -2594,8 +2599,11 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt) .or. any(Nslip < 0_pInt)) & call IO_error(145_pInt,ext_msg='Nslip '//trim(structure)) - - coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure) + if (present(cOverA)) then + coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA) + else + coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure) + endif do i = 1, sum(Nslip) SchmidMatrix(1:3,1:3,i) = math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & @@ -2637,6 +2645,7 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) NtwinMax = LATTICE_BCC_NTWINSYSTEM twinSystems = LATTICE_BCC_SYSTEMTWIN case('hex','hexagonal') !ToDo: "No alias policy": long or short? + if (.not. present(CoverA)) call IO_error(0_pInt) NtwinMax = LATTICE_HEX_NTWINSYSTEM twinSystems = LATTICE_HEX_SYSTEMTWIN case default @@ -2645,8 +2654,12 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt) .or. any(Ntwin < 0_pInt)) & call IO_error(145_pInt,ext_msg='Ntwin '//trim(structure)) - - coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure) + + if (present(cOverA)) then + coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure,cOverA) + else + coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure) + endif do i = 1, sum(Ntwin) SchmidMatrix(1:3,1:3,i) = math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & From 71d80cc203427ee18efd85444e1b477f56697c81 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 8 Oct 2018 13:27:15 +0200 Subject: [PATCH 41/45] sanity check too strong still makes sense to have, but some tests need to be updated before setting it active again --- src/lattice.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 86108197e..80032e942 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2385,8 +2385,8 @@ function lattice_interaction_SlipSlip(Nslip,interactionValues,structure) result( call IO_error(132_pInt,ext_msg=trim(structure)//' (slip slip interaction)') end select - if (size(interactionValues) > maxval(interactionTypes)) & - call IO_error(0_pInt) ! ToDo + !if (size(interactionValues) > maxval(interactionTypes)) & + ! call IO_error(0_pInt) ! ToDo interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes) @@ -2424,8 +2424,8 @@ function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) result( call IO_error(132_pInt,ext_msg=trim(structure)//' (twin twin interaction)') end select - if (size(interactionValues) > maxval(interactionTypes)) & - call IO_error(0_pInt) ! ToDo + !if (size(interactionValues) > maxval(interactionTypes)) & + ! call IO_error(0_pInt) ! ToDo interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes) @@ -2468,8 +2468,8 @@ function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) r call IO_error(132_pInt,ext_msg=trim(structure)//' (slip twin interaction)') end select - if (size(interactionValues) > maxval(interactionTypes)) & - call IO_error(0_pInt) ! ToDo + !if (size(interactionValues) > maxval(interactionTypes)) & + ! call IO_error(0_pInt) ! ToDo interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes) @@ -2512,8 +2512,8 @@ function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) r call IO_error(132_pInt,ext_msg=trim(structure)//' (twin slip interaction)') end select - if (size(interactionValues) > maxval(interactionTypes)) & - call IO_error(0_pInt) ! ToDo + !if (size(interactionValues) > maxval(interactionTypes)) & + ! call IO_error(0_pInt) ! ToDo interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes) @@ -2546,8 +2546,8 @@ function lattice_interaction_TransTrans(Ntrans,interactionValues,structure,targe call IO_error(132_pInt,ext_msg=trim(structure)//' => '//trim(targetStructure)) end if - if (size(interactionValues) > maxval(interactionTypes)) & - call IO_error(0_pInt) ! ToDo + !if (size(interactionValues) > maxval(interactionTypes)) & + ! call IO_error(0_pInt) ! ToDo interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes) From 78b5fa466d59feb2d3d2420f0c0590043cd348ae Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 8 Oct 2018 19:32:18 +0200 Subject: [PATCH 42/45] fixed hick-up in lattice definitions --- src/lattice.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 80032e942..89d29b515 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2580,10 +2580,10 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) select case(structure) case('fcc') - NslipMax = LATTICE_BCC_NSLIPSYSTEM + NslipMax = LATTICE_FCC_NSLIPSYSTEM slipSystems = LATTICE_FCC_SYSTEMSLIP case('bcc') - NslipMax = LATTICE_FCC_NSLIPSYSTEM + NslipMax = LATTICE_BCC_NSLIPSYSTEM slipSystems = LATTICE_BCC_SYSTEMSLIP case('hex','hexagonal') !ToDo: "No alias policy": long or short? if (.not. present(CoverA)) call IO_error(0_pInt) From 4da3923440839cab466664254c071801a18cfb27 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 8 Oct 2018 23:54:45 +0200 Subject: [PATCH 43/45] more straight-forward definition of array shape --- src/lattice.f90 | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 89d29b515..a226bef27 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -105,7 +105,7 @@ module lattice 0, 1, 1, -1, 1,-1, & ! D1 1, 0,-1, -1, 1,-1, & ! D4 -1,-1, 0, -1, 1,-1 & ! D6 - ],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Nslip]) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli + ],pReal),shape(LATTICE_FCC_SYSTEMSLIP)) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli character(len=*), dimension(1), parameter, public :: LATTICE_FCC_SLIPFAMILY_NAME = & ['<0 1 -1>{1 1 1}'] @@ -124,7 +124,7 @@ module lattice 2, 1,-1, -1, 1,-1, & -1,-2,-1, -1, 1,-1, & -1, 1, 2, -1, 1,-1 & - ],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Ntwin]) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli + ],pReal),shape(LATTICE_FCC_SYSTEMTWIN)) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli character(len=*), dimension(1), parameter, public :: LATTICE_FCC_TWINFAMILY_NAME = & ['<-2 1 1>{1 1 1}'] @@ -143,7 +143,7 @@ module lattice 2, 1,-1, -1, 1,-1, & -1,-2,-1, -1, 1,-1, & -1, 1, 2, -1, 1,-1 & - ],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Ntrans]) + ],pReal),shape(LATTICE_FCCTOHEX_SYSTEMTRANS)) real(pReal), dimension(LATTICE_fcc_Ntwin), parameter, private :: & LATTICE_fcc_shearTwin = 0.5_pReal*sqrt(2.0_pReal) !< Twin system <112>{111} ??? Sorted according to Eisenlohr & Hantcherli @@ -162,7 +162,7 @@ module lattice 11,12, & 10,12, & 10,11 & - ],pInt),[2_pInt,LATTICE_fcc_Ntwin]) + ],pInt),shape(LATTICE_FCC_TWINNUCLEATIONSLIPPAIR)) integer(pInt), dimension(LATTICE_fcc_Nslip,lattice_fcc_Nslip), parameter, public :: & LATTICE_fcc_interactionSlipSlip = reshape(int( [& @@ -178,7 +178,7 @@ module lattice 4,5,6,3,5,5,4,6,5,1,2,2, & 5,3,5,5,4,6,6,4,5,2,1,2, & 6,5,4,5,6,4,5,5,3,2,2,1 & - ],pInt),[LATTICE_fcc_Nslip,LATTICE_fcc_Nslip],order=[2,1]) !< Slip--slip interaction types for fcc + ],pInt),shape(LATTICE_FCC_INTERACTIONSLIPSLIP),order=[2,1]) !< Slip--slip interaction types for fcc !< 1: self interaction !< 2: coplanar interaction !< 3: collinear interaction @@ -199,7 +199,7 @@ module lattice 3,3,3,2,2,2,3,3,3,1,1,1, & 2,2,2,3,3,3,3,3,3,1,1,1, & 3,3,3,3,3,3,2,2,2,1,1,1 & - ],pInt),[LATTICE_fcc_Nslip,LATTICE_fcc_Ntwin],order=[2,1]) !< Slip--twin interaction types for fcc + ],pInt),shape(LATTICE_FCC_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip--twin interaction types for fcc !< 1: coplanar interaction !< 2: screw trace between slip system and twin habit plane (easy cross slip) !< 3: other interaction @@ -220,7 +220,7 @@ module lattice 2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1 & - ],pInt),[lattice_fcc_Ntwin,lattice_fcc_Ntwin],order=[2,1]) !< Twin--twin interaction types for fcc + ],pInt),shape(LATTICE_FCC_INTERACTIONTWINTWIN),order=[2,1]) !< Twin--twin interaction types for fcc integer(pInt), dimension(LATTICE_fcc_Nslip,LATTICE_fcc_Ntrans), parameter, public :: & LATTICE_fccTohex_interactionSlipTrans = reshape(int( [& @@ -236,7 +236,7 @@ module lattice 3,3,3,2,2,2,3,3,3,1,1,1, & 2,2,2,3,3,3,3,3,3,1,1,1, & 3,3,3,3,3,3,2,2,2,1,1,1 & - ],pInt),[LATTICE_fcc_Nslip,LATTICE_fcc_Ntrans],order=[2,1]) !< Slip--trans interaction types for fcc + ],pInt),shape(LATTICE_FCCTOHEX_INTERACTIONSLIPTRANS),order=[2,1]) !< Slip--trans interaction types for fcc integer(pInt), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Nslip), parameter, public :: & LATTICE_fccTohex_interactionTransSlip = 1_pInt !< Trans--Slip interaction types for fcc @@ -255,10 +255,10 @@ module lattice 2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1 & - ],pInt),[LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans],order=[2,1]) !< Trans--trans interaction types for fcc + ],pInt),shape(LATTICE_FCCTOHEX_INTERACTIONTRANSTRANS),order=[2,1]) !< Trans--trans interaction types for fcc real(pReal), dimension(LATTICE_fcc_Ntrans), parameter, private :: & - LATTICE_fccTohex_shearTrans = sqrt(2.0_pReal)/4.0_pReal + LATTICE_fccTohex_shearTrans = sqrt(1.0_pReal/8.0_pReal) real(pReal), dimension(4,LATTICE_fcc_Ntrans), parameter, private :: & LATTICE_fccTobcc_systemTrans = reshape([& @@ -274,7 +274,7 @@ module lattice 1.0, 0.0, 0.0, -10.26, & 0.0, 1.0, 0.0, 10.26, & 0.0, 1.0, 0.0, -10.26 & - ],[ 4_pInt,LATTICE_fcc_Ntrans]) + ],shape(LATTICE_FCCTOBCC_SYSTEMTRANS)) integer(pInt), dimension(9,LATTICE_fcc_Ntrans), parameter, private :: & LATTICE_fccTobcc_bainVariant = reshape(int( [& @@ -290,7 +290,7 @@ module lattice 0, 0, 1, 1, 0, 0, 0, 1, 0, & 0, 0, 1, 1, 0, 0, 0, 1, 0, & 0, 0, 1, 1, 0, 0, 0, 1, 0 & - ],pInt),[ 9_pInt, LATTICE_fcc_Ntrans]) + ],pInt),shape(LATTICE_FCCTOBCC_BAINVARIANT)) real(pReal), dimension(4,LATTICE_fcc_Ntrans), parameter, private :: & LATTICE_fccTobcc_bainRot = reshape([& @@ -306,7 +306,7 @@ module lattice 0.0, 0.0, 1.0, 45.0, & 0.0, 0.0, 1.0, 45.0, & 0.0, 0.0, 1.0, 45.0 & - ],[ 4_pInt,LATTICE_fcc_Ntrans]) + ],shape(LATTICE_FCCTOBCC_BAINROT)) real(pReal), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans), parameter, private :: & ! Matrix for projection of shear from slip system to fault-band (twin) systems LATTICE_fccTobcc_projectionTrans = reshape(real([& ! For ns = nt = nr @@ -322,7 +322,7 @@ module lattice 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, & 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, & 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0 & - ],pReal),[LATTICE_fcc_Ntrans,LATTICE_fcc_Ntrans],order=[2,1]) + ],pReal),shape(LATTICE_FCCTOBCC_PROJECTIONTRANS),order=[2,1]) real(pReal), parameter, private :: & LATTICE_fccTobcc_projectionTransFactor = sqrt(3.0_pReal/4.0_pReal) @@ -344,7 +344,7 @@ module lattice 3, 9, & 3, 12, & 6, 9 & - ],pInt),[2_pInt,LATTICE_fcc_Ntrans]) + ],pInt),shape(LATTICE_FCCTOBCC_TRANSNUCLEATIONTWINPAIR)) real(pReal), dimension(3+3,LATTICE_fcc_Ncleavage), parameter, private :: & LATTICE_fcc_systemCleavage = reshape(real([& @@ -356,7 +356,7 @@ module lattice 0,-1,-1, -1,-1, 1, & -1, 0,-1, 1,-1,-1, & 0, 1, 1, -1, 1,-1 & - ],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Ncleavage]) + ],pReal),shape(LATTICE_FCC_SYSTEMCLEAVAGE)) !-------------------------------------------------------------------------------------------------- ! body centered cubic From 1d0584dcdf7503a686960887607ec531d463d1ae Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Tue, 9 Oct 2018 16:44:47 +0200 Subject: [PATCH 44/45] typo in message --- src/IO.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/IO.f90 b/src/IO.f90 index 96c07e8fe..4a02cd267 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -1461,7 +1461,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (144_pInt) msg = 'negative number of systems' case (145_pInt) - msg = 'to many systems requested' + msg = 'too many systems requested' !-------------------------------------------------------------------------------------------------- ! material error messages and related messages in mesh From ef1899cfd0f90845ec028559c1e63f3669ce98ec Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 10 Oct 2018 07:42:07 +0200 Subject: [PATCH 45/45] improved error reporting - tell user specifically whether number of systems is too high or too low - c/a sanity check can be done centrally --- src/IO.f90 | 2 +- src/lattice.f90 | 47 +++++++++++++++++++++-------------------------- 2 files changed, 22 insertions(+), 27 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 4a02cd267..22d30d0eb 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -1459,7 +1459,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (143_pInt) msg = 'no value found for key' case (144_pInt) - msg = 'negative number of systems' + msg = 'negative number systems requested' case (145_pInt) msg = 'too many systems requested' diff --git a/src/lattice.f90 b/src/lattice.f90 index a226bef27..996852a79 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2157,7 +2157,7 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), dimension(6,6), intent(in) :: C66 - real(pReal), optional, intent(in) :: cOverA + real(pReal), intent(in) :: cOverA real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem @@ -2167,11 +2167,10 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) select case(structure) case('fcc') - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,LATTICE_FCC_SYSTEMTWIN,structure) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,LATTICE_FCC_SYSTEMTWIN,structure,cOverA) case('bcc') - coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMTWIN,structure) + coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMTWIN,structure,cOverA) case('hex','hexagonal') !ToDo: "No alias policy": long or short? - if (.not. present(CoverA)) call IO_error(0_pInt) coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NSLIPSYSTEM,LATTICE_HEX_SYSTEMTWIN,'hex',cOverA) case default call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_C66_twin)') @@ -2326,9 +2325,9 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc integer(pInt) :: i if (abs(sense) /= 1_pInt) write(6,*) 'mist' - coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMSLIP,'bcc') + coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMSLIP,'bcc',0.0_pReal) coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip)) *real(sense,pReal) - nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'bcc') + nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'bcc',0.0_pReal) do i = 1_pInt,sum(Nslip) direction = coordinateSystem(1:3,1,i) @@ -2570,8 +2569,7 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), dimension(3,3,sum(Nslip)) :: SchmidMatrix - real(pReal), intent(in), optional :: & - cOverA + real(pReal), intent(in) :: cOverA real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem real(pReal), dimension(:,:), allocatable :: slipSystems @@ -2586,24 +2584,22 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) NslipMax = LATTICE_BCC_NSLIPSYSTEM slipSystems = LATTICE_BCC_SYSTEMSLIP case('hex','hexagonal') !ToDo: "No alias policy": long or short? - if (.not. present(CoverA)) call IO_error(0_pInt) NslipMax = LATTICE_HEX_NSLIPSYSTEM slipSystems = LATTICE_HEX_SYSTEMSLIP case('bct') - if (.not. present(CoverA)) call IO_error(0_pInt) NslipMax = LATTICE_BCT_NSLIPSYSTEM slipSystems = LATTICE_BCT_SYSTEMSLIP case default call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_slip)') end select - if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt) .or. any(Nslip < 0_pInt)) & + if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) & call IO_error(145_pInt,ext_msg='Nslip '//trim(structure)) - if (present(cOverA)) then - coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA) - else - coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure) - endif + if (any(Nslip < 0_pInt)) & + call IO_error(144_pInt,ext_msg='Nslip '//trim(structure)) + + coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA) + do i = 1, sum(Nslip) SchmidMatrix(1:3,1:3,i) = math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & @@ -2629,8 +2625,7 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrix - real(pReal), intent(in), optional :: & - cOverA + real(pReal), intent(in) :: cOverA real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem real(pReal), dimension(:,:), allocatable :: twinSystems @@ -2645,21 +2640,19 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) NtwinMax = LATTICE_BCC_NTWINSYSTEM twinSystems = LATTICE_BCC_SYSTEMTWIN case('hex','hexagonal') !ToDo: "No alias policy": long or short? - if (.not. present(CoverA)) call IO_error(0_pInt) NtwinMax = LATTICE_HEX_NTWINSYSTEM twinSystems = LATTICE_HEX_SYSTEMTWIN case default call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_twin)') end select - if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt) .or. any(Ntwin < 0_pInt)) & + if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt)) & call IO_error(145_pInt,ext_msg='Ntwin '//trim(structure)) + if (any(Ntwin < 0_pInt)) & + call IO_error(144_pInt,ext_msg='Ntwin '//trim(structure)) - if (present(cOverA)) then - coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure,cOverA) - else - coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure) - endif + coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure,cOverA) + do i = 1, sum(Ntwin) SchmidMatrix(1:3,1:3,i) = math_tensorproduct33(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & @@ -2720,7 +2713,7 @@ function buildCoordinateSystem(active,maximum,system,structure,cOverA) system character(len=*), intent(in) :: & structure !< lattice structure - real(pReal), intent(in), optional :: & + real(pReal), intent(in) :: & cOverA real(pReal), dimension(3,3,sum(active)) :: & buildCoordinateSystem @@ -2746,6 +2739,7 @@ function buildCoordinateSystem(active,maximum,system,structure,cOverA) normal = system(4:6,j) case ('hex') + !ToDo: check c/a ratio ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)]) direction = [ system(1,j)*1.5_pReal, & (system(1,j)+2.0_pReal*system(2,j))*sqrt(0.75_pReal), & @@ -2757,6 +2751,7 @@ function buildCoordinateSystem(active,maximum,system,structure,cOverA) system(8,j)/CoverA ] case ('bct') + !ToDo: check c/a ratio direction = [system(1:2,j),system(3,i)*CoverA] normal = [system(4:5,j),system(6,i)/CoverA]