From 0a7d4f61acfce1de1c7deb2b7962802641b0163b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Sep 2020 16:24:08 +0200 Subject: [PATCH] Need only logical mask 'merge' substitutes multiplication with float mask --- src/grid/DAMASK_grid.f90 | 105 +++++++++---------- src/grid/grid_mech_FEM.f90 | 26 +++-- src/grid/grid_mech_spectral_basic.f90 | 26 +++-- src/grid/grid_mech_spectral_polarisation.f90 | 30 +++--- src/grid/spectral_utilities.f90 | 15 ++- 5 files changed, 96 insertions(+), 106 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 54d080744..d33d136f7 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -24,7 +24,7 @@ program DAMASK_grid use grid_damage_spectral use grid_thermal_spectral use results - + implicit none !-------------------------------------------------------------------------------------------------- @@ -88,7 +88,7 @@ program DAMASK_grid mech_updateCoords procedure(grid_mech_spectral_basic_restartWrite), pointer :: & mech_restartWrite - + external :: & quit class (tNode), pointer :: & @@ -97,10 +97,10 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) - - call CPFEM_initAll + + call CPFEM_initAll print'(/,a)', ' <<<+- DAMASK_spectral init -+>>>'; flush(6) - + print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' @@ -123,16 +123,16 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! assign mechanics solver depending on selected type - + debug_grid => config_debug%get('grid',defaultVal=emptyList) - select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) + select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) case ('Basic') mech_init => grid_mech_spectral_basic_init mech_forward => grid_mech_spectral_basic_forward mech_solution => grid_mech_spectral_basic_solution mech_updateCoords => grid_mech_spectral_basic_updateCoords mech_restartWrite => grid_mech_spectral_basic_restartWrite - + case ('Polarisation') if(debug_grid%contains('basic')) & call IO_warning(42, ext_msg='debug Divergence') @@ -141,7 +141,7 @@ program DAMASK_grid mech_solution => grid_mech_spectral_polarisation_solution mech_updateCoords => grid_mech_spectral_polarisation_updateCoords mech_restartWrite => grid_mech_spectral_polarisation_restartWrite - + case ('FEM') if(debug_grid%contains('basic')) & call IO_warning(42, ext_msg='debug Divergence') @@ -150,24 +150,24 @@ program DAMASK_grid mech_solution => grid_mech_FEM_solution mech_updateCoords => grid_mech_FEM_updateCoords mech_restartWrite => grid_mech_FEM_restartWrite - + case default call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver'))) - + end select - + !-------------------------------------------------------------------------------------------------- -! reading information from load case file and to sanity checks +! reading information from load case file and to sanity checks fileContent = IO_readlines(trim(interface_loadFile)) if(size(fileContent) == 0) call IO_error(307,ext_msg='No load case specified') - + allocate (loadCases(0)) ! array of load cases do currentLoadCase = 1, size(fileContent) line = fileContent(currentLoadCase) if (IO_isBlank(line)) cycle chunkPos = IO_stringPos(line) - + do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase select case (IO_lc(IO_stringValue(line,chunkPos,i))) case('l','fdot','dotf','f') @@ -180,7 +180,7 @@ program DAMASK_grid enddo if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile)) ! error message for incomplete loadcase - + newLoadCase%stress%myType='p' field = 1 newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default @@ -192,7 +192,7 @@ program DAMASK_grid field = field + 1 newLoadCase%ID(field) = FIELD_DAMAGE_ID endif damageActive - + call newLoadCase%rot%fromEulers(real([0.0,0.0,0.0],pReal)) readIn: do i = 1, chunkPos(1) select case (IO_lc(IO_stringValue(line,chunkPos,i))) @@ -210,18 +210,16 @@ program DAMASK_grid temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a * if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo - newLoadCase%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) ! logical mask in 3x3 notation - newLoadCase%deformation%maskFloat = merge(ones,zeros,newLoadCase%deformation%maskLogical) ! float (1.0/0.0) mask in 3x3 notation - newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation + newLoadCase%deformation%mask = transpose(reshape(temp_maskVector,[ 3,3])) ! mask in 3x3 notation + newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation case('p','stress', 's') temp_valueVector = 0.0_pReal do j = 1, 9 temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo - newLoadCase%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) - newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical) - newLoadCase%stress%values = math_9to33(temp_valueVector) + newLoadCase%stress%mask = transpose(reshape(temp_maskVector,[ 3,3])) + newLoadCase%stress%values = math_9to33(temp_valueVector) case('t','time','delta') ! increment time newLoadCase%time = IO_floatValue(line,chunkPos,i+1) case('n','incs','increments') ! number of increments @@ -258,9 +256,9 @@ program DAMASK_grid call newLoadCase%rot%fromMatrix(math_9to33(temp_valueVector)) end select enddo readIn - + newLoadCase%followFormerTrajectory = merge(.true.,.false.,currentLoadCase > 1) ! by default, guess from previous load case - + reportAndCheck: if (worldrank == 0) then write (loadcase_string, '(i0)' ) currentLoadCase print'(/,a,i0)', ' load case: ', currentLoadCase @@ -268,8 +266,8 @@ program DAMASK_grid print*, ' drop guessing along trajectory' if (newLoadCase%deformation%myType == 'l') then do j = 1, 3 - if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & - any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined + if (any(newLoadCase%deformation%mask(j,1:3) .eqv. .true.) .and. & + any(newLoadCase%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined enddo print*, ' velocity gradient:' else if (newLoadCase%deformation%myType == 'f') then @@ -278,20 +276,19 @@ program DAMASK_grid print*, ' deformation gradient rate:' endif do i = 1, 3; do j = 1, 3 - if(newLoadCase%deformation%maskLogical(i,j)) then + if(newLoadCase%deformation%mask(i,j)) then write(6,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j) else write(6,'(2x,12a)',advance='no') ' * ' endif enddo; write(6,'(/)',advance='no') enddo - if (any(newLoadCase%stress%maskLogical .eqv. & - newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only - if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) & - .and. (math_I3<1))) errorID = 838 ! no rotation is allowed by stress BC + if (any(newLoadCase%stress%mask .eqv. newLoadCase%deformation%mask)) errorID = 831 ! exclusive or masking only + if (any(newLoadCase%stress%mask .and. transpose(newLoadCase%stress%mask) .and. (math_I3<1))) & + errorID = 838 ! no rotation is allowed by stress BC print*, ' stress / GPa:' do i = 1, 3; do j = 1, 3 - if(newLoadCase%stress%maskLogical(i,j)) then + if(newLoadCase%stress%mask(i,j)) then write(6,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal else write(6,'(2x,12a)',advance='no') ' * ' @@ -325,13 +322,13 @@ program DAMASK_grid select case (loadCases(1)%ID(field)) case(FIELD_MECH_ID) call mech_init - + case(FIELD_THERMAL_ID) call grid_thermal_spectral_init - + case(FIELD_DAMAGE_ID) call grid_damage_spectral_init - + end select enddo @@ -348,16 +345,16 @@ program DAMASK_grid '.sta',form='FORMATTED', position='APPEND', status='OLD') endif writeHeader endif - + writeUndeformed: if (interface_restartInc < 1) then print'(/,a)', ' ... writing initial configuration to file ........................' call CPFEM_results(0,0.0_pReal) endif writeUndeformed - + loadCaseLooping: do currentLoadCase = 1, size(loadCases) time0 = time ! load case start time guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc - + incLooping: do inc = 1, loadCases(currentLoadCase)%incs totalIncsCounter = totalIncsCounter + 1 @@ -382,13 +379,13 @@ program DAMASK_grid endif endif timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step - + skipping: if (totalIncsCounter <= interface_restartInc) then ! not yet at restart inc? time = time + timeinc ! just advance time, skip already performed calculation guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference else skipping stepFraction = 0 ! fraction scaled by stepFactor**cutLevel - + subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time time = time + timeinc ! forward target time @@ -417,7 +414,7 @@ program DAMASK_grid deformation_BC = loadCases(currentLoadCase)%deformation, & stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rot) - + case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) end select @@ -438,9 +435,9 @@ program DAMASK_grid case(FIELD_DAMAGE_ID) solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld) end select - + if (.not. solres(field)%converged) exit ! no solution found - + enddo stagIter = stagIter + 1 stagIterate = stagIter < stagItMax & @@ -474,17 +471,17 @@ program DAMASK_grid if (worldrank == 0) close(statUnit) call quit(0) ! quit endif - + enddo subStepLooping - + cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc - + if (all(solres(:)%converged)) then print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged' else print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged' endif; flush(6) - + if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency print'(1/,a)', ' ... writing results to file ......................................' flush(6) @@ -495,17 +492,17 @@ program DAMASK_grid call CPFEM_restartWrite endif endif skipping - + enddo incLooping - + enddo loadCaseLooping - - + + !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation print'(/,a)', ' ###########################################################################' if (worldrank == 0) close(statUnit) - + call quit(0) ! no complains ;) - + end program DAMASK_grid diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 9efd88505..9d9c962e1 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -25,8 +25,6 @@ module grid_mech_FEM implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types type(tSolutionParams) :: params type :: tNumerics @@ -282,7 +280,7 @@ function grid_mech_FEM_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -326,7 +324,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, !-------------------------------------------------------------------------------------------------- ! set module wide available data - params%stress_mask = stress_BC%maskFloat + params%stress_mask = stress_BC%mask params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -340,20 +338,20 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, else C_volAvgLastInc = C_volAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + F_aimDot = F_aimDot & + + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + F_aimDot = F_aimDot & + + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif if (guess) then @@ -374,9 +372,9 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc if (stress_BC%myType=='p') then - P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) elseif (stress_BC%myType=='pdot') then !UNTESTED - P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) endif call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) @@ -543,7 +541,7 @@ subroutine formResidual(da_local,x_local, & !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc - err_BC = maxval(abs(params%stress_mask * (P_av - P_aim))) ! mask = 0.0 when no stress bc + err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index cd795f005..fd5e63663 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -24,8 +24,6 @@ module grid_mech_spectral_basic implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types type(tSolutionParams) :: params type :: tNumerics @@ -247,7 +245,7 @@ function grid_mech_spectral_basic_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) !-------------------------------------------------------------------------------------------------- @@ -291,7 +289,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo !-------------------------------------------------------------------------------------------------- ! set module wide available data - params%stress_mask = stress_BC%maskFloat + params%stress_mask = stress_BC%mask params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -305,20 +303,20 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + F_aimDot = F_aimDot & + + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + F_aimDot = F_aimDot & + + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif Fdot = utilities_calculateRate(guess, & @@ -333,9 +331,9 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc if (stress_BC%myType=='p') then - P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) elseif (stress_BC%myType=='pdot') then !UNTESTED - P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) endif F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average @@ -491,7 +489,7 @@ subroutine formResidual(in, F, & ! stress BC handling deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - P_aim))) ! mask = 0.0 when no stress bc + err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 4dc9d6d7c..91632f3b9 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -24,8 +24,6 @@ module grid_mech_spectral_polarisation implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types type(tSolutionParams) :: params type :: tNumerics @@ -275,7 +273,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) if(num%update_gamma) then call utilities_updateGamma(C_minMaxAvg) C_scale = C_minMaxAvg @@ -325,7 +323,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc !-------------------------------------------------------------------------------------------------- ! set module wide available data - params%stress_mask = stress_BC%maskFloat + params%stress_mask = stress_BC%mask params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -341,20 +339,20 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + F_aimDot = F_aimDot & + + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + F_aimDot = F_aimDot & + + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif Fdot = utilities_calculateRate(guess, & @@ -373,9 +371,9 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc if (stress_BC%myType=='p') then - P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) elseif (stress_BC%myType=='pdot') then !UNTESTED - P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) endif F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average @@ -593,9 +591,9 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc - err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & - -params%rotation_BC%rotate(F_av)) + & - params%stress_mask * (P_av-P_aim))) ! mask = 0.0 for no bc + err_BC = maxval(abs(merge(P_av-P_aim, & + math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)),& + params%stress_mask))) ! calculate divergence tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 52df30c39..896337bf6 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -82,10 +82,9 @@ module spectral_utilities end type tSolutionState type, public :: tBoundaryCondition !< set of parameters defining a boundary condition - real(pReal), dimension(3,3) :: values = 0.0_pReal, & - maskFloat = 0.0_pReal - logical, dimension(3,3) :: maskLogical = .false. - character(len=pStringLen) :: myType = 'None' + real(pReal), dimension(3,3) :: values = 0.0_pReal + logical, dimension(3,3) :: mask = .false. + character(len=pStringLen) :: myType = 'None' end type tBoundaryCondition type, public :: tLoadCase @@ -101,11 +100,11 @@ module spectral_utilities integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) end type tLoadCase - type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase - real(pReal), dimension(3,3) :: stress_mask, stress_BC + type, public :: tSolutionParams + real(pReal), dimension(3,3) :: stress_BC + logical, dimension(3,3) :: stress_mask type(rotation) :: rotation_BC - real(pReal) :: timeinc - real(pReal) :: timeincOld + real(pReal) :: timeinc, timeincOld end type tSolutionParams type :: tNumerics