223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

detabbed file for better readability

This commit is contained in:
Christoph Kords 2010-02-19 13:44:38 +00:00
parent 8c5852dedf
commit aab5598d96
1 changed files with 132 additions and 137 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

269
code/crystallite.f90
View File

@ -347,8 +347,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
storedConvergenceFlag storedConvergenceFlag
logical, dimension(3,3) :: mask logical, dimension(3,3) :: mask
logical forceLocalStiffnessCalculation logical forceLocalStiffnessCalculation ! flag indicating that stiffness calculation is always done locally
forceLocalStiffnessCalculation = .true. ! flag indicating that stiffness calculation is always done locally forceLocalStiffnessCalculation = .true.
! ------ initialize to starting condition ------ ! ------ initialize to starting condition ------
@ -547,7 +548,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_todo = crystallite_todo .and. crystallite_onTrack ! continue with non-broken grains crystallite_todo = crystallite_todo .and. crystallite_onTrack ! continue with non-broken grains
if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & ! any non-local is broken? if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & ! any non-local is broken?
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -635,7 +636,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged? if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged?
crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! all non-local not converged crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! all non-local not converged
crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged
@ -662,11 +663,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! call IO_warning(600,e,i,g) ! call IO_warning(600,e,i,g)
invFp = math_inv3x3(crystallite_partionedFp0(:,:,g,i,e)) invFp = math_inv3x3(crystallite_partionedFp0(:,:,g,i,e))
Fe_guess = math_mul33x33(crystallite_partionedF(:,:,g,i,e),invFp) Fe_guess = math_mul33x33(crystallite_partionedF(:,:,g,i,e),invFp)
Tstar = math_Mandel6to33( & Tstar = math_Mandel6to33( math_mul66x6( 0.5_pReal*constitutive_homogenizedC(g,i,e), &
math_mul66x6( 0.5_pReal*constitutive_homogenizedC(g,i,e), & math_Mandel33to6( math_mul33x33(transpose(Fe_guess),Fe_guess) - math_I3 ) ) )
math_Mandel33to6( math_mul33x33(transpose(Fe_guess),Fe_guess) - math_I3 ) &
) &
)
crystallite_P(:,:,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp))) crystallite_P(:,:,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp)))
endif endif
enddo enddo
@ -690,15 +688,15 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
storedState(1:mySizeState,g,i,e) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state, ... storedState(1:mySizeState,g,i,e) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state, ...
storedDotState(1:mySizeDotState,g,i,e) = constitutive_dotState(g,i,e)%p ! ... dotStates, ... storedDotState(1:mySizeDotState,g,i,e) = constitutive_dotState(g,i,e)%p ! ... dotStates, ...
enddo; enddo; enddo enddo; enddo; enddo
storedTemperature = crystallite_Temperature ! ... Temperature, ... storedTemperature = crystallite_Temperature ! ... Temperature, ...
storedF = crystallite_subF ! ... and kinematics storedF = crystallite_subF ! ... and kinematics
storedFp = crystallite_Fp storedFp = crystallite_Fp
storedInvFp = crystallite_invFp storedInvFp = crystallite_invFp
storedFe = crystallite_Fe storedFe = crystallite_Fe
storedLp = crystallite_Lp storedLp = crystallite_Lp
storedTstar_v = crystallite_Tstar_v storedTstar_v = crystallite_Tstar_v
storedP = crystallite_P storedP = crystallite_P
storedConvergenceFlag = crystallite_converged storedConvergenceFlag = crystallite_converged
if (all(crystallite_localConstitution) .or. theInc < 2 .or. forceLocalStiffnessCalculation) then ! all grains have local constitution, so local convergence of perturbed grain is sufficient if (all(crystallite_localConstitution) .or. theInc < 2 .or. forceLocalStiffnessCalculation) then ! all grains have local constitution, so local convergence of perturbed grain is sufficient
@ -783,11 +781,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo ! perturbation direction enddo ! perturbation direction
select case(pert_method) select case(pert_method)
case (1) case (1)
crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,1) crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,1)
case (2) case (2)
crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,2) crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,2)
case (3) case (3)
crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_perturbation(:,:,:,:,1)+dPdF_perturbation(:,:,:,:,2)) crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_perturbation(:,:,:,:,1)+dPdF_perturbation(:,:,:,:,2))
end select end select
else ! grain did not converge else ! grain did not converge
crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback
@ -818,8 +816,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
NiterationState = 0_pInt NiterationState = 0_pInt
crystallite_todo = .true. crystallite_todo = .true.
do while ( any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) & do while ( any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) .and. NiterationState < nState)
.and. NiterationState < nState)
NiterationState = NiterationState + 1_pInt NiterationState = NiterationState + 1_pInt
do ee = FEsolving_execElem(1),FEsolving_execElem(2) do ee = FEsolving_execElem(1),FEsolving_execElem(2)
@ -828,7 +825,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do gg = 1,myNgrains do gg = 1,myNgrains
if (crystallite_todo(gg,ii,ee)) & if (crystallite_todo(gg,ii,ee)) &
crystallite_onTrack(gg,ii,ee) = crystallite_integrateStress(gg,ii,ee) ! stress integration crystallite_onTrack(gg,ii,ee) = crystallite_integrateStress(gg,ii,ee) ! stress integration
enddo; enddo; enddo enddo; enddo; enddo
crystallite_todo = crystallite_todo .and. crystallite_onTrack ! continue with non-broken grains crystallite_todo = crystallite_todo .and. crystallite_onTrack ! continue with non-broken grains
if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & ! any non-local is broken? if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & ! any non-local is broken?
@ -863,9 +860,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_converged(gg,ii,ee) = crystallite_stateConverged(gg,ii,ee) & crystallite_converged(gg,ii,ee) = crystallite_stateConverged(gg,ii,ee) &
.and. crystallite_temperatureConverged(gg,ii,ee) .and. crystallite_temperatureConverged(gg,ii,ee)
endif endif
enddo enddo; enddo; enddo
enddo
enddo
if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged? if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged?
crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! all non-local not converged crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! all non-local not converged
@ -889,14 +884,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
constitutive_state(gg,ii,ee)%p = storedState(1:mySizeState,gg,ii,ee) constitutive_state(gg,ii,ee)%p = storedState(1:mySizeState,gg,ii,ee)
constitutive_dotState(gg,ii,ee)%p = storedDotState(1:mySizeDotState,gg,ii,ee) constitutive_dotState(gg,ii,ee)%p = storedDotState(1:mySizeDotState,gg,ii,ee)
enddo; enddo; enddo enddo; enddo; enddo
crystallite_Temperature = storedTemperature crystallite_Temperature = storedTemperature
crystallite_subF = storedF crystallite_subF = storedF
crystallite_Fp = storedFp crystallite_Fp = storedFp
crystallite_invFp = storedInvFp crystallite_invFp = storedInvFp
crystallite_Fe = storedFe crystallite_Fe = storedFe
crystallite_Lp = storedLp crystallite_Lp = storedLp
crystallite_Tstar_v = storedTstar_v crystallite_Tstar_v = storedTstar_v
crystallite_P = storedP crystallite_P = storedP
!$OMP CRITICAL (out) !$OMP CRITICAL (out)
debug_StiffnessStateLoopDistribution(NiterationState) = debug_StiffnessstateLoopDistribution(NiterationState) + 1 debug_StiffnessStateLoopDistribution(NiterationState) = debug_StiffnessstateLoopDistribution(NiterationState) + 1
@ -1377,28 +1372,28 @@ LpLoop: do
subroutine crystallite_orientations() subroutine crystallite_orientations()
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pInt, & use prec, only: pInt, &
pReal pReal
use math, only: math_pDecomposition, & use math, only: math_pDecomposition, &
math_RtoEuler, & math_RtoEuler, &
math_misorientation, & math_misorientation, &
inDeg inDeg
use FEsolving, only: FEsolving_execElem, & use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP FEsolving_execIP
use IO, only: IO_warning use IO, only: IO_warning
use material, only: material_phase, & use material, only: material_phase, &
homogenization_Ngrains, & homogenization_Ngrains, &
phase_constitution, & phase_constitution, &
phase_constitutionInstance phase_constitutionInstance
use mesh, only: mesh_element, & use mesh, only: mesh_element, &
mesh_ipNeighborhood, & mesh_ipNeighborhood, &
FE_NipNeighbors FE_NipNeighbors
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, & use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, &
constitutive_phenopowerlaw_structure constitutive_phenopowerlaw_structure
use constitutive_dislotwin, only: constitutive_dislotwin_label, & use constitutive_dislotwin, only: constitutive_dislotwin_label, &
constitutive_dislotwin_structure constitutive_dislotwin_structure
use constitutive_nonlocal, only: constitutive_nonlocal_label, & use constitutive_nonlocal, only: constitutive_nonlocal_label, &
constitutive_nonlocal_structure constitutive_nonlocal_structure
implicit none implicit none
@ -1407,97 +1402,97 @@ implicit none
!*** output variables ***! !*** output variables ***!
!*** local variables ***! !*** local variables ***!
integer(pInt) e, & ! element index integer(pInt) e, & ! element index
i, & ! integration point index i, & ! integration point index
g, & ! grain index g, & ! grain index
n, & ! neighbor index n, & ! neighbor index
myPhase, & ! phase myPhase, & ! phase
myStructure, & ! lattice structure myStructure, & ! lattice structure
neighboring_e, & ! element index of my neighbor neighboring_e, & ! element index of my neighbor
neighboring_i, & ! integration point index of my neighbor neighboring_i, & ! integration point index of my neighbor
neighboringPhase, & ! phase of my neighbor neighboringPhase, & ! phase of my neighbor
neighboringStructure, & ! lattice structure of my neighbor neighboringStructure, & ! lattice structure of my neighbor
symmetryType ! type of crystal symmetry symmetryType ! type of crystal symmetry
real(pReal), dimension(3,3) :: U, R, & ! polar decomposition of Fe real(pReal), dimension(3,3) :: U, R, & ! polar decomposition of Fe
netRotation ! net rotation between two orientations netRotation ! net rotation between two orientations
logical error logical error
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e)) do g = 1,homogenization_Ngrains(mesh_element(3,e))
! calculate orientation in terms of rotation matrix and euler angles ! calculate orientation in terms of rotation matrix and euler angles
call math_pDecomposition(crystallite_Fe(:,:,g,i,e), U, R, error) ! polar decomposition of Fe call math_pDecomposition(crystallite_Fe(:,:,g,i,e), U, R, error) ! polar decomposition of Fe
if (error) then if (error) then
call IO_warning(650, e, i, g) call IO_warning(650, e, i, g)
crystallite_R(:,:,g,i,e) = 0.0_pReal crystallite_R(:,:,g,i,e) = 0.0_pReal
crystallite_eulerangles(:,g,i,e) = (/400.0, 400.0, 400.0/) ! fake orientation crystallite_eulerangles(:,g,i,e) = (/400.0, 400.0, 400.0/) ! fake orientation
else else
crystallite_R(:,:,g,i,e) = transpose(R) crystallite_R(:,:,g,i,e) = transpose(R)
crystallite_eulerangles(:,g,i,e) = math_RtoEuler(crystallite_R(:,:,g,i,e)) * inDeg crystallite_eulerangles(:,g,i,e) = math_RtoEuler(crystallite_R(:,:,g,i,e)) * inDeg
endif endif
enddo enddo
enddo enddo
enddo enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
if (homogenization_Ngrains(mesh_element(3,e)) == 1_pInt) then ! if single grain ip if (homogenization_Ngrains(mesh_element(3,e)) == 1_pInt) then ! if single grain ip
myPhase = material_phase(1,i,e) ! get my crystal structure myPhase = material_phase(1,i,e) ! get my crystal structure
select case (phase_constitution(myPhase)) select case (phase_constitution(myPhase))
case (constitutive_phenopowerlaw_label) case (constitutive_phenopowerlaw_label)
myStructure = constitutive_phenopowerlaw_structure(phase_constitutionInstance(myPhase)) myStructure = constitutive_phenopowerlaw_structure(phase_constitutionInstance(myPhase))
case (constitutive_dislotwin_label) case (constitutive_dislotwin_label)
myStructure = constitutive_dislotwin_structure(phase_constitutionInstance(myPhase)) myStructure = constitutive_dislotwin_structure(phase_constitutionInstance(myPhase))
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
myStructure = constitutive_nonlocal_structure(phase_constitutionInstance(myPhase)) myStructure = constitutive_nonlocal_structure(phase_constitutionInstance(myPhase))
case default case default
myStructure = '' myStructure = ''
end select end select
do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors
neighboring_e = mesh_ipNeighborhood(1,n,i,e) neighboring_e = mesh_ipNeighborhood(1,n,i,e)
neighboring_i = mesh_ipNeighborhood(2,n,i,e) neighboring_i = mesh_ipNeighborhood(2,n,i,e)
if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists
neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's crystal structure neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's crystal structure
if (myPhase == neighboringPhase) then ! if my neighbor has same phase like me if (myPhase == neighboringPhase) then ! if my neighbor has same phase like me
select case (myStructure) ! get type of symmetry select case (myStructure) ! get type of symmetry
case (1_pInt, 2_pInt) ! fcc and bcc: case (1_pInt, 2_pInt) ! fcc and bcc:
symmetryType = 1_pInt ! -> cubic symmetry symmetryType = 1_pInt ! -> cubic symmetry
case (3_pInt) ! hex: case (3_pInt) ! hex:
symmetryType = 2_pInt ! -> hexagonal symmetry symmetryType = 2_pInt ! -> hexagonal symmetry
case default case default
symmetryType = 0_pInt symmetryType = 0_pInt
end select end select
call math_misorientation( crystallite_misorientation(1:3,n,1,i,e), & call math_misorientation( crystallite_misorientation(1:3,n,1,i,e), &
crystallite_misorientation(4,n,1,i,e), & crystallite_misorientation(4,n,1,i,e), &
netRotation, & netRotation, &
crystallite_R(:,:,1,i,e), & crystallite_R(:,:,1,i,e), &
crystallite_R(:,:,1,neighboring_i,neighboring_e), & crystallite_R(:,:,1,neighboring_i,neighboring_e), &
symmetryType) ! calculate misorientation symmetryType) ! calculate misorientation
else ! for neighbor with different phase else ! for neighbor with different phase
crystallite_misorientation(4,n,1,i,e) = 400.0_pReal ! set misorientation angle to 400 crystallite_misorientation(4,n,1,i,e) = 400.0_pReal ! set misorientation angle to 400
endif endif
else ! no existing neighbor else ! no existing neighbor
crystallite_misorientation(4,n,1,i,e) = 0.0_pReal ! set misorientation angle to zero crystallite_misorientation(4,n,1,i,e) = 0.0_pReal ! set misorientation angle to zero
endif endif
enddo enddo
endif endif
enddo enddo
enddo enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
endsubroutine endsubroutine