diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 0b842d2ca..621e72fab 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -403,7 +403,7 @@ return endfunction -subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,disorientation,ipc,ip,el) +subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el) !********************************************************************* !* This function calculates from state needed variables * !* INPUT: * @@ -432,7 +432,6 @@ integer(pInt), intent(in) :: ipc,ip,el real(pReal), intent(in) :: Temperature real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp -real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: disorientation select case (phase_constitution(material_phase(ipc,ip,el))) @@ -449,7 +448,7 @@ real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: disorientation call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el) case (constitutive_nonlocal_label) - call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Tstar_v, Fe, Fp, disorientation, ipc, ip, el) + call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Tstar_v, Fe, Fp, ipc, ip, el) end select @@ -509,7 +508,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ip endsubroutine -subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fe, Fp, Temperature, disorientation, subdt, ipc, ip, el) +subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, ipc, ip, el) !********************************************************************* !* This subroutine contains the constitutive equation for * !* calculating the rate of change of microstructure * @@ -548,8 +547,6 @@ implicit none integer(pInt), intent(in) :: ipc, ip, el real(pReal), intent(in) :: Temperature, & subdt -real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: & - disorientation real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & Fe, & Fp @@ -579,7 +576,7 @@ select case (phase_constitution(material_phase(ipc,ip,el))) constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) case (constitutive_nonlocal_label) - call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, disorientation, subdt, & + call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, & constitutive_state, constitutive_subState0, constitutive_relevantState, subdt, ipc, ip, el) end select diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index c5c485990..444ac9dbd 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -81,6 +81,8 @@ real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_ constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_v, & ! dislocation velocity constitutive_nonlocal_rhoDotFlux ! dislocation convection term +real(pReal), dimension(:,:,:,:,:,:), allocatable :: constitutive_nonlocal_compatibility, & ! slip system compatibility between me and my neighbors + constitutive_nonlocal_transmissivity ! transmissivity between me and my neighbors real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance constitutive_nonlocal_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance constitutive_nonlocal_interactionMatrixSlipSlip ! interaction matrix of the different slip systems for each instance @@ -116,7 +118,8 @@ use IO, only: IO_lc, & IO_intValue, & IO_error use mesh, only: mesh_NcpElems, & - mesh_maxNips + mesh_maxNips, & + FE_maxNipNeighbors use material, only: homogenization_maxNgrains, & phase_constitution, & phase_constitutionInstance, & @@ -422,6 +425,12 @@ constitutive_nonlocal_v = 0.0_pReal allocate(constitutive_nonlocal_rhoDotFlux(maxTotalNslip, 8, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems)) constitutive_nonlocal_rhoDotFlux = 0.0_pReal +allocate(constitutive_nonlocal_compatibility(maxTotalNslip, 2, maxTotalNslip, FE_maxNipNeighbors, mesh_maxNips, mesh_NcpElems)) +constitutive_nonlocal_compatibility = 0.0_pReal + +allocate(constitutive_nonlocal_transmissivity(maxTotalNslip, 2, maxTotalNslip, FE_maxNipNeighbors, mesh_maxNips, mesh_NcpElems)) +constitutive_nonlocal_transmissivity = 0.0_pReal + do i = 1,maxNinstance myStructure = constitutive_nonlocal_structure(i) ! lattice structure of this instance @@ -742,7 +751,7 @@ endfunction !********************************************************************* !* calculates quantities characterizing the microstructure * !********************************************************************* -subroutine constitutive_nonlocal_microstructure(state, Temperature, Tstar_v, Fe, Fp, disorientation, g, ip, el) +subroutine constitutive_nonlocal_microstructure(state, Temperature, Tstar_v, Fe, Fp, g, ip, el) use prec, only: pReal, & pInt, & @@ -769,6 +778,7 @@ use mesh, only: mesh_NcpElems, & mesh_ipCenterOfGravity use material, only: homogenization_maxNgrains, & material_phase, & + phase_localConstitution, & phase_constitutionInstance use lattice, only: lattice_Sslip, & lattice_Sslip_v, & @@ -791,8 +801,6 @@ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) Fp ! plastic deformation gradient real(pReal), dimension(6), intent(in) :: & Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation -real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: & - disorientation ! crystal disorientation between me and my neighbor as quaternion !*** input/output variables @@ -825,6 +833,8 @@ real(pReal), dimension(3,3) :: sigma, & ! dislocation stre invPositionDifference ! inverse of a 3x3 matrix containing finite differences of pairs of position vectors real(pReal), dimension(6) :: transmissivity, & ! transmissivity factor for each interface Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress in Mandel notation +real(pReal), dimension(2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & + rhoExcess ! central excess density real(pReal), dimension(6,2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & neighboring_rhoExcess ! excess density for each neighbor, dislo character and slip system real(pReal), dimension(6,3) :: neighboring_position ! position vector of each neighbor when seen from the centreal material point's lattice frame @@ -879,28 +889,49 @@ F = math_mul33x33(Fe(:,:,g,ip,el), Fp(:,:,g,ip,el)) invFe = math_inv3x3(Fe(:,:,g,ip,el)) nu = constitutive_nonlocal_nu(myInstance) +forall (s = 1:ns, c = 1:2) & + rhoExcess(c,s) = state(g,ip,el)%p((2*c-2)*ns+s) + abs(state(g,ip,el)%p((2*c+2)*ns+s)) & + - state(g,ip,el)%p((2*c-1)*ns+s) - abs(state(g,ip,el)%p((2*c+3)*ns+s)) do n = 1,6 neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) - if ( neighboring_ip == 0 .or. constitutive_nonlocal_transmissivity(disorientation(:,n)) < 1.0_pReal ) then - neighboring_el = el + if ( neighboring_ip == 0 .or. neighboring_el == 0 ) then ! at free surfaces ... + neighboring_el = el ! ... use central values instead of neighboring values neighboring_ip = ip - neighboring_F = F + neighboring_position(n,:) = 0.0_pReal + neighboring_rhoExcess(n,:,:) = rhoExcess + elseif (.not. phase_localConstitution(material_phase(1,neighboring_ip,neighboring_el))) then ! for neighbors with local constitution + neighboring_el = el ! ... use central values instead of neighboring values + neighboring_ip = ip + neighboring_position(n,:) = 0.0_pReal + neighboring_rhoExcess(n,:,:) = rhoExcess + elseif (myStructure /= & + constitutive_nonlocal_structure(phase_constitutionInstance(material_phase(1,neighboring_ip,neighboring_el)))) then ! for neighbors with different crystal structure + neighboring_el = el ! ... use central values instead of neighboring values + neighboring_ip = ip + neighboring_position(n,:) = 0.0_pReal + neighboring_rhoExcess(n,:,:) = rhoExcess else - neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el)) + forall (s = 1:ns, c = 1:2) & + neighboring_rhoExcess(n,c,s) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*ns+s) & + + abs(state(g,neighboring_ip,neighboring_el)%p((2*c+2)*ns+s)) & + - state(g,neighboring_ip,neighboring_el)%p((2*c-1)*ns+s) & + - abs(state(g,neighboring_ip,neighboring_el)%p((2*c+3)*ns+s)) + if ( any( neighboring_rhoExcess(n,:,:)*rhoExcess < 0.0_pReal & + .and. abs(neighboring_rhoExcess(n,:,:)) > 1.0_pReal ) ) then ! at grain boundary (=significant change of sign in any excess density) ... + neighboring_el = el ! ... use central values instead of neighboring values + neighboring_ip = ip + neighboring_position(n,:) = 0.0_pReal + neighboring_rhoExcess(n,:,:) = rhoExcess + else + neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el)) + neighboring_position(n,:) = & + 0.5_pReal * math_mul33x3( math_mul33x33(invFe,neighboring_F) + Fp(:,:,g,ip,el), & + mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el) ) + endif endif - neighboring_position(n,:) = & - 0.5_pReal * math_mul33x3( math_mul33x33(invFe,neighboring_F) + Fp(:,:,g,ip,el), & - mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el) ) - - forall (s = 1:ns, c = 1:2) & - neighboring_rhoExcess(n,c,s) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*ns+s) & - + abs(state(g,neighboring_ip,neighboring_el)%p((2*c+2)*ns+s)) & - - state(g,neighboring_ip,neighboring_el)%p((2*c-1)*ns+s) & - - abs(state(g,neighboring_ip,neighboring_el)%p((2*c+3)*ns+s)) - enddo invPositionDifference = math_inv3x3(neighboring_position((/1,3,5/),:) - neighboring_position((/2,4,6/),:)) @@ -1192,7 +1223,7 @@ endsubroutine !********************************************************************* !* rate of change of microstructure * !********************************************************************* -subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, previousTstar_v, Fe, Fp, Temperature, disorientation, dt_previous, & +subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, previousTstar_v, Fe, Fp, Temperature, dt_previous, & state, previousState, relevantState, timestep, g,ip,el) use prec, only: pReal, & @@ -1245,8 +1276,6 @@ real(pReal), intent(in) :: Temperature, & ! temperat dt_previous ! time increment between previous and current state real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation previousTstar_v ! previous 2nd Piola-Kirchhoff stress in Mandel notation -real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: & - disorientation ! crystal disorientation between me and my neighbor as quaternion real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & Fe, & ! elastic deformation gradient Fp ! plastic deformation gradient @@ -1265,6 +1294,7 @@ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), in integer(pInt) myInstance, & ! current instance of this constitution myStructure, & ! current lattice structure ns, & ! short notation for the total number of active slip systems + neighboring_n, & ! neighbor index of myself when seen from my neighbor neighboring_el, & ! element number of my neighbor neighboring_ip, & ! integration point of my neighbor c, & ! character of dislocation @@ -1273,7 +1303,9 @@ integer(pInt) myInstance, & ! current opposite_ip, & ! ip of my opposite neighbor opposite_el, & ! element index of my opposite neighbor t, & ! type of dislocation + topp, & ! type of dislocation with opposite sign to t s, & ! index of my current slip system + s2, & sLattice, & ! index of my current slip system according to lattice order i real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),10) :: & @@ -1318,7 +1350,8 @@ real(pReal), dimension(3) :: surfaceNormal, & ! surface surfaceNormal_currentconf ! surface normal in current configuration real(pReal) area, & ! area of the current interface detFe, & ! determinant of elastic defornmation gradient - transmissivity, & ! transmissivity of interfaces for dislocation flux + transmissivity, & ! transmissivity of dislocation flux for different slip systems in neighboring material points for a specific charcter of dislocations + compatibility, & ! compatibility of different slip systems in neighboring material points for a specific charcter of dislocations lineLength, & ! dislocation line length leaving the current interface D ! self diffusion logical, dimension(3) :: periodicSurfaceFlux ! flag indicating periodic fluxes at surfaces when surface normal points mainly in x, y and z direction respectively (in reference configuration) @@ -1365,7 +1398,7 @@ previousTdislocation_v = previousState(g,ip,el)%p(12*ns+1:12*ns+6) !*** sanity check for timestep if (timestep <= 0.0_pReal) then ! if illegal timestep... - dotState(1,ip,el)%p(1:10*ns) = 0.0_pReal ! ...return without doing anything (-> zero dotState) + dotState(g,ip,el)%p = 0.0_pReal ! ...return without doing anything (-> zero dotState) return endif @@ -1470,6 +1503,14 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el) opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el) + if ( neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt ) then ! if neighbor exists ... + do neighboring_n = 1,FE_NipNeighbors(mesh_element(2,neighboring_el)) ! find neighboring index that points from my neighbor to myself + if ( el == mesh_ipNeighborhood(1,opposite_n,ip,neighboring_el) & ! special case if no neighbor at all... + .and. ip == mesh_ipNeighborhood(2,opposite_n,ip,neighboring_el) ) & + exit ! ...exit without any flux calculation + enddo + endif + if ( neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt ) then ! if neighbor exists, average deformation gradient neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el)) Favg = 0.5_pReal * (F + neighboring_F) @@ -1481,9 +1522,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) surfaceNormal = math_mul33x3(transpose(Fe(:,:,g,ip,el)), surfaceNormal_currentconf) / detFe ! ... and lattice configuration area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal) surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) ! normalize the surface normal to unit length - - transmissivity = constitutive_nonlocal_transmissivity(disorientation(:,n)) - + if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then ! if neighbor exists... if ( .not. phase_localConstitution(material_phase(1,neighboring_ip,neighboring_el))) then ! ... and is of nonlocal constitution... forall (t = 1:4) & ! ... then calculate neighboring flux density @@ -1502,21 +1541,26 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) endif do s = 1,ns - do t = 1,4 + do t = 1,4 + c = (t + 1) / 2 ! dislocation character if ( fluxdensity(s,t) * math_mul3x3(m(:,s,t), surfaceNormal) > 0.0_pReal ) then ! outgoing flux + transmissivity = sum(constitutive_nonlocal_transmissivity(:,c,s,n,ip,el)) ! overall transmissivity between my system s and all neighboring systems s2 for this dislocation character lineLength = fluxdensity(s,t) * math_mul3x3(m(:,s,t), surfaceNormal) * area ! line length that wants to leave this interface rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract positive dislocation flux that leaves the material point rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) & * sign(1.0_pReal, fluxdensity(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point -! if (selectiveDebugger .and. s==1) & -! write(6,'(a22,i2,a15,i2,a3,e20.10)') 'outgoing flux of type ',t,' to neighbor ',n,' : ', & -! -lineLength / mesh_ipVolume(ip,el) else ! incoming flux - lineLength = neighboring_fluxdensity(s,t) * math_mul3x3(m(:,s,t), surfaceNormal) * area ! line length that wants to leave this interface - rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) * transmissivity ! subtract negative dislocation flux that enters the material point -! if (selectiveDebugger .and. s==1) & -! write(6,'(a22,i2,a15,i2,a3,e20.10)') 'incoming flux of type ',t,' from neighbor ',n,' : ', & -! -lineLength / mesh_ipVolume(ip,el) * transmissivity + do s2 = 1,ns + compatibility = constitutive_nonlocal_compatibility(s,c,s2,neighboring_n,neighboring_ip,neighboring_el) ! compatibility of system s2 of my neighbor to system s in my material point + transmissivity = constitutive_nonlocal_transmissivity(s,c,s2,neighboring_n,neighboring_ip,neighboring_el) + lineLength = neighboring_fluxdensity(s2,t) * math_mul3x3(m(:,s,t), surfaceNormal) * area ! line length that wants to leave this interface + if (compatibility > 0.0_pReal) then ! compatible with equally signed dislocation density + rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) * transmissivity * compatibility ! subtract negative dislocation flux that enters the material point + elseif (compatibility < 0.0_pReal) then ! compatible with inversely signed dislocation density + topp = t + mod(t,2) - mod(t+1,2) + rhoDotFlux(s,topp) = rhoDotFlux(s,topp) - lineLength / mesh_ipVolume(ip,el) * transmissivity * abs(compatibility) ! subtract negative dislocation flux that enters the material point + endif + enddo endif enddo enddo @@ -1664,40 +1708,137 @@ endsubroutine !********************************************************************* -!* transmissivity of IP interface * +!* TRANSMISSIVITY AND COMPATIBILITY UPDATE * +!* Transmissivity is defined as absolute minimum of the cosine of * +!* the angle between the slip directions and the cosine of the angle * +!* between the slip plane normals. Compatibility is defined as * +!* product of cosine of the angle between the slip plane normals and * +!* signed cosine of the angle between the slip directions * !********************************************************************* -function constitutive_nonlocal_transmissivity(disorientation) +subroutine constitutive_nonlocal_updateCompatibility(orientation,i,e) -use prec, only: pReal, & - pInt -use math, only: math_QuaternionToAxisAngle +use prec, only: pReal, & + pInt +use math, only: math_QuaternionDisorientation, & + math_mul3x3, & + math_qRot +use material, only: material_phase, & + phase_constitution, & + phase_localConstitution, & + phase_constitutionInstance, & + homogenization_maxNgrains +use mesh, only: mesh_element, & + mesh_ipNeighborhood, & + FE_NipNeighbors, & + mesh_maxNips, & + mesh_NcpElems +use lattice, only: lattice_sn, & + lattice_sd, & + lattice_st, & + lattice_maxNslip +use debug, only: debugger, & + debug_e, debug_i, debug_g, & + verboseDebugger, & + selectiveDebugger implicit none !* input variables -real(pReal), dimension(4), intent(in) :: disorientation ! disorientation as quaternion +integer(pInt), intent(in) :: i, & ! ip index + e ! element index +real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + orientation ! crystal orientation in quaternions !* output variables -real(pReal) constitutive_nonlocal_transmissivity ! transmissivity of an IP interface for dislocations !* local variables -real(pReal) disorientationAngle -real(pReal), dimension(3) :: disorientationAxis -real(pReal), dimension(4) :: disorientationAxisAngle +integer(pInt) n, & ! neighbor index + neighboring_e, & ! element index of my neighbor + neighboring_i, & ! integration point index of my neighbor + myPhase, & ! phase + neighboringPhase, & + myInstance, & ! instance of constitution + neighboringInstance, & + myStructure, & ! lattice structure + neighboringStructure, & + myNSlipSystems, & ! number of active slip systems + neighboringNSlipSystems, & + c, & ! dislocation character index (1=edge, 2=screw) + s1, & ! slip system index (me) + s2 ! slip system index (my neighbor) +integer(pInt), dimension(lattice_maxNslip) :: mySlipSystems, & ! slip system numbering according to lattice + neighboringSlipSystems +real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor +real(pReal), dimension(3,lattice_maxNslip) :: myNormals, & ! slip plane normals + neighboringNormals +real(pReal), dimension(3,lattice_maxNslip,2) :: mySlipDirections, & ! slip directions for positive edge and screws + neighboringSlipDirections -disorientationAxisAngle = math_QuaternionToAxisAngle(disorientation) +myPhase = material_phase(1,i,e) +myInstance = phase_constitutionInstance(myPhase) +myStructure = constitutive_nonlocal_structure(myInstance) +myNSlipSystems = constitutive_nonlocal_totalNslip(myInstance) +mySlipSystems(1:myNSlipSystems) = constitutive_nonlocal_slipSystemLattice(1:myNSlipSystems,myInstance) +myNormals = lattice_sn(:, mySlipSystems, myStructure) +mySlipDirections(:,:,1) = lattice_sd(:, mySlipSystems, myStructure) ! direction of positive edges +mySlipDirections(:,:,2) = lattice_st(:, mySlipSystems, myStructure) ! direction of positive screws -disorientationAxis = disorientationAxisAngle(1:3) -disorientationAngle = disorientationAxisAngle(4) - -if (disorientationAngle < 5.0_pReal) then - constitutive_nonlocal_transmissivity = 1.0_pReal -else - constitutive_nonlocal_transmissivity = 0.5_pReal -endif +do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors + neighboring_e = mesh_ipNeighborhood(1,n,i,e) + neighboring_i = mesh_ipNeighborhood(2,n,i,e) -endfunction + if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists + neighboringPhase = material_phase(1,neighboring_i,neighboring_e) + + if (.not. phase_localConstitution(neighboringPhase)) then ! neighbor got also nonlocal constitution + neighboringInstance = phase_constitutionInstance(neighboringPhase) + neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) + neighboringNSlipSystems = constitutive_nonlocal_totalNslip(neighboringInstance) + neighboringSlipSystems(1:neighboringNSlipSystems) = constitutive_nonlocal_slipSystemLattice(1:neighboringNSlipSystems,& + neighboringInstance) + neighboringNormals = lattice_sn(:, neighboringSlipSystems, neighboringStructure) + neighboringSlipDirections(:,:,1) = lattice_sd(:, neighboringSlipSystems, neighboringStructure) ! direction of positive edges + neighboringSlipDirections(:,:,2) = lattice_st(:, neighboringSlipSystems, neighboringStructure) ! direction of positive screws + + if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me + absoluteMisorientation = math_QuaternionDisorientation( orientation(:,1,i,e), & + orientation(:,1,neighboring_i,neighboring_e), 0_pInt) + + do s1 = 1,myNSlipSystems ! loop through my slip systems + do c = 1,2 ! loop through edge and screw character + do s2 = 1,neighboringNSlipSystems ! loop through my neighbors' slip systems + constitutive_nonlocal_transmissivity(s2,c,s1,n,i,e) = & + min(abs(math_mul3x3(mySlipDirections(:,s1,c), & + math_qRot(absoluteMisorientation, neighboringSlipDirections(:,s2,c)))), & + abs(math_mul3x3(myNormals(:,s1), math_qRot(absoluteMisorientation, neighboringNormals(:,s2)))) ) + constitutive_nonlocal_compatibility(s2,c,s1,n,i,e) = & + abs(math_mul3x3(myNormals(:,s1), math_qRot(absoluteMisorientation, neighboringNormals(:,s2)))) & + * math_mul3x3(mySlipDirections(:,s1,c), math_qRot(absoluteMisorientation, neighboringSlipDirections(:,s2,c))) + enddo + if (any(abs(constitutive_nonlocal_compatibility(:,c,s1,n,i,e)) > 0.0_pReal)) then + constitutive_nonlocal_compatibility(:,c,s1,n,i,e) = constitutive_nonlocal_compatibility(:,c,s1,n,i,e) & + / sum(abs(constitutive_nonlocal_compatibility(:,c,s1,n,i,e))) ! normalize to a total of one + endif + enddo + + enddo + else ! neighbor has different crystal structure + constitutive_nonlocal_transmissivity(:,:,:,n,i,e) = 0.0_pReal ! no transmissivity... + constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 0.0_pReal ! ...and compatibility + endif + else ! neighbor has local constitution + constitutive_nonlocal_transmissivity(:,:,:,n,i,e) = 1.0_pReal ! assume perfectly transmissive... + constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 1.0_pReal ! ...and compatible + endif + else ! no neighbor present + constitutive_nonlocal_transmissivity(:,:,:,n,i,e) = 1.0_pReal ! perfectly transmissive... + constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 1.0_pReal ! ...and compatible + endif + +enddo + +endsubroutine diff --git a/code/crystallite.f90 b/code/crystallite.f90 index dcf308b82..abb7a81a0 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -184,7 +184,7 @@ subroutine crystallite_init(Temperature) allocate(crystallite_statedamper(gMax,iMax,eMax)); crystallite_statedamper = 1.0_pReal allocate(crystallite_symmetryID(gMax,iMax,eMax)); crystallite_symmetryID = 0.0_pReal !NEW allocate(crystallite_orientation(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal - allocate(crystallite_orientation0(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal + allocate(crystallite_orientation0(4,gMax,iMax,eMax)); crystallite_orientation0 = 0.0_pReal allocate(crystallite_rotation(4,gMax,iMax,eMax)); crystallite_rotation = 0.0_pReal allocate(crystallite_disorientation(4,nMax,gMax,iMax,eMax)); crystallite_disorientation = 0.0_pReal allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true. @@ -342,7 +342,7 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0) write(6,'(a35,x,7(i5,x))') 'crystallite_subF: ', shape(crystallite_subF) write(6,'(a35,x,7(i5,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0) - write(6,'(a35,x,7(i5,x))') 'crystallite_symmetryID: ', shape(crystallite_symmetryID) !NEW + write(6,'(a35,x,7(i5,x))') 'crystallite_symmetryID: ', shape(crystallite_symmetryID) write(6,'(a35,x,7(i5,x))') 'crystallite_subF0: ', shape(crystallite_subF0) write(6,'(a35,x,7(i5,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) write(6,'(a35,x,7(i5,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0) @@ -483,7 +483,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) logical forceLocalStiffnessCalculation ! flag indicating that stiffness calculation is always done locally forceLocalStiffnessCalculation = .true. - ! --+>> INITIALIZE TO STARTING CONDITION <<+-- crystallite_subStep = 0.0_pReal @@ -564,7 +563,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress - ! canÕt restore dotState here, since not yet calculated in first cutback after initialization + ! can�t restore dotState here, since not yet calculated in first cutback after initialization if (debugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,'(a78,f10.8)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& @@ -907,7 +906,7 @@ endif do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -919,10 +918,9 @@ RK4dotTemperature = 0.0_pReal !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g,i,e) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState @@ -968,8 +966,8 @@ do n = 1,4 !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -980,7 +978,7 @@ do n = 1,4 !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) if (crystallite_integrateStress(mode,g,i,e,timeStepFraction(n))) then ! fraction of original times step if (n == 4) then ! final integration step if (mode==1 .and. verboseDebugger .and. e == debug_e .and. i == debug_i .and. g == debug_g) then @@ -1020,10 +1018,9 @@ do n = 1,4 !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), & timeStepFraction(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) @@ -1201,8 +1198,8 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -1213,10 +1210,9 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g,i,e) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState @@ -1268,8 +1264,8 @@ do n = 1,5 !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -1279,7 +1275,7 @@ do n = 1,5 !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) if (crystallite_todo(g,i,e)) then if (.not. crystallite_integrateStress(mode,g,i,e,c(n))) then ! fraction of original time step if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... @@ -1299,11 +1295,11 @@ do n = 1,5 !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), c(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep + c(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState @@ -1402,8 +1398,8 @@ relTemperatureResiduum = 0.0_pReal if (crystallite_todo(g,i,e)) then sizeDotState = constitutive_sizeDotState(g,i,e) if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif endif enddo; enddo; enddo @@ -1415,7 +1411,7 @@ relTemperatureResiduum = 0.0_pReal !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) sizeDotState = constitutive_sizeDotState(g,i,e) if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then @@ -1544,8 +1540,8 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -1556,11 +1552,10 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g,i,e) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) @@ -1602,8 +1597,8 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -1613,7 +1608,7 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) if (crystallite_todo(g,i,e)) then if (.not. crystallite_integrateStress(mode,g,i,e)) then if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... @@ -1633,11 +1628,10 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g,i,e) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState @@ -1800,8 +1794,8 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -1812,10 +1806,9 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g,i,e) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState @@ -1865,8 +1858,8 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -1876,7 +1869,7 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) if (crystallite_todo(g,i,e)) then if (crystallite_integrateStress(mode,g,i,e)) then crystallite_converged(g,i,e) = .true. @@ -1984,8 +1977,8 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -1995,13 +1988,12 @@ endif !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) if (crystallite_todo(g,i,e)) then constitutive_previousDotState2(g,i,e)%p = 0.0_pReal constitutive_previousDotState(g,i,e)%p = 0.0_pReal - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g, i, e) endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -2039,7 +2031,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) if (crystallite_todo(g,i,e)) then crystallite_todo(g,i,e) = crystallite_integrateStress(mode,g,i,e) if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if broken non-local... @@ -2062,13 +2054,12 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) + selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) if (crystallite_todo(g,i,e)) then constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! wind forward dotStates constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & - crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & + crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g, i, e) endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -2117,8 +2108,8 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) !$OMP PARALLEL DO do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & - crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -2656,7 +2647,9 @@ use FEsolving, only: FEsolving_execElem, & use IO, only: IO_warning use material, only: material_phase, & homogenization_Ngrains, & - phase_constitution + phase_constitution, & + phase_localConstitution, & + phase_constitutionInstance use mesh, only: mesh_element, & mesh_ipNeighborhood, & FE_NipNeighbors @@ -2664,7 +2657,8 @@ use debug, only: debugger, & debug_e, debug_i, debug_g, & verboseDebugger, & selectiveDebugger -use constitutive_nonlocal, only: constitutive_nonlocal_label +use constitutive_nonlocal, only: constitutive_nonlocal_structure, & + constitutive_nonlocal_updateCompatibility implicit none @@ -2677,75 +2671,87 @@ integer(pInt) e, & ! element index i, & ! integration point index g, & ! grain index n, & ! neighbor index - myPhase, & ! phase neighboring_e, & ! element index of my neighbor neighboring_i, & ! integration point index of my neighbor - neighboringPhase, & ! phase of my neighbor - neighboringStructure ! lattice structure of my neighbor + myPhase, & ! phase + neighboringPhase, & + myInstance, & ! instance of constitution + neighboringInstance, & + myStructure, & ! lattice structure + neighboringStructure real(pReal), dimension(3,3) :: U, R logical error +! --- CALCULATE ORIENTATION AND LATTICE ROTATION --- !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) - - ! 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 call IO_warning(650, e, i, g) - crystallite_orientation(:,g,i,e) = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! fake orientation + crystallite_orientation(:,g,i,e) = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! fake orientation else crystallite_orientation(:,g,i,e) = math_RtoQuaternion(transpose(R)) endif - + crystallite_rotation(:,g,i,e) = & - math_QuaternionDisorientation( math_qConj(crystallite_orientation(:,g,i,e)), & ! calculate grainrotation + math_QuaternionDisorientation( math_qConj(crystallite_orientation(:,g,i,e)), & ! calculate grainrotation math_qConj(crystallite_orientation0(:,g,i,e)), & - 0_pInt ) ! we don't want symmetry here + 0_pInt ) ! we don't want symmetry here enddo enddo enddo !$OMPEND PARALLEL DO + +! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL --- +! --- we use crystallite_orientation from above, so need a seperate loop + !$OMP PARALLEL DO - ! Another loop for nonlocal material which uses the orientations from the first one. do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) selectiveDebugger = (e == debug_e .and. i == debug_i) - myPhase = material_phase(1,i,e) ! get my crystal structure - if (phase_constitution(myPhase) == constitutive_nonlocal_label) then ! if nonlocal model + myPhase = material_phase(1,i,e) ! get my phase + if (.not. phase_localConstitution(myPhase)) then ! if nonlocal model + myInstance = phase_constitutionInstance(myPhase) + myStructure = constitutive_nonlocal_structure(myInstance) ! get my crystal structure - do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors - + + ! --- calculate disorientation between me and my neighbor --- + + do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors neighboring_e = mesh_ipNeighborhood(1,n,i,e) neighboring_i = mesh_ipNeighborhood(2,n,i,e) - - 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 - if (myPhase == neighboringPhase) then ! if my neighbor has same phase like me - - crystallite_disorientation(:,n,1,i,e) = & - math_QuaternionDisorientation( crystallite_orientation(:,1,i,e), & - crystallite_orientation(:,1,neighboring_i,neighboring_e), & - crystallite_symmetryID(1,i,e)) ! calculate disorientation - - else ! for neighbor with different phase - crystallite_disorientation(:,n,1,i,e) = (/0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal/) ! 180 degree rotation about 100 axis - + if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists + neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase + if (.not. phase_localConstitution(neighboringPhase)) then ! neighbor got also nonlocal constitution + neighboringInstance = phase_constitutionInstance(neighboringPhase) + neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure + if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me + crystallite_disorientation(:,n,1,i,e) = & + math_QuaternionDisorientation( crystallite_orientation(:,1,i,e), & + crystallite_orientation(:,1,neighboring_i,neighboring_e), & + crystallite_symmetryID(1,i,e)) ! calculate disorientation + else ! for neighbor with different phase + crystallite_disorientation(:,n,1,i,e) = (/0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal/) ! 180 degree rotation about 100 axis + endif + else ! for neighbor with local constitution + crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity endif - else ! no existing neighbor - crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity - endif - if (verboseDebugger .and. selectiveDebugger) then - !$OMP CRITICAL (write2out) - write(6,'(a27,i2,a3,4(f12.5,x))') 'disorientation to neighbor ',n,' : ',crystallite_disorientation(:,n,1,i,e) - !$OMP END CRITICAL (write2out) + else ! no existing neighbor + crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity endif enddo + + + ! --- calculate compatibility and transmissivity between me and my neighbor --- + + call constitutive_nonlocal_updateCompatibility(crystallite_orientation,i,e) + endif enddo enddo