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

Merge branch 'development' into 12-fixOrientationSampling

This commit is contained in:
Martin Diehl 2018-03-13 04:12:03 +01:00
parent 2ceb3aeea0 03c643c433
commit 302a5017d4
11 changed files with 68 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
LICENSE
View File

@ -1,4 +1,4 @@
Copyright 2011-17 Max-Planck-Institut für Eisenforschung GmbH Copyright 2011-18 Max-Planck-Institut für Eisenforschung GmbH
DAMASK is free software: you can redistribute it and/or modify DAMASK is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by

2
VERSION
View File

@ -1 +1 @@
v2.0.1-1082-g8d81e12 v2.0.1-1104-g88c46d5

2
processing/post/addCalculation.py
View File

@ -2,7 +2,7 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,re,sys,collections import os,re,sys,collections
import math # noqa import math,scipy,scipy.linalg # noqa
import numpy as np import numpy as np
from optparse import OptionParser from optparse import OptionParser
import damask import damask

0
src/DAMASK_spectral.f90 Executable file → Normal file
View File

0
src/crystallite.f90 Executable file → Normal file
View File

21
src/kinematics_thermal_expansion.f90
View File

@ -174,8 +174,12 @@ pure function kinematics_thermal_expansion_initialStrain(ipc, ip, el)
offset = thermalMapping(homog)%p(ip,el) offset = thermalMapping(homog)%p(ip,el)
kinematics_thermal_expansion_initialStrain = & kinematics_thermal_expansion_initialStrain = &
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase)) * & (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * &
lattice_thermalExpansion33(1:3,1:3,phase) lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * &
lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * &
lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient
end function kinematics_thermal_expansion_initialStrain end function kinematics_thermal_expansion_initialStrain
@ -215,9 +219,16 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc,
TDot = temperatureRate(homog)%p(offset) TDot = temperatureRate(homog)%p(offset)
TRef = lattice_referenceTemperature(phase) TRef = lattice_referenceTemperature(phase)
Li = TDot* & Li = TDot * ( &
lattice_thermalExpansion33(1:3,1:3,phase)/ & lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient
(1.0_pReal + lattice_thermalExpansion33(1:3,1:3,phase)*(T - TRef)) + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient
) / &
(1.0_pReal \
+ lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. &
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. &
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. &
)
dLi_dTstar3333 = 0.0_pReal dLi_dTstar3333 = 0.0_pReal
end subroutine kinematics_thermal_expansion_LiAndItsTangent end subroutine kinematics_thermal_expansion_LiAndItsTangent

49
src/lattice.f90
View File

@ -1082,9 +1082,10 @@ module lattice
lattice_nu, & lattice_nu, &
lattice_trans_mu, & lattice_trans_mu, &
lattice_trans_nu lattice_trans_nu
real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & ! with higher-order parameters (e.g. temperature-dependent)
lattice_thermalExpansion33
real(pReal), dimension(:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_thermalConductivity33, & lattice_thermalConductivity33, &
lattice_thermalExpansion33, &
lattice_damageDiffusion33, & lattice_damageDiffusion33, &
lattice_vacancyfluxDiffusion33, & lattice_vacancyfluxDiffusion33, &
lattice_vacancyfluxMobility33, & lattice_vacancyfluxMobility33, &
@ -1380,8 +1381,8 @@ subroutine lattice_init
allocate(lattice_C3333(3,3,3,3,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_C66(6,6,Nphases), source=0.0_pReal)
allocate(lattice_trans_C3333(3,3,3,3,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) allocate(lattice_thermalConductivity33 (3,3,Nphases), source=0.0_pReal)
allocate(lattice_thermalExpansion33 (3,3,Nphases), source=0.0_pReal)
allocate(lattice_damageDiffusion33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_damageDiffusion33 (3,3,Nphases), source=0.0_pReal)
allocate(lattice_vacancyfluxDiffusion33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_vacancyfluxDiffusion33 (3,3,Nphases), source=0.0_pReal)
allocate(lattice_vacancyfluxMobility33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_vacancyfluxMobility33 (3,3,Nphases), source=0.0_pReal)
@ -1545,11 +1546,17 @@ subroutine lattice_init
case ('thermal_conductivity33') case ('thermal_conductivity33')
lattice_thermalConductivity33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) lattice_thermalConductivity33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt)
case ('thermal_expansion11') case ('thermal_expansion11')
lattice_thermalExpansion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T)
lattice_thermalExpansion33(1,1,i-1_pInt,section) = IO_floatValue(line,chunkPos,i)
enddo
case ('thermal_expansion22') case ('thermal_expansion22')
lattice_thermalExpansion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T)
lattice_thermalExpansion33(2,2,i-1_pInt,section) = IO_floatValue(line,chunkPos,i)
enddo
case ('thermal_expansion33') case ('thermal_expansion33')
lattice_thermalExpansion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T)
lattice_thermalExpansion33(3,3,i-1_pInt,section) = IO_floatValue(line,chunkPos,i)
enddo
case ('specific_heat') case ('specific_heat')
lattice_specificHeat(section) = IO_floatValue(line,chunkPos,2_pInt) lattice_specificHeat(section) = IO_floatValue(line,chunkPos,2_pInt)
case ('vacancyformationenergy') case ('vacancyformationenergy')
@ -1756,22 +1763,24 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc)
end select end select
end select end select
lattice_thermalConductivity33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& forall (i = 1_pInt:3_pInt) &
lattice_thermalConductivity33(1:3,1:3,myPhase)) lattice_thermalExpansion33 (1:3,1:3,i,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_thermalExpansion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_thermalExpansion33 (1:3,1:3,i,myPhase))
lattice_thermalExpansion33(1:3,1:3,myPhase))
lattice_DamageDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_thermalConductivity33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_DamageDiffusion33(1:3,1:3,myPhase)) lattice_thermalConductivity33 (1:3,1:3,myPhase))
lattice_vacancyfluxDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_DamageDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_vacancyfluxDiffusion33(1:3,1:3,myPhase)) lattice_DamageDiffusion33 (1:3,1:3,myPhase))
lattice_vacancyfluxMobility33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_vacancyfluxDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_vacancyfluxMobility33(1:3,1:3,myPhase)) lattice_vacancyfluxDiffusion33 (1:3,1:3,myPhase))
lattice_PorosityDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_vacancyfluxMobility33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_PorosityDiffusion33(1:3,1:3,myPhase)) lattice_vacancyfluxMobility33 (1:3,1:3,myPhase))
lattice_PorosityDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_PorosityDiffusion33 (1:3,1:3,myPhase))
lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase)) lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase))
lattice_hydrogenfluxMobility33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_hydrogenfluxMobility33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_hydrogenfluxMobility33(1:3,1:3,myPhase)) lattice_hydrogenfluxMobility33 (1:3,1:3,myPhase))
select case(lattice_structure(myPhase)) select case(lattice_structure(myPhase))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

10
src/math.f90
View File

@ -2566,8 +2566,7 @@ pure function math_rotate_forward33(tensor,rot_tensor)
real(pReal), dimension(3,3) :: math_rotate_forward33 real(pReal), dimension(3,3) :: math_rotate_forward33
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
math_rotate_forward33 = math_mul33x33(rot_tensor,& math_rotate_forward33 = math_mul33x33(rot_tensor,math_mul33x33(tensor,transpose(rot_tensor)))
math_mul33x33(tensor,math_transpose33(rot_tensor)))
end function math_rotate_forward33 end function math_rotate_forward33
@ -2581,8 +2580,7 @@ pure function math_rotate_backward33(tensor,rot_tensor)
real(pReal), dimension(3,3) :: math_rotate_backward33 real(pReal), dimension(3,3) :: math_rotate_backward33
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
math_rotate_backward33 = math_mul33x33(math_transpose33(rot_tensor),& math_rotate_backward33 = math_mul33x33(transpose(rot_tensor),math_mul33x33(tensor,rot_tensor))
math_mul33x33(tensor,rot_tensor))
end function math_rotate_backward33 end function math_rotate_backward33
@ -2603,8 +2601,8 @@ pure function math_rotate_forward3333(tensor,rot_tensor)
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt; do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt; do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt
do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt; do o = 1_pInt,3_pInt; do p = 1_pInt,3_pInt do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt; do o = 1_pInt,3_pInt; do p = 1_pInt,3_pInt
math_rotate_forward3333(i,j,k,l) = math_rotate_forward3333(i,j,k,l) & math_rotate_forward3333(i,j,k,l) = math_rotate_forward3333(i,j,k,l) &
+ rot_tensor(m,i) * rot_tensor(n,j) & + rot_tensor(i,m) * rot_tensor(j,n) &
* rot_tensor(o,k) * rot_tensor(p,l) * tensor(m,n,o,p) * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
end function math_rotate_forward3333 end function math_rotate_forward3333

14
src/spectral_mech_Basic.f90
View File

@ -301,7 +301,6 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
grid3 grid3
use math, only: & use math, only: &
math_rotate_backward33, & math_rotate_backward33, &
math_transpose33, &
math_mul3333xx33 math_mul3333xx33
use debug, only: & use debug, only: &
debug_level, & debug_level, &
@ -349,9 +348,9 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', math_transpose33(F_aim) ' deformation gradient aim =', transpose(F_aim)
flush(6) flush(6)
endif newIteration endif newIteration
@ -516,13 +515,10 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg C_minMaxAvgLastInc = C_minMaxAvg
if (guess) then ! QUESTION: better with a = L ? x:y F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aimDot = stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old ! initialize with correction based on last inc
else
F_aimDot = 0.0_pReal
endif ! components of deformation_BC%maskFloat always start out with zero
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F

20
src/spectral_mech_Polarisation.f90
View File

@ -333,7 +333,6 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
IO_intOut IO_intOut
use math, only: & use math, only: &
math_rotate_backward33, & math_rotate_backward33, &
math_transpose33, &
math_mul3333xx33, & math_mul3333xx33, &
math_invSym3333, & math_invSym3333, &
math_mul33x33 math_mul33x33
@ -402,9 +401,9 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', math_transpose33(F_aim) ' deformation gradient aim =', transpose(F_aim)
flush(6) flush(6)
endif newIteration endif newIteration
@ -548,7 +547,6 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_mul3333xx33, & math_mul3333xx33, &
math_transpose33, &
math_rotate_backward33 math_rotate_backward33
use numerics, only: & use numerics, only: &
worldrank worldrank
@ -630,11 +628,7 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg C_minMaxAvgLastInc = C_minMaxAvg
if (guess) then ! QUESTION: better with a = L ? x:y F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aimDot = stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old ! initialize with correction based on last inc
else
F_aimDot = 0.0_pReal
endif
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
@ -676,16 +670,16 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
math_mul3333xx33(C_scale,& math_mul3333xx33(C_scale,&
math_mul33x33(math_transpose33(F_lambda33),& math_mul33x33(transpose(F_lambda33),&
F_lambda33)-math_I3))*0.5_pReal)& F_lambda33)-math_I3))*0.5_pReal)&
+ math_I3 + math_I3
F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
endif endif
nullify(F) nullify(F)
nullify(F_tau) nullify(F_tau)
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine Polarisation_forward end subroutine Polarisation_forward

10
src/spectral_utilities.f90 Executable file → Normal file
View File

@ -70,7 +70,6 @@ module spectral_utilities
! derived types ! derived types
type, public :: tSolutionState !< return type of solution from spectral solver variants type, public :: tSolutionState !< return type of solution from spectral solver variants
logical :: converged = .true. logical :: converged = .true.
logical :: regrid = .false.
logical :: stagConverged = .true. logical :: stagConverged = .true.
logical :: termIll = .false. logical :: termIll = .false.
integer(pInt) :: iterationsNeeded = 0_pInt integer(pInt) :: iterationsNeeded = 0_pInt
@ -782,7 +781,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
if(debugGeneral) then if(debugGeneral) then
write(6,'(/,a)') ' ... updating masked compliance ............................................' write(6,'(/,a)') ' ... updating masked compliance ............................................'
write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',&
transpose(temp99_Real)/1.e9_pReal transpose(temp99_Real)*1.0e-9_pReal
flush(6) flush(6)
endif endif
k = 0_pInt ! calculate reduced stiffness k = 0_pInt ! calculate reduced stiffness
@ -837,7 +836,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
endif endif
if(debugGeneral) then if(debugGeneral) then
write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') & write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') &
' Masked Compliance (load) / GPa =', transpose(temp99_Real*1.e-9_pReal) ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal
flush(6) flush(6)
endif endif
utilities_maskedCompliance = math_Plain99to3333(temp99_Real) utilities_maskedCompliance = math_Plain99to3333(temp99_Real)
@ -935,7 +934,6 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
debug_reset, & debug_reset, &
debug_info debug_info
use math, only: & use math, only: &
math_transpose33, &
math_rotate_forward33, & math_rotate_forward33, &
math_det33 math_det33
use mesh, only: & use mesh, only: &
@ -994,10 +992,10 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if (debugRotation) & if (debugRotation) &
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
math_transpose33(P_av)*1.e-6_pReal transpose(P_av)*1.e-6_pReal
P_av = math_rotate_forward33(P_av,rotation_BC) P_av = math_rotate_forward33(P_av,rotation_BC)
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
math_transpose33(P_av)*1.e-6_pReal transpose(P_av)*1.e-6_pReal
flush(6) flush(6)
max_dPdF = 0.0_pReal max_dPdF = 0.0_pReal