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

numpy.MaskedArray behavior

This commit is contained in:
Martin Diehl 2021-07-19 22:56:34 +02:00
parent 1c1dc9383e
commit 03b7532cc5
5 changed files with 36 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
src/grid/DAMASK_grid.f90
View File

@ -231,14 +231,14 @@ program DAMASK_grid
endif endif
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if (loadCases(l)%deformation%mask(i,j)) then if (loadCases(l)%deformation%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j)
else
write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' write(IO_STDOUT,'(2x,12a)',advance='no') ' x '
else
write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j)
endif endif
enddo; write(IO_STDOUT,'(/)',advance='no') enddo; write(IO_STDOUT,'(/)',advance='no')
enddo enddo
if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831 if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831
if (any(loadCases(l)%stress%mask .and. transpose(loadCases(l)%stress%mask) .and. (math_I3<1))) & if (any(.not.(loadCases(l)%stress%mask .or. transpose(loadCases(l)%stress%mask)) .and. (math_I3<1))) &
errorID = 838 ! no rotation is allowed by stress BC errorID = 838 ! no rotation is allowed by stress BC
if (loadCases(l)%stress%myType == 'P') print*, ' P / MPa:' if (loadCases(l)%stress%myType == 'P') print*, ' P / MPa:'
@ -247,9 +247,9 @@ program DAMASK_grid
if (loadCases(l)%stress%myType /= '') then if (loadCases(l)%stress%myType /= '') then
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if (loadCases(l)%stress%mask(i,j)) then if (loadCases(l)%stress%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pReal
else
write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' write(IO_STDOUT,'(2x,12a)',advance='no') ' x '
else
write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pReal
endif endif
enddo; write(IO_STDOUT,'(/)',advance='no') enddo; write(IO_STDOUT,'(/)',advance='no')
enddo enddo
@ -484,15 +484,15 @@ subroutine getMaskedTensor(values,mask,tensor)
values = 0.0 values = 0.0
if (tensor%length == 9) then ! temporary support for deprecated 1D tensor if (tensor%length == 9) then ! temporary support for deprecated 1D tensor
do i = 1,9 do i = 1,9
mask((i-1)/3+1,mod(i-1,3)+1) = tensor%get_asString(i) /= 'x' mask((i-1)/3+1,mod(i-1,3)+1) = tensor%get_asString(i) == 'x'
if (mask((i-1)/3+1,mod(i-1,3)+1)) values((i-1)/3+1,mod(i-1,3)+1) = tensor%get_asFloat(i) if (.not. mask((i-1)/3+1,mod(i-1,3)+1)) values((i-1)/3+1,mod(i-1,3)+1) = tensor%get_asFloat(i)
enddo enddo
else else
do i = 1,3 do i = 1,3
row => tensor%get(i) row => tensor%get(i)
do j = 1,3 do j = 1,3
mask(i,j) = row%get_asString(j) /= 'x' ! ToDo change to np.masked behavior mask(i,j) = row%get_asString(j) == 'x'
if (mask(i,j)) values(i,j) = row%get_asFloat(j) if (.not. mask(i,j)) values(i,j) = row%get_asFloat(j)
enddo enddo
enddo enddo
endif endif

16
src/grid/grid_mech_FEM.f90
View File

@ -324,7 +324,7 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution)
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_FEM_solution end function grid_mechanical_FEM_solution
@ -363,20 +363,20 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
else else
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components F_aimDot = merge(merge(.0_pReal,(F_aim-F_aim_lastInc)/Delta_t_old,stress_BC%mask),.0_pReal,guess) ! estimate deformation rate for prescribed stress components
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
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) + merge(.0_pReal,matmul(deformation_BC%values, F_aim_lastInc),deformation_BC%mask)
elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask) + merge(.0_pReal,deformation_BC%values,deformation_BC%mask)
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask) + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask)
endif endif
if (guess) then if (guess) then
@ -397,9 +397,9 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * Delta_t F_aim = F_aim_lastInc + F_aimDot * Delta_t
if (stress_BC%myType=='P') P_aim = P_aim & if (stress_BC%myType=='P') P_aim = P_aim &
+ merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t
if (stress_BC%myType=='dot_P') P_aim = P_aim & if (stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t
call VecAXPY(solution_current,Delta_t,solution_rate,ierr); CHKERRQ(ierr) call VecAXPY(solution_current,Delta_t,solution_rate,ierr); CHKERRQ(ierr)
@ -574,7 +574,7 @@ subroutine formResidual(da_local,x_local, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask)))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual

16
src/grid/grid_mech_spectral_basic.f90
View File

@ -277,7 +277,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution)
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_spectral_basic_solution end function grid_mechanical_spectral_basic_solution
@ -314,20 +314,20 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components F_aimDot = merge(merge(.0_pReal,(F_aim-F_aim_lastInc)/Delta_t_old,stress_BC%mask),.0_pReal,guess) ! estimate deformation rate for prescribed stress components
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
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) + merge(.0_pReal,matmul(deformation_BC%values, F_aim_lastInc),deformation_BC%mask)
elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask) + merge(.0_pReal,deformation_BC%values,deformation_BC%mask)
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask) + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask)
endif endif
Fdot = utilities_calculateRate(guess, & Fdot = utilities_calculateRate(guess, &
@ -342,9 +342,9 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * Delta_t F_aim = F_aim_lastInc + F_aimDot * Delta_t
if (stress_BC%myType=='P') P_aim = P_aim & if (stress_BC%myType=='P') P_aim = P_aim &
+ merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t
if (stress_BC%myType=='dot_P') P_aim = P_aim & if (stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
@ -499,7 +499,7 @@ subroutine formResidual(in, F, &
! stress BC handling ! stress BC handling
deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
F_aim = F_aim - deltaF_aim F_aim = F_aim - deltaF_aim
err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask)))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme ! updated deformation gradient using fix point algorithm of basic scheme

18
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -309,7 +309,7 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_spectral_polarisation_solution end function grid_mechanical_spectral_polarisation_solution
@ -350,20 +350,20 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components F_aimDot = merge(merge(.0_pReal,(F_aim-F_aim_lastInc)/Delta_t_old,stress_BC%mask),.0_pReal,guess) ! estimate deformation rate for prescribed stress components
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
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) + merge(.0_pReal,matmul(deformation_BC%values, F_aim_lastInc),deformation_BC%mask)
elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask) + merge(.0_pReal,deformation_BC%values,deformation_BC%mask)
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask) + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask)
endif endif
Fdot = utilities_calculateRate(guess, & Fdot = utilities_calculateRate(guess, &
@ -382,9 +382,9 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * Delta_t F_aim = F_aim_lastInc + F_aimDot * Delta_t
if(stress_BC%myType=='P') P_aim = P_aim & if(stress_BC%myType=='P') P_aim = P_aim &
+ merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t
if(stress_BC%myType=='dot_P') P_aim = P_aim & if(stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),& rotation_BC%rotate(F_aim,active=.true.)),&
@ -598,8 +598,8 @@ subroutine formResidual(in, FandF_tau, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
err_BC = maxval(abs(merge(P_av-P_aim, & err_BC = maxval(abs(merge(math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)), &
math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)),& P_av-P_aim, &
params%stress_mask))) params%stress_mask)))
! calculate divergence ! calculate divergence
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal

4
src/grid/spectral_utilities.f90
View File

@ -88,7 +88,7 @@ module spectral_utilities
type, public :: tBoundaryCondition !< set of parameters defining a boundary condition type, public :: tBoundaryCondition !< set of parameters defining a boundary condition
real(pReal), dimension(3,3) :: values = 0.0_pReal real(pReal), dimension(3,3) :: values = 0.0_pReal
logical, dimension(3,3) :: mask = .false. logical, dimension(3,3) :: mask = .true.
character(len=:), allocatable :: myType character(len=:), allocatable :: myType
end type tBoundaryCondition end type tBoundaryCondition
@ -684,7 +684,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
logical :: errmatinv logical :: errmatinv
character(len=pStringLen):: formatString character(len=pStringLen):: formatString
mask_stressVector = reshape(transpose(mask_stress), [9]) mask_stressVector = .not. reshape(transpose(mask_stress), [9])
size_reduced = count(mask_stressVector) size_reduced = count(mask_stressVector)
if(size_reduced > 0) then if(size_reduced > 0) then
temp99_real = math_3333to99(rot_BC%rotate(C)) temp99_real = math_3333to99(rot_BC%rotate(C))