diff --git a/dislocation_based_subroutine/CPFEM.f90 b/dislocation_based_subroutine/CPFEM.f90 deleted file mode 100644 index 591801fe6..000000000 --- a/dislocation_based_subroutine/CPFEM.f90 +++ /dev/null @@ -1,488 +0,0 @@ - -!############################################################## - MODULE CPFEM -!############################################################## -! *** CPFEM engine *** -! - use prec, only: pReal,pInt - implicit none -! -! **************************************************************** -! *** General variables for the material behaviour calculation *** -! **************************************************************** - real(pReal) CPFEM_Tp - real(pReal), dimension (:,:,:), allocatable :: CPFEM_stress_all - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jacobi_all - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_all - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_all - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_results - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ini_ori - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_sigma_old - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_sigma_new - real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old - real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_old - integer(pInt) :: CPFEM_inc_old = 0_pInt - integer(pInt) :: CPFEM_subinc_old = 1_pInt - integer(pInt) :: CPFEM_Nresults = 3_pInt - logical :: CPFEM_first_call = .true. - - CONTAINS - -!********************************************************* -!*** allocate the arrays defined in module CPFEM *** -!*** and initialize them *** -!********************************************************* - SUBROUTINE CPFEM_init() -! - use prec, only: pReal,pInt - use math, only: math_EulertoR - use mesh - use constitutive -! - implicit none - - integer(pInt) e,i,g -! -! *** mpie.marc parameters *** - CPFEM_Tp = 0.0_pReal - allocate(CPFEM_ffn_all (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn_all = 0.0_pReal - allocate(CPFEM_ffn1_all (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_all = 0.0_pReal - allocate(CPFEM_stress_all( 6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_all = 0.0_pReal - allocate(CPFEM_jacobi_all(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jacobi_all = 0.0_pReal -! -! *** User defined results !!! MISSING incorporate consti_Nresults *** - allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) - CPFEM_results = 0.0_pReal -! -! *** Second Piola-Kirchoff stress tensor at (t=t0) and (t=t1) *** - allocate(CPFEM_sigma_old(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_sigma_old = 0.0_pReal - allocate(CPFEM_sigma_new(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_sigma_new = 0.0_pReal -! -! *** Plastic deformation gradient at (t=t0) and (t=t1) *** - allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) - forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:constitutive_maxNgrains) & - CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(constitutive_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation - allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal -! -! *** Old jacobian (consistent tangent) *** - allocate(CPFEM_jaco_old(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_old = 0.0_pReal -! -! *** Output to MARC output file *** - write(6,*) - write(6,*) 'Arrays allocated:' - write(6,*) 'CPFEM_ffn_all: ', shape(CPFEM_ffn_all) - write(6,*) 'CPFEM_ffn1_all: ', shape(CPFEM_ffn1_all) - write(6,*) 'CPFEM_stress_all: ', shape(CPFEM_stress_all) - write(6,*) 'CPFEM_jacobi_all: ', shape(CPFEM_jacobi_all) - write(6,*) 'CPFEM_results: ', shape(CPFEM_results) - write(6,*) 'CPFEM_sigma_old: ', shape(CPFEM_sigma_old) - write(6,*) 'CPFEM_sigma_new: ', shape(CPFEM_sigma_new) - write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old) - write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new) - write(6,*) 'CPFEM_jaco_old: ', shape(CPFEM_jaco_old) - write(6,*) - call flush(6) - return - - END SUBROUTINE -! -! -!*********************************************************************** -!*** perform initialization at first call, update variables and *** -!*** call the actual material model *** -!*********************************************************************** - SUBROUTINE CPFEM_general(ffn, ffn1, Tp, CPFEM_inc, CPFEM_subinc, CPFEM_cn, CPFEM_dt, CPFEM_en, CPFEM_in) -! - use prec, only: pReal,pInt - use math, only: math_init - use mesh, only: mesh_init,mesh_FEasCP - use constitutive, only: constitutive_init,constitutive_state_old,constitutive_state_new - implicit none -! - real(pReal) ffn(3,3), ffn1(3,3), Tp, CPFEM_dt - integer(pInt) CPFEM_inc, CPFEM_subinc, CPFEM_cn, CPFEM_en, CPFEM_in, cp_en -! -! initialization step - if (CPFEM_first_call) then -! three dimensional stress state ? - call math_init() - call mesh_init() - call constitutive_init() - call crystal_init() - call CPFEM_init() - CPFEM_first_call = .false. - endif - if (CPFEM_inc==CPFEM_inc_old) then ! not a new increment -! case of a new subincrement:update starting with subinc 2 - if (CPFEM_subinc > CPFEM_subinc_old) then - CPFEM_sigma_old = CPFEM_sigma_new - CPFEM_Fp_old = CPFEM_Fp_new - constitutive_state_old = constitutive_state_new - CPFEM_subinc_old = CPFEM_subinc - endif - else ! new increment - CPFEM_sigma_old = CPFEM_sigma_new - CPFEM_Fp_old = CPFEM_Fp_new - constitutive_state_old = constitutive_state_new - CPFEM_inc_old = CPFEM_inc - CPFEM_subinc_old = 1_pInt - endif - - cp_en = mesh_FEasCP('elem',CPFEM_en) - CPFEM_Tp = Tp - CPFEM_ffn_all(:,:,CPFEM_in, cp_en) = ffn - CPFEM_ffn1_all(:,:,CPFEM_in, cp_en) = ffn1 - call CPFEM_stressIP(CPFEM_cn, CPFEM_dt, cp_en, CPFEM_in) - return - - END SUBROUTINE - - -!********************************************************** -!*** calculate the material behaviour at IP level *** -!********************************************************** - SUBROUTINE CPFEM_stressIP(& - CPFEM_cn,& ! Cycle number - CPFEM_dt,& ! Time increment (dt) - cp_en,& ! Element number - CPFEM_in) ! Integration point number - - use prec, only: pReal,pInt,ijaco,nCutback - use math, only: math_pDecomposition,math_RtoEuler - use IO, only: IO_error - use mesh, only: mesh_element - use constitutive -! - implicit none - - integer(pInt), parameter :: i_now = 1_pInt,i_then = 2_pInt - character(len=128) msg - integer(pInt) CPFEM_cn,cp_en,CPFEM_in,grain,i - logical updateJaco,error - real(pReal) CPFEM_dt,dt,t,volfrac - real(pReal), dimension(6) :: cs,Tstar_v - real(pReal), dimension(6,6) :: cd - real(pReal), dimension(3,3) :: Fe,U,R,deltaFg - real(pReal), dimension(3,3,2) :: Fg,Fp - real(pReal), dimension(constitutive_maxNstatevars,2) :: state - - updateJaco = (mod(CPFEM_cn,ijaco)==0) ! update consistent tangent every ijaco'th iteration - - CPFEM_stress_all(:,CPFEM_in,cp_en) = 0.0_pReal ! average Cauchy stress - if (updateJaco) CPFEM_jaco_old(:,:,CPFEM_in,cp_en) = 0.0_pReal ! average consistent tangent - -! -------------- grain loop ----------------- - do grain = 1,texture_Ngrains(mesh_element(4,cp_en)) -! ------------------------------------------- - - i = 0_pInt ! cutback counter - state(:,i_now) = constitutive_state_old(:,grain,CPFEM_in,cp_en) - Fg(:,:,i_now) = CPFEM_ffn_all(:,:,CPFEM_in,cp_en) - Fp(:,:,i_now) = CPFEM_Fp_old(:,:,grain,CPFEM_in,cp_en) - - deltaFg = CPFEM_ffn1_all(:,:,CPFEM_in,cp_en)-CPFEM_ffn_all(:,:,CPFEM_in,cp_en) - dt = CPFEM_dt - - Tstar_v = CPFEM_sigma_old(:,grain,CPFEM_in,cp_en) ! use last result as initial guess - Fg(:,:,i_then) = Fg(:,:,i_now) - state(:,i_then) = 0.0_pReal ! state_old as initial guess - t = 0.0_pReal - -! ------- crystallite integration ----------- - do -! ------------------------------------------- - if (t+dt < CPFEM_dt) then ! intermediate solution - t = t+dt ! next time inc - Fg(:,:,i_then) = Fg(:,:,i_then)+deltaFg ! corresponding Fg - else ! full step solution - t = CPFEM_dt ! final time - Fg(:,:,i_then) = CPFEM_ffn1_all(:,:,CPFEM_in,cp_en) ! final Fg - endif - - call CPFEM_stressCrystallite(msg,cs,cd,Tstar_v,Fp(:,:,i_then),Fe,state(:,i_then),& - dt,cp_en,CPFEM_in,grain,updateJaco .and. t==CPFEM_dt,& - Fg(:,:,i_now),Fg(:,:,i_then),Fp(:,:,i_now),state(:,i_now)) - if (msg == 'ok') then ! solution converged - if (t == CPFEM_dt) exit ! reached final "then" - else ! solution not found - i = i+1_pInt ! inc cutback counter -! write(6,*) 'ncut:', i - if (i > nCutback) then ! limit exceeded? - write(6,*) 'cutback limit --> '//msg - write(6,*) 'Grain: ',grain - write(6,*) 'Integration point: ',CPFEM_in - write(6,*) 'Element: ',mesh_element(1,cp_en) - call IO_error(600) - return ! byebye - else - t = t-dt ! rewind time - Fg(:,:,i_then) = Fg(:,:,i_then)-deltaFg ! rewind Fg - dt = 0.5_pReal*dt ! cut time-step in half - deltaFg = 0.5_pReal*deltaFg ! cut Fg-step in half - endif - endif - enddo ! crystallite integration (cutback loop) - -! ---- update crystallite matrices at t = t1 ---- - CPFEM_Fp_new(:,:,grain,CPFEM_in,cp_en) = Fp(:,:,i_then) - constitutive_state_new(:,grain,CPFEM_in,cp_en) = state(:,i_then) - CPFEM_sigma_new(:,grain,CPFEM_in,cp_en) = Tstar_v -! ---- update results plotted in MENTAT ---- - call math_pDecomposition(Fe,U,R,error) ! polar decomposition - if (error) then - write(6,*) 'polar decomposition' - write(6,*) 'Grain: ',grain - write(6,*) 'Integration point: ',CPFEM_in - write(6,*) 'Element: ',mesh_element(1,cp_en) - call IO_error(600) - return - endif - CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R)) ! orientation - CPFEM_results(4:3+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en) = & - constitutive_post_results(Tstar_v,state(:,i_then),CPFEM_dt,CPFEM_Tp,grain,CPFEM_in,cp_en) - -! ---- contribute to IP result ---- - volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en) - CPFEM_stress_all(:,CPFEM_in,cp_en) = CPFEM_stress_all(:,CPFEM_in,cp_en)+volfrac*cs ! average Cauchy stress - if (updateJaco) CPFEM_jaco_old(:,:,CPFEM_in,cp_en) = CPFEM_jaco_old(:,:,CPFEM_in,cp_en)+volfrac*cd ! average consistent tangent - - enddo ! grain loop - - return - END SUBROUTINE - - -!******************************************************************** -! Calculates the stress for a single component -! it is based on the paper by Kalidindi et al.: -! J. Mech. Phys, Solids Vol. 40, No. 3, pp. 537-569, 1992 -! it is modified to use anisotropic elasticity matrix -!******************************************************************** - subroutine CPFEM_stressCrystallite(& - msg,& ! return message - cs,& ! Cauchy stress vector - dcs_de,& ! consistent tangent - Tstar_v,& ! second Piola-Kirchoff stress tensor - Fp_new,& ! new plastic deformation gradient - Fe_new,& ! new "elastic" deformation gradient - state_new,& ! new state variable array -! - dt,& ! time increment - cp_en,& ! element number - CPFEM_in,& ! integration point number - grain,& ! grain number - updateJaco,& ! boolean to calculate Jacobi matrix - Fg_old,& ! old global deformation gradient - Fg_new,& ! new global deformation gradient - Fp_old,& ! old plastic deformation gradient - state_old) ! old state variable array - - use prec, only: pReal,pInt,pert_e - use constitutive, only: constitutive_Nstatevars - use math, only: math_Mandel6to33,mapMandel - implicit none - - character(len=*) msg - logical updateJaco - integer(pInt) cp_en,CPFEM_in,grain,i - real(pReal) dt - real(pReal), dimension(3,3) :: Fg_old,Fg_new,Fg_pert,Fp_old,Fp_new,Fp_pert,Fe_new,Fe_pert,E_pert - real(pReal), dimension(6) :: cs,Tstar_v,Tstar_v_pert - real(pReal), dimension(6,6) :: dcs_de - real(pReal), dimension(constitutive_Nstatevars(grain,CPFEM_in,cp_en)) :: state_old,state_new,state_pert - - call CPFEM_timeIntegration(msg,Fp_new,Fe_new,Tstar_v,state_new, & ! def gradients and PK2 at end of time step - dt,cp_en,CPFEM_in,grain,Fg_new,Fp_old,state_old) - if (msg /= 'ok') return - cs = CPFEM_CauchyStress(Tstar_v,Fe_new) ! Cauchy stress - - if (updateJaco) then ! consistent tangent using numerical perturbation of Fg - do i = 1,6 ! Fg component - E_pert = 0.0_pReal - E_pert(mapMandel(1,i),mapMandel(2,i)) = E_pert(mapMandel(1,i),mapMandel(2,i)) + pert_e/2.0_pReal - E_pert(mapMandel(2,i),mapMandel(1,i)) = E_pert(mapMandel(2,i),mapMandel(1,i)) + pert_e/2.0_pReal - - Fg_pert = Fg_new+matmul(E_pert,Fg_old) ! perturbated Fg - Tstar_v_pert = Tstar_v ! initial guess from end of time step - state_pert = state_new ! initial guess from end of time step - call CPFEM_timeIntegration(msg,Fp_pert,Fe_pert,Tstar_v_pert,state_pert, & - dt,cp_en,CPFEM_in,grain,Fg_pert,Fp_old,state_old) - if (msg /= 'ok') then - msg = 'consistent tangent --> '//msg - return - endif -! Remark: (perturbated) Cauchy stress is Mandel hence dcs_de(:,4:6) is too large by sqrt(2) - dcs_de(:,i) = (CPFEM_CauchyStress(Tstar_v_pert,Fe_pert)-cs)/pert_e - enddo - endif - - return - - END SUBROUTINE - - -!*********************************************************************** -!*** fully-implicit two-level time integration *** -!*********************************************************************** - SUBROUTINE CPFEM_timeIntegration(& - msg,& ! return message - Fp_new,& ! new plastic deformation gradient - Fe_new,& ! new "elastic" deformation gradient - Tstar_v,& ! 2nd PK stress (taken as initial guess if /= 0) - state_new,& ! current microstructure at end of time inc (taken as guess if /= 0) -! - dt,& ! time increment - cp_en,& ! element number - CPFEM_in,& ! integration point number - grain,& ! grain number - Fg_new,& ! new total def gradient - Fp_old,& ! former plastic def gradient - state_old) ! former microstructure - - use prec - use constitutive, only: constitutive_Nstatevars,& - constitutive_HomogenizedC,constitutive_dotState,constitutive_LpAndItsTangent,& - constitutive_Microstructure - use math - implicit none - - character(len=*) msg - integer(pInt) cp_en, CPFEM_in, grain - integer(pInt) iState,iStress,dummy, i,j,k,l,m - real(pReal) dt,det, p_hydro - real(pReal), dimension(6) :: Tstar_v,dTstar_v,Rstress, E, T_elastic, Rstress_old - real(pReal), dimension(6,6) :: C_66,Jacobi,invJacobi - real(pReal), dimension(3,3) :: Fg_new,Fp_old,Fp_new,Fe_new,invFp_old,invFp_new,Lp,A,B,AB - real(pReal), dimension(3,3,3,3) :: dLp, LTL - real(pReal), dimension(constitutive_Nstatevars(grain, CPFEM_in, cp_en)) :: state_old,state_new,dstate,Rstate,RstateS - logical failed - - msg = 'ok' ! error-free so far - - call math_invert3x3(Fp_old,invFp_old,det,failed) ! inversion of Fp - if (failed) then - msg = 'inversion Fp_old' - return - endif - - C_66 = constitutive_HomogenizedC(grain, CPFEM_in, cp_en) - A = matmul(Fg_new,invFp_old) ! actually Fe - A = matmul(transpose(A), A) - -! former state guessed, if none specified - if (all(state_new == 0.0_pReal)) state_new = state_old - RstateS = state_new - iState = 0_pInt - - Rstress = Tstar_v - Rstress_old=Rstress - -state: do ! outer iteration: state - iState = iState+1 - if (iState > nState) then - msg = 'limit state iteration' - return - endif - call constitutive_Microstructure(state_new,CPFEM_Tp,grain,CPFEM_in,cp_en) - iStress = 0_pInt -stress: do ! inner iteration: stress - iStress = iStress+1 - if (iStress > nStress) then ! too many loops required - msg = 'limit stress iteration' - return - endif - p_hydro=(Tstar_v(1)+Tstar_v(2)+Tstar_v(3))/3.0_pReal - forall(i=1:3) Tstar_v(i)=Tstar_v(i)-p_hydro - call constitutive_LpAndItsTangent(Lp,dLp,Tstar_v,state_new,CPFEM_Tp,grain,CPFEM_in,cp_en) - B = math_I3-dt*Lp -! B = B / math_det3x3(B)**(1.0_pReal/3.0_pReal) - AB = matmul(A,B) - T_elastic= 0.5_pReal*matmul(C_66,math_Mandel33to6(matmul(transpose(B),AB)-math_I3)) - p_hydro=(T_elastic(1)+T_elastic(2)+T_elastic(3))/3.0_pReal - forall(i=1:3) T_elastic(i)=T_elastic(i)-p_hydro - Rstress = Tstar_v - T_elastic -! step size control: if residuum does not improve redo iteration with reduced step size - if(maxval(abs(Rstress)) > maxval(abs(Rstress_old)) .and. & - maxval(abs(Rstress)) > 1.0e-6 .and. iStress > 1) then -! write(6,*) 'Hallo', iStress - Tstar_v=Tstar_v+0.5*dTstar_v - dTstar_v=0.5*dTstar_v - cycle - endif - if (iStress > 1 .and. (maxval(abs(Tstar_v)) < 1.0e-3_pReal .or. maxval(abs(Rstress/maxval(abs(Tstar_v)))) < tol_Stress)) exit stress - -! update stress guess using inverse of dRes/dTstar (Newton--Raphson) - LTL = 0.0_pReal - do i=1,3 - do j=1,3 - do k=1,3 - do l=1,3 - do m=1,3 - LTL(i,j,k,l) = LTL(i,j,k,l) + dLp(j,i,m,k)*AB(m,l) + AB(m,i)*dLp(m,j,k,l) - enddo - enddo - enddo - enddo - enddo - Jacobi = math_identity2nd(6) + 0.5_pReal*dt*matmul(C_66,math_Mandel3333to66(LTL)) - j = 0_pInt - call math_invert6x6(Jacobi,invJacobi,dummy,failed) - do while (failed .and. j <= nReg) - forall (i=1:6) Jacobi(i,i) = 1.05_pReal*maxval(Jacobi(i,:)) ! regularization - call math_invert6x6(Jacobi,invJacobi,dummy,failed) - j = j+1 - enddo - if (failed) then - msg = 'regularization Jacobi' - return - endif - dTstar_v = matmul(invJacobi,Rstress) ! correction to Tstar - Rstress_old=Rstress - Tstar_v = Tstar_v-dTstar_v -! write(999,*) Tstar_v, dTstar_v, Rstress - - enddo stress -! write(6,*) 'istress', istress - Tstar_v = 0.5_pReal*matmul(C_66,math_Mandel33to6(matmul(transpose(B),AB)-math_I3)) - dstate = dt*constitutive_dotState(Tstar_v,state_new,CPFEM_Tp,grain,CPFEM_in,cp_en) ! evolution of microstructure - Rstate = state_new - (state_old+dstate) - RstateS = 0.0_pReal - forall (i=1:constitutive_Nstatevars(grain,CPFEM_in,cp_en), state_new(i)/=0.0_pReal) & - RstateS(i) = Rstate(i)/state_new(i) - state_new = state_old+dstate - if (maxval(abs(RstateS)) < tol_State) exit state - - enddo state -! write(6,*) 'istate', istate -! write(999,*) 'Tstar_v raus', Tstar_v -! write(999,*) - - invFp_new = matmul(invFp_old,B) - call math_invert3x3(invFp_new,Fp_new,det,failed) - if (failed) then - msg = 'inversion Fp_new' - return - endif - Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! det = det(InvFp_new) !! - Fe_new = matmul(Fg_new,invFp_new) - return - END SUBROUTINE - - - FUNCTION CPFEM_CauchyStress(PK_v,Fe) -!*********************************************************************** -!*** Cauchy stress calculation *** -!*********************************************************************** - use prec, only: pReal,pInt - use math, only: math_Mandel33to6,math_Mandel6to33,math_det3x3 - implicit none -! *** Subroutine parameters *** - real(pReal) PK_v(6), Fe(3,3), CPFEM_CauchyStress(6) - - CPFEM_CauchyStress = math_Mandel33to6(matmul(matmul(Fe,math_Mandel6to33(PK_v)),transpose(Fe))/math_det3x3(Fe)) - return - END FUNCTION - - - END MODULE diff --git a/dislocation_based_subroutine/IO.f90 b/dislocation_based_subroutine/IO.f90 deleted file mode 100644 index bd20f0476..000000000 --- a/dislocation_based_subroutine/IO.f90 +++ /dev/null @@ -1,555 +0,0 @@ - -!############################################################## - MODULE IO -!############################################################## - - CONTAINS -!--------------------------- -! function IO_open_file(unit,relPath) -! function IO_open_inputFile(unit) -! function IO_hybridIA(Nast,ODFfileName) -! private function hybridIA_reps(dV_V,steps,C) -! function IO_stringPos(line,maxN) -! function IO_stringValue(line,positions,pos) -! function IO_floatValue(line,positions,pos) -! function IO_intValue(line,positions,pos) -! function IO_fixedStringValue(line,ends,pos) -! function IO_fixedFloatValue(line,ends,pos) -! function IO_fixedFloatNoEValue(line,ends,pos) -! function IO_fixedIntValue(line,ends,pos) -! function IO_continousTntValues(unit,maxN) -! function IO_lc(line) -! subroutine IO_lcInplace(line) -! subroutine IO_error(ID) -!--------------------------- - - - -!******************************************************************** -! open existing file to given unit -! path to file is relative to working directory -!******************************************************************** - logical FUNCTION IO_open_file(unit,relPath) - - use prec, only: pInt - implicit none - - character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ - character(len=*) relPath - integer(pInt) unit - character(256) path - - inquire(6, name=path) ! determine outputfile - open(unit,status='old',err=100,file=path(1:scan(path,pathSep,back=.true.))//relPath) - IO_open_file = .true. - return -100 IO_open_file = .false. - return - - END FUNCTION - - -!******************************************************************** -! open FEM inputfile to given unit -!******************************************************************** - logical FUNCTION IO_open_inputFile(unit) - - use prec, only: pReal, pInt - implicit none - - character(256) outName - integer(pInt) unit, extPos - character(3) ext - - inquire(6, name=outName) ! determine outputfileName - extPos = len_trim(outName)-2 - if(outName(extPos:extPos+2)=='out') then - ext='dat' ! MARC - else - ext='inp' ! ABAQUS - end if - open(unit,status='old',err=100,file=outName(1:extPos-1)//ext) - IO_open_inputFile = .true. - return -100 IO_open_inputFile = .false. - return - - END FUNCTION - - -!******************************************************************** -! hybrid IA repetition counter -!******************************************************************** - FUNCTION hybridIA_reps(dV_V,steps,C) - - use prec, only: pReal, pInt - implicit none - - integer(pInt), intent(in), dimension(3) :: steps - integer(pInt) hybridIA_reps, phi1,Phi,phi2 - real(pReal), intent(in), dimension(steps(3),steps(2),steps(1)) :: dV_V - real(pReal), intent(in) :: C - - hybridIA_reps = 0_pInt - do phi1=1,steps(1) - do Phi =1,steps(2) - do phi2=1,steps(3) - hybridIA_reps = hybridIA_reps+nint(C*dV_V(phi2,Phi,phi1), pInt) - end do - end do - end do - return - - END FUNCTION - - -!******************************************************************** -! hybrid IA sampling of ODFfile -!******************************************************************** - FUNCTION IO_hybridIA(Nast,ODFfileName) - - use prec, only: pReal, pInt - use math, only: inRad - - implicit none - - character(len=*) ODFfileName - character(len=80) line - character(len=*), parameter :: fileFormat = '(A80)' - integer(pInt) i,j,bin,Nast,NnonZero,Nset,Nreps,reps,phi1,Phi,phi2 - integer(pInt), dimension(7) :: pos - integer(pInt), dimension(3) :: steps - integer(pInt), dimension(:), allocatable :: binSet - real(pReal) center,sum_dV_V,prob,dg_0,C,lowerC,upperC,rnd - real(pReal), dimension(3) :: limits,deltas - real(pReal), dimension(:,:,:), allocatable :: dV_V - real(pReal), dimension(3,Nast) :: IO_hybridIA - - if (.not. IO_open_file(999,ODFfileName)) goto 100 - -!--- parse header of ODF file --- -!--- limits in phi1, Phi, phi2 --- - read(999,fmt=fileFormat,end=100) line - pos = IO_stringPos(line,3) - if (pos(1).ne.3) goto 100 - do i=1,3 - limits(i) = IO_intValue(line,pos,i)*inRad - end do - -!--- deltas in phi1, Phi, phi2 --- - read(999,fmt=fileFormat,end=100) line - pos = IO_stringPos(line,3) - if (pos(1).ne.3) goto 100 - do i=1,3 - deltas(i) = IO_intValue(line,pos,i)*inRad - end do - steps = nint(limits/deltas,pInt) - allocate(dV_V(steps(3),steps(2),steps(1))) - -!--- box boundary/center at origin? --- - read(999,fmt=fileFormat,end=100) line - if (index(IO_lc(line),'bound')>0) then - center = 0.5_pReal - else - center = 0.0_pReal - end if - -!--- skip blank line --- - read(999,fmt=fileFormat,end=100) line - - sum_dV_V = 0.0_pReal - dV_V = 0.0_pReal - dg_0 = deltas(1)*deltas(3)*2.0_pReal*sin(deltas(2)/2.0_pReal) - NnonZero = 0_pInt - - do phi1=1,steps(1) - do Phi=1,steps(2) - do phi2=1,steps(3) - read(999,fmt='(F)',end=100) prob - if (prob > 0.0_pReal) then - NnonZero = NnonZero+1 - sum_dV_V = sum_dV_V+prob - else - prob = 0.0_pReal - end if - dV_V(phi2,Phi,phi1) = prob*dg_0*sin((Phi-1.0_pReal+center)*deltas(2)) - end do - end do - end do - - dV_V = dV_V/sum_dV_V ! normalize to 1 - -!--- now fix bounds --- - Nset = max(Nast,NnonZero) - lowerC = 0.0_pReal - upperC = real(Nset, pReal) - - do while (hybridIA_reps(dV_V,steps,upperC) < Nset) - lowerC = upperC - upperC = upperC*2.0_pReal - end do -!--- binary search for best C --- - do - C = (upperC+lowerC)/2.0_pReal - Nreps = hybridIA_reps(dV_V,steps,C) - if (abs(upperC-lowerC) < upperC*1.0e-14_pReal) then - C = upperC - Nreps = hybridIA_reps(dV_V,steps,C) - exit - elseif (Nreps < Nset) then - lowerC = C - elseif (Nreps > Nset) then - upperC = C - else - exit - end if - end do - - allocate(binSet(Nreps)) - bin = 0 ! bin counter - i = 1 ! set counter - do phi1=1,steps(1) - do Phi=1,steps(2) - do phi2=1,steps(3) - reps = nint(C*dV_V(phi2,Phi,phi1), pInt) - binSet(i:i+reps-1) = bin - bin = bin+1 ! advance bin - i = i+reps ! advance set - end do - end do - end do - - do i=1,Nast - if (i < Nast) then - call random_number(rnd) - j = nint(rnd*(Nast-i)+i+0.5_pReal,pInt) - else - j = i - end if - bin = binSet(j) - IO_hybridIA(1,i) = deltas(1)*(mod(bin/(steps(3)*steps(2)),steps(1))+center) ! phi1 - IO_hybridIA(2,i) = deltas(2)*(mod(bin/ steps(3) ,steps(2))+center) ! Phi - IO_hybridIA(3,i) = deltas(3)*(mod(bin ,steps(3))+center) ! phi2 - binSet(j) = binSet(i) - end do - close(999) - return - -! on error -100 IO_hybridIA = -1 - close(999) - return - - END FUNCTION - - -!******************************************************************** -! locate at most N space-separated parts in line -! return array containing number of parts found and -! their left/right positions to be used by IO_xxxVal -!******************************************************************** - FUNCTION IO_stringPos (line,N) - - use prec, only: pReal,pInt - implicit none - - character(len=*) line - character(len=*), parameter :: sep=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces - integer(pInt) N, part - integer(pInt) IO_stringPos(1+N*2) - - IO_stringPos = -1 - IO_stringPos(1) = 0 - part = 1 - do while ((N<1 .or. part<=N) .and. verify(line(IO_stringPos(part*2-1)+1:),sep)>0) - IO_stringPos(part*2) = IO_stringPos(part*2-1)+verify(line(IO_stringPos(part*2-1)+1:),sep) - IO_stringPos(part*2+1) = IO_stringPos(part*2)+scan(line(IO_stringPos(part*2):),sep)-2 - part = part+1 - end do - IO_stringPos(1) = part-1 - return - - END FUNCTION - - -!******************************************************************** -! read string value at pos from line -!******************************************************************** - FUNCTION IO_stringValue (line,positions,pos) - - use prec, only: pReal,pInt - implicit none - - character(len=*) line - integer(pInt) positions(*),pos - character(len=1+positions(pos*2+1)-positions(pos*2)) IO_stringValue - - if (positions(1) < pos) then - IO_stringValue = '' - else - IO_stringValue = line(positions(pos*2):positions(pos*2+1)) - endif - return - - END FUNCTION - - -!******************************************************************** -! read string value at pos from fixed format line -!******************************************************************** - FUNCTION IO_fixedStringValue (line,ends,pos) - - use prec, only: pReal,pInt - implicit none - - character(len=*) line - integer(pInt) ends(*),pos - character(len=ends(pos+1)-ends(pos)) IO_fixedStringValue - - IO_fixedStringValue = line(ends(pos)+1:ends(pos+1)) - return - - END FUNCTION - - -!******************************************************************** -! read float value at pos from line -!******************************************************************** - FUNCTION IO_floatValue (line,positions,pos) - - use prec, only: pReal,pInt - implicit none - - character(len=*) line - real(pReal) IO_floatValue - integer(pInt) positions(*),pos - - if (positions(1) >= pos) then - read(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT=*) IO_floatValue - return - endif -100 IO_floatValue = huge(1.0_pReal) - return - - END FUNCTION - - -!******************************************************************** -! read float value at pos from fixed format line -!******************************************************************** - FUNCTION IO_fixedFloatValue (line,ends,pos) - - use prec, only: pReal,pInt - implicit none - - character(len=*) line - real(pReal) IO_fixedFloatValue - integer(pInt) ends(*),pos - - read(UNIT=line(ends(pos-1)+1:ends(pos)),ERR=100,FMT=*) IO_fixedFloatValue - return -100 IO_fixedFloatValue = huge(1.0_pReal) - return - - END FUNCTION - - -!******************************************************************** -! read float x.y+z value at pos from format line line -!******************************************************************** - FUNCTION IO_fixedNoEFloatValue (line,ends,pos) - - use prec, only: pReal,pInt - implicit none - - character(len=*) line - real(pReal) IO_fixedNoEFloatValue,base - integer(pInt) ends(*),pos,pos_exp,expon - - pos_exp = scan(line(ends(pos)+1:ends(pos+1)),'+-',back=.true.) - if (pos_exp > 1) then - read(UNIT=line(ends(pos)+1:ends(pos)+pos_exp-1),ERR=100,FMT='(F)') base - read(UNIT=line(ends(pos)+pos_exp:ends(pos+1)),ERR=100,FMT='(I)') expon - else - read(UNIT=line(ends(pos)+1:ends(pos+1)),ERR=100,FMT=*) base - expon = 0_pInt - endif - IO_fixedNoEFloatValue = base*10.0_pReal**expon - return -100 IO_fixedNoEFloatValue = huge(1.0_pReal) - return - - END FUNCTION - - -!******************************************************************** -! read int value at pos from line -!******************************************************************** - FUNCTION IO_intValue (line,positions,pos) - - use prec, only: pReal,pInt - implicit none - - character(len=*) line - integer(pInt) IO_intValue - integer(pInt) positions(*),pos - - if (positions(1) >= pos) then - read(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT='(I)') IO_intValue - return - endif -100 IO_intValue = huge(1_pInt) - return - - END FUNCTION - - -!******************************************************************** -! read int value at pos from fixed format line -!******************************************************************** - FUNCTION IO_fixedIntValue (line,ends,pos) - - use prec, only: pReal,pInt - implicit none - - character(len=*) line - integer(pInt) IO_fixedIntValue - integer(pInt) ends(*),pos - - read(UNIT=line(ends(pos)+1:ends(pos+1)),ERR=100,FMT='(I)') IO_fixedIntValue - return -100 IO_fixedIntValue = huge(1_pInt) - return - - END FUNCTION - - -!******************************************************************** -! change character in line to lower case -!******************************************************************** - FUNCTION IO_lc (line) - - use prec, only: pInt - implicit none - - character (len=*) line - character (len=len(line)) IO_lc - integer(pInt) i - - IO_lc = line - do i=1,len(line) - if(64') then - part=tag(2:len_trim(tag)-1) - exit - elseif (tag(1:1)=='[') then - count=count+1 - endif -enddo -100 return - -end subroutine - - -character(len=80) function constitutive_assignNGaussAndFiber(file) -!********************************************************************* -!********************************************************************* -use prec, only: pInt -use IO, only: IO_stringPos,IO_stringValue,IO_lc -implicit none - -!* Definition of variables -character(len=80) line,tag -integer(pInt) file,section -integer(pInt), dimension(3) :: positions - -constitutive_assignNGaussAndFiber='' -section = 0_pInt - -do - read(file,'(a80)',END=100) line - positions=IO_stringPos(line,1) - tag=IO_lc(IO_stringValue(line,positions,1)) - if (tag(1:1)=='#' .OR. positions(1)==0) then ! skip comment and empty lines - cycle - elseif (tag(1:1)=='<'.AND.tag(len_trim(tag):len_trim(tag))=='>') then - constitutive_assignNGaussAndFiber=tag(2:len_trim(tag)-1) - exit - elseif (tag(1:1)=='[') then - section=section+1 - texture_NGauss(section) = 0_pInt - texture_NFiber(section) = 0_pInt - elseif (tag=='(gauss)') then - texture_NGauss(section)=texture_NGauss(section)+1 - elseif (tag=='(fiber)') then - texture_NFiber(section)=texture_NFiber(section)+1 - endif -enddo -100 return - -end function - - -character(len=80) function constitutive_Parse_UnknownPart(file) -!********************************************************************* -!* read an unknown "part" from the input file until * -!* the next part is reached * -!* INPUT: * -!* - file : file ID * -!********************************************************************* -use prec, only: pInt -use IO, only: IO_stringPos,IO_stringValue,IO_lc -implicit none - -!* Definition of variables -character(len=80) line,tag -integer(pInt), parameter :: maxNchunks = 1 -integer(pInt) file -integer(pInt), dimension(1+2*maxNchunks) :: positions - -constitutive_parse_unknownPart='' - -do - read(file,'(a80)',END=100) line - positions=IO_stringPos(line,maxNchunks) - tag=IO_lc(IO_stringValue(line,positions,1)) - if (tag(1:1)=='#' .OR. positions(1)==0) then ! skip comment and empty lines - cycle - elseif (tag(1:1)=='<'.AND.tag(len_trim(tag):len_trim(tag))=='>') then - constitutive_Parse_UnknownPart=tag(2:len_trim(tag)-1) - exit - endif -enddo -100 return - -end function - - -character(len=80) function constitutive_Parse_MaterialPart(file) -!********************************************************************* -!* This function reads a material "part" from the input file until * -!* the next part is reached * -!* INPUT: * -!* - file : file ID * -!********************************************************************* -use prec, only: pInt,pReal -use IO -use crystal, only: crystal_MaxMaxNslipOfStructure -implicit none - -!* Definition of variables -character(len=80) line,tag -integer(pInt), parameter :: maxNchunks = 2 -integer(pInt) i,file,section -integer(pInt), dimension(1+2*maxNchunks) :: positions - -section = 0 -constitutive_parse_materialPart = '' - -do while(.true.) - read(file,'(a80)',END=100) line - positions=IO_stringPos(line,maxNchunks) ! parse leading chunks - tag=IO_lc(IO_stringValue(line,positions,1)) - if (tag(1:1)=='#' .OR. positions(1)==0) then ! skip comment and empty lines - cycle - elseif (tag(1:1)=='<'.AND.tag(len_trim(tag):len_trim(tag))=='>') then - constitutive_parse_materialPart=tag(2:len_trim(tag)-1) - exit - elseif (tag(1:1)=='[') then - section=section+1 - else - if (section>0) then - select case(tag) - case ('crystal_structure') - material_CrystalStructure(section)=IO_intValue(line,positions,2) - case ('nslip') - material_Nslip(section)=IO_intValue(line,positions,2) - case ('c11') - material_C11(section)=IO_floatValue(line,positions,2) - case ('c12') - material_C12(section)=IO_floatValue(line,positions,2) - case ('c13') - material_C13(section)=IO_floatValue(line,positions,2) - case ('c33') - material_C33(section)=IO_floatValue(line,positions,2) - case ('c44') - material_C44(section)=IO_floatValue(line,positions,2) - case ('rho0') !* conversion in 1/mm˛ - material_rho0(section)=IO_floatValue(line,positions,2)/1.0e6_pReal - case ('interaction_coefficient') - do i=1,crystal_MaxMaxNslipOfStructure - material_SlipIntCoeff(i,section)=IO_floatValue(line,positions,i+1) - enddo - case ('bg') !* conversion in mm - material_bg(section)=IO_floatValue(line,positions,2)*1.0e3_pReal - case ('Qedge') !* conversion in mJ/Kelvin - material_Qedge(section)=IO_floatValue(line,positions,2)*1.0e3_pReal - case ('tau0') - material_tau0(section)=IO_floatValue(line,positions,2) - case ('c1') - material_c1(section)=IO_floatValue(line,positions,2) - case ('c2') - material_c2(section)=IO_floatValue(line,positions,2) - case ('c3') - material_c3(section)=IO_floatValue(line,positions,2) - case ('c4') - material_c4(section)=IO_floatValue(line,positions,2) - case ('c5') - material_c5(section)=IO_floatValue(line,positions,2) - end select - endif - endif -enddo -100 return - -end function - - -character(len=80) function constitutive_Parse_TexturePart(file) -!********************************************************************* -!* This function reads a texture "part" from the input file until * -!* the next part is reached * -!* INPUT: * -!* - file : file ID * -!********************************************************************* -use prec, only: pInt -use IO -use math, only: inRad -implicit none - -!* Definition of variables -character(len=80) line,tag -integer(pInt), parameter :: maxNchunks = 13 ! may be more than 10 chunks ..? -integer(pInt) file,section,gaussCount,fiberCount,i -integer(pInt), dimension(1+2*maxNchunks) :: positions - -section = 0 -gaussCount = 0 -fiberCount = 0 -constitutive_parse_texturePart = '' - -do while(.true.) - read(file,'(a80)',END=100) line - positions=IO_stringPos(line,maxNchunks) ! parse leading chunks - tag=IO_lc(IO_stringValue(line,positions,1)) - if (tag(1:1)=='#' .OR. positions(1)==0) then ! skip comment and empty lines - cycle - elseif (tag(1:1)=='<'.AND.tag(len_trim(tag):len_trim(tag))=='>') then - constitutive_parse_texturePart=tag(2:len_trim(tag)-1) - exit - elseif (tag(1:1)=='[') then - section=section+1 - gaussCount=0 - fiberCount=0 - else - if (section>0) then - select case(tag) - case ('hybridIA') - texture_ODFfile(section)=IO_stringValue(line,positions,2) - case ('(gauss)') - gaussCount=gaussCount+1 - do i=2,10,2 - tag=IO_lc(IO_stringValue(line,positions,i)) - select case (tag) - case('phi1') - texture_Gauss(1,gaussCount,section)=IO_floatValue(line,positions,i+1)*inRad - case('phi') - texture_Gauss(2,gaussCount,section)=IO_floatValue(line,positions,i+1)*inRad - case('phi2') - texture_Gauss(3,gaussCount,section)=IO_floatValue(line,positions,i+1)*inRad - case('scatter') - texture_Gauss(4,gaussCount,section)=IO_floatValue(line,positions,i+1)*inRad - case('fraction') - texture_Gauss(5,gaussCount,section)=IO_floatValue(line,positions,i+1) - end select - enddo - case ('(fiber)') - fiberCount=fiberCount+1 - do i=2,12,2 - tag=IO_lc(IO_stringValue(line,positions,i)) - select case (tag) - case('alpha1') - texture_fiber(1,fiberCount,section)=IO_floatValue(line,positions,i+1)*inRad - case('alpha2') - texture_fiber(2,fiberCount,section)=IO_floatValue(line,positions,i+1)*inRad - case('beta1') - texture_fiber(3,fiberCount,section)=IO_floatValue(line,positions,i+1)*inRad - case('beta2') - texture_fiber(4,fiberCount,section)=IO_floatValue(line,positions,i+1)*inRad - case('scatter') - texture_fiber(5,fiberCount,section)=IO_floatValue(line,positions,i+1)*inRad - case('fraction') - texture_fiber(6,fiberCount,section)=IO_floatValue(line,positions,i+1) - end select - enddo - case ('ngrains') - texture_Ngrains(section)=IO_intValue(line,positions,2) - case ('symmetry') - texture_symmetry(section)=IO_stringValue(line,positions,2) - end select - endif - endif -enddo -100 return - -end function - - -subroutine constitutive_Parse_MatTexDat(filename) -!********************************************************************* -!* This function reads the material and texture input file * -!* INPUT: * -!* - filename : name of input file * -!********************************************************************* -use prec, only: pReal,pInt -use IO, only: IO_error, IO_open_file -use math, only: math_Mandel3333to66, math_Voigt66to3333 -use crystal, only: crystal_MaxMaxNslipOfStructure -implicit none - -!* Definition of variables -character(len=*) filename -character(len=80) part,formerPart -integer(pInt) sectionCount,i,j,k, fileunit - -! set fileunit -fileunit=200 -!----------------------------- -!* First reading: number of materials and textures -!----------------------------- -!* determine material_maxN and texture_maxN from last respective parts -if(IO_open_file(fileunit,filename)==.false.) goto 100 -part = '_dummy_' -do while (part/='') - formerPart = part - call constitutive_CountSections(fileunit,sectionCount,part) - select case (formerPart) - case ('materials') - material_maxN = sectionCount - case ('textures') - texture_maxN = sectionCount - end select -enddo -!* Array allocation -allocate(material_CrystalStructure(material_maxN)) ; material_CrystalStructure=0_pInt -allocate(material_Nslip(material_maxN)) ; material_Nslip=0_pInt -allocate(material_C11(material_maxN)) ; material_C11=0.0_pReal -allocate(material_C12(material_maxN)) ; material_C12=0.0_pReal -allocate(material_C13(material_maxN)) ; material_C13=0.0_pReal -allocate(material_C33(material_maxN)) ; material_C33=0.0_pReal -allocate(material_C44(material_maxN)) ; material_C44=0.0_pReal -allocate(material_Gmod(material_maxN)) ; material_Gmod=0.0_pReal -allocate(material_Cslip_66(6,6,material_maxN)) ; material_Cslip_66=0.0_pReal -allocate(material_rho0(material_maxN)) ; material_rho0=0.0_pReal -allocate(material_SlipIntCoeff(crystal_MaxMaxNslipOfStructure,material_maxN)) ; material_SlipIntCoeff=0.0_pReal -allocate(material_bg(material_maxN)) ; material_bg=0.0_pReal -allocate(material_Qedge(material_maxN)) ; material_Qedge=0.0_pReal -allocate(material_tau0(material_maxN)) ; material_tau0=0.0_pReal -allocate(material_c1(material_maxN)) ; material_c1=0.0_pReal -allocate(material_c2(material_maxN)) ; material_c2=0.0_pReal -allocate(material_c3(material_maxN)) ; material_c3=0.0_pReal -allocate(material_c4(material_maxN)) ; material_c4=0.0_pReal -allocate(material_c5(material_maxN)) ; material_c5=0.0_pReal -allocate(texture_ODFfile(texture_maxN)) ; texture_ODFfile='' -allocate(texture_Ngrains(texture_maxN)) ; texture_Ngrains=0_pInt -allocate(texture_symmetry(texture_maxN)) ; texture_symmetry='' -allocate(texture_NGauss(texture_maxN)) ; texture_NGauss=0_pInt -allocate(texture_NFiber(texture_maxN)) ; texture_NFiber=0_pInt -allocate(texture_NRandom(texture_maxN)) ; texture_NRandom=0_pInt - -!----------------------------- -!* Second reading: number of Gauss and Fiber -!----------------------------- -rewind(fileunit) -part = '_dummy_' -do while (part/='') - select case (part) - case ('textures') - part = constitutive_assignNGaussAndFiber(fileunit) - case default - part = constitutive_Parse_UnknownPart(fileunit) - end select -enddo -!* Array allocation -texture_maxNGauss=maxval(texture_NGauss) -texture_maxNFiber=maxval(texture_NFiber) -allocate(texture_Gauss(5,texture_maxNGauss,texture_maxN)) ; texture_Gauss=0.0_pReal -allocate(texture_Fiber(6,texture_maxNFiber,texture_maxN)) ; texture_Fiber=0.0_pReal - -!----------------------------- -!* Third reading: materials and textures are stored -!----------------------------- -rewind(fileunit) -part='_dummy_' -do while (part/='') - select case (part) - case ('materials') - part=constitutive_Parse_MaterialPart(fileunit) - case ('textures') - part=constitutive_Parse_TexturePart(fileunit) - case default - part=constitutive_Parse_UnknownPart(fileunit) - end select -enddo -close(fileunit) - - -!* Construction of the elasticity matrices -do i=1,material_maxN - material_Gmod(i)=material_C44(i) - select case (material_CrystalStructure(i)) - case(1:2) ! cubic(s) - do k=1,3 - do j=1,3 - material_Cslip_66(k,j,i)=material_C12(i) - enddo - material_Cslip_66(k,k,i)=material_C11(i) - material_Cslip_66(k+3,k+3,i)=material_C44(i) - enddo - case(3) ! hcp - material_Cslip_66(1,1,i)=material_C11(i) - material_Cslip_66(2,2,i)=material_C11(i) - material_Cslip_66(3,3,i)=material_C33(i) - material_Cslip_66(1,2,i)=material_C12(i) - material_Cslip_66(2,1,i)=material_C12(i) - material_Cslip_66(1,3,i)=material_C13(i) - material_Cslip_66(3,1,i)=material_C13(i) - material_Cslip_66(2,3,i)=material_C13(i) - material_Cslip_66(3,2,i)=material_C13(i) - material_Cslip_66(4,4,i)=material_C44(i) - material_Cslip_66(5,5,i)=material_C44(i) - material_Cslip_66(6,6,i)=0.5_pReal*(material_C11(i)-material_C12(i)) - end select - material_Cslip_66(:,:,i) = math_Mandel3333to66(math_Voigt66to3333(material_Cslip_66(:,:,i))) -enddo - - -! MISSING some consistency checks may be..? -! if ODFfile present then set NGauss NFiber =0 -return -100 call IO_error(200) ! corrupt materials_textures file -end subroutine - - -subroutine constitutive_Assignment() -!********************************************************************* -!* This subroutine assign material parameters according to ipc,ip,el * -!********************************************************************* -use prec, only: pReal,pInt -use math, only: math_sampleGaussOri,math_sampleFiberOri,math_sampleRandomOri,math_symmetricEulers,math_EulerToR -use mesh, only: mesh_NcpElems,FE_Nips,FE_mapElemtype,mesh_maxNips,mesh_element -use IO, only: IO_hybridIA -use crystal, only: crystal_MaxNslipOfStructure,crystal_SlipIntType,crystal_sn,crystal_st -implicit none - -!* Definition of variables -integer(pInt) e,i,j,k,l,m,o,g,s -integer(pInt) matID,texID -real(pReal) K_inter,x,y -integer(pInt), dimension(:,:,:), allocatable :: hybridIA_population -integer(pInt), dimension(texture_maxN) :: Ncomponents,Nsym,multiplicity,sumVolfrac,ODFmap,sampleCount -real(pReal), dimension(3,4*(1+texture_maxNGauss+texture_maxNfiber)) :: Euler -real(pReal), dimension(4*(1+texture_maxNGauss+texture_maxNfiber)) :: texVolfrac - -! process textures -o = 0_pInt ! ODF counter -ODFmap = 0_pInt ! blank mapping -sampleCount = 0_pInt ! count orientations assigned per texture - -do texID=1,texture_maxN - if (texture_ODFfile(texID)=='') then - sumVolfrac(texID) = sum(texture_gauss(5,:,texID))+sum(texture_fiber(6,:,texID)) - if (sumVolfrac(texID)<1.0_pReal) texture_NRandom(texID) = 1_pInt ! check whether random component missing - select case (texture_symmetry(texID)) ! set symmetry factor - case ('orthotropic') - Nsym(texID) = 4_pInt - case ('monoclinic') - Nsym(texID) = 2_pInt - case default - Nsym(texID) = 1_pInt - end select - Ncomponents(texID) = texture_NGauss(texID)+texture_NFiber(texID)+texture_NRandom(texID) - else ! hybrid IA - o = o+1 - ODFmap(texID) = o ! remember mapping - Ncomponents(texID) = 1_pInt ! single "component" - Nsym(texID) = 1_pInt ! no symmetry (use full ODF instead) - endif -! adjust multiplicity and number of grains per IP of components - multiplicity(texID) = max(1_pInt,texture_Ngrains(texID)/Ncomponents(texID)/Nsym(texID)) - if (mod(texture_Ngrains(texID),Ncomponents(texID)*Nsym(texID)) /= 0_pInt) then - texture_Ngrains(texID) = multiplicity(texID)*Ncomponents(texID)*Nsym(texID) - write (6,*) 'changed Ngrains to',texture_Ngrains(texID),' for texture',texID - endif -enddo - -!* publish globals -constitutive_maxNgrains = maxval(texture_Ngrains) -constitutive_maxNstatevars = maxval(material_Nslip) + 0_pInt -constitutive_maxNresults = 1_pInt - -!* calc texture_totalNgrains -allocate(texture_totalNgrains(texture_maxN)) ; texture_totalNgrains=0_pInt -do i=1,mesh_NcpElems - texID = mesh_element(4,i) - texture_totalNgrains(texID) = texture_totalNgrains(texID) + FE_Nips(FE_mapElemtype(mesh_element(2,i)))*texture_Ngrains(texID) -enddo - -! generate hybridIA samplings for ODFfile textures to later draw from these populations -allocate(hybridIA_population(3,maxval(texture_totalNgrains,ODFmap /= 0),o)) -do texID = 1,texture_maxN - if (ODFmap(texID) > 0) & - hybridIA_population(:,:,ODFmap(texID)) = IO_hybridIA(texture_totalNgrains(texID),texture_ODFfile(texID)) -enddo - -!* Array allocation -allocate(constitutive_Ngrains(mesh_maxNips,mesh_NcpElems)) ; constitutive_Ngrains=0_pInt -allocate(constitutive_matID(constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_matID=0_pInt -allocate(constitutive_texID(constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_texID=0_pInt -allocate(constitutive_MatVolFrac(constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_MatVolFrac=0.0_pReal -allocate(constitutive_TexVolFrac(constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_TexVolFrac=0.0_pReal -allocate(constitutive_EulerAngles(3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_EulerAngles=0.0_pReal -allocate(constitutive_Nresults(constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_Nresults=0_pInt -allocate(constitutive_results(constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) -constitutive_results=0.0_pReal -allocate(constitutive_Nstatevars(constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_Nstatevars=0_pInt -allocate(constitutive_state_old(constitutive_maxNstatevars,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) -constitutive_state_old=0.0_pReal -allocate(constitutive_state_new(constitutive_maxNstatevars,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) -constitutive_state_new=0.0_pReal -allocate(constitutive_Pforest(constitutive_maxNstatevars,constitutive_maxNstatevars,material_maxN)) -constitutive_Pforest=0.0_pReal -allocate(constitutive_Pparallel(constitutive_maxNstatevars,constitutive_maxNstatevars,material_maxN)) -constitutive_Pparallel=0.0_pReal -allocate(constitutive_rho_p(constitutive_maxNstatevars)) ; constitutive_rho_p=0.0_pReal -allocate(constitutive_rho_f(constitutive_maxNstatevars)) ; constitutive_rho_f=0.0_pReal -allocate(constitutive_rho_m(constitutive_maxNstatevars)) ; constitutive_rho_m=0.0_pReal -allocate(constitutive_passing_stress(constitutive_maxNstatevars)) ; constitutive_passing_stress=0.0_pReal -allocate(constitutive_jump_width(constitutive_maxNstatevars)) ; constitutive_jump_width=0.0_pReal -allocate(constitutive_activation_volume(constitutive_maxNstatevars)) ; constitutive_activation_volume=0.0_pReal -allocate(constitutive_g0_slip(constitutive_maxNstatevars)) ; constitutive_g0_slip=0.0_pReal - -!* Assignment of all grains in all IPs of all cp-elements -do e=1,mesh_NcpElems - matID=mesh_element(3,e) - texID=mesh_element(4,e) - do i=1,FE_Nips(FE_mapElemtype(mesh_element(2,e))) - g = 0_pInt ! grain counter - do m = 1,multiplicity(texID) - o = 0_pInt ! component counter - if (texture_ODFfile(texID)=='') then - do k = 1,texture_nGauss(texID) ! *** gauss *** - o = o+1 - Euler(:,o) = math_sampleGaussOri(texture_Gauss(1:3,k,texID),texture_Gauss(4,k,texID)) - texVolFrac(o) = texture_Gauss(5,k,texID) - enddo - do k = 1,texture_nFiber(texID) ! *** fiber *** - o = o+1 - Euler(:,o) = math_sampleFiberOri(texture_Fiber(1:2,k,texID),texture_Fiber(3:4,k,texID),texture_Fiber(5,k,texID)) - texVolFrac(o) = texture_Fiber(6,k,texID) - enddo - do k = 1,texture_nRandom(texID) ! *** random *** - o = o+1 - Euler(:,o) = math_sampleRandomOri() - texVolfrac(o) = 1.0_pReal-sumVolfrac(texID) - enddo - else ! *** hybrid IA *** - o = 1 ! only singular orientation, i.e. single "component" - Euler(:,o) = hybridIA_population(:,1+sampleCount(texID),ODFmap(texID)) - texVolfrac(o) = 1.0_pReal - endif - if (Nsym(texID) > 1) then ! symmetry generates additional orientations - forall (k=1:o) - Euler(:,1+o+(Nsym(texID)-1)*(k-1):3+o+(Nsym(texID)-1)*(k-1)) = & - math_symmetricEulers(texture_symmetry(texID),Euler(:,k)) - texVolfrac(1+o+(Nsym(texID)-1)*(k-1):3+o+(Nsym(texID)-1)*(k-1)) = texVolfrac(k) - end forall - endif - do s = 1,Nsym(texID)*o ! loop over orientations to be assigned to ip (ex multiplicity) - g = g+1 ! next "grain" - sampleCount(texID) = sampleCount(texID)+1 ! next member of population - constitutive_matID(g,i,e) = matID ! copy matID of element - constitutive_texID(g,i,e) = texID ! copy texID of element - constitutive_MatVolFrac(g,i,e) = 1.0_pReal ! singular material (so far) - constitutive_TexVolFrac(g,i,e) = texVolfrac(s)/multiplicity(texID)/Nsym(texID) - constitutive_Nstatevars(g,i,e) = material_Nslip(matID) ! number of state variables (i.e. tau_c of each slip system) - constitutive_Nresults(g,i,e) = 0 ! number of constitutive results - constitutive_EulerAngles(:,g,i,e) = Euler(:,s) ! store initial orientation - forall (l=1:constitutive_Nstatevars(g,i,e)) ! initialize state variables - constitutive_state_old(l,g,i,e) = material_rho0(matID) - constitutive_state_new(l,g,i,e) = material_rho0(matID) - end forall - enddo ! components - enddo ! multiplicity - enddo ! ip -enddo ! cp_element - - -!* Construction of the hardening matrices -do i=1,material_maxN -!* Iteration over the systems - do j=1,constitutive_maxNstatevars - do k=1,constitutive_maxNstatevars -!* Hardening type * - select case (crystal_SlipIntType(j,k,i)) - case (0) - K_inter=0.0_pReal - case (1) - K_inter=material_SlipIntCoeff(1,i) - case (2) - K_inter=material_SlipIntCoeff(2,i) - case (3) - K_inter=material_SlipIntCoeff(3,i) - case (4) - K_inter=material_SlipIntCoeff(4,i) - case (5) - K_inter=material_SlipIntCoeff(5,i) - case (6) - K_inter=material_SlipIntCoeff(6,i) - end select -!* Projection of the dislocation * - x=dot_product(crystal_sn(:,j,i),crystal_st(:,k,i)) - y=1.0_pReal-x**(2.0_pReal) -!* Interaction matrix * - constitutive_Pforest(j,k,i)=abs(x)*K_inter - if (y>0.0_pReal) then - constitutive_Pparallel(j,k,i)=sqrt(y)*K_inter - else - constitutive_Pparallel(j,k,i)=0.0_pReal - endif - enddo - enddo -enddo - -end subroutine - - -function constitutive_HomogenizedC(ipc,ip,el) -!********************************************************************* -!* This function returns the homogenized elacticity matrix * -!* INPUT: * -!* - ipc : component-ID of current integration point * -!* - ip : current integration point * -!* - el : current element * -!********************************************************************* -use prec, only: pReal,pInt -implicit none - -!* Definition of variables -integer(pInt) ipc,ip,el -real(pReal), dimension(6,6) :: constitutive_homogenizedC - -!* Homogenization scheme -constitutive_homogenizedC=constitutive_MatVolFrac(ipc,ip,el)*material_Cslip_66(:,:,constitutive_matID(ipc,ip,el)) - -return -end function - - -subroutine constitutive_Microstructure(state,Tp,ipc,ip,el) -!********************************************************************* -!* This function calculates from state needed variables * -!* INPUT: * -!* - state : state variables * -!* - Tp : temperature * -!* - ipc : component-ID of current integration point * -!* - ip : current integration point * -!* - el : current element * -!********************************************************************* -use prec, only: pReal,pInt -implicit none - -!* Definition of variables -integer(pInt) ipc,ip,el -integer(pInt) matID,i -real(pReal) Tp -real(pReal), dimension(constitutive_Nstatevars(ipc,ip,el)) :: state - -!* Get the material-ID from the triplet(ipc,ip,el) -matID = constitutive_matID(ipc,ip,el) - -!* Quantities derivated from state -constitutive_rho_f=matmul(constitutive_Pforest(1:constitutive_Nstatevars(ipc,ip,el),& - 1:constitutive_Nstatevars(ipc,ip,el),matID),state) -constitutive_rho_p=matmul(constitutive_Pparallel(1:constitutive_Nstatevars(ipc,ip,el),& - 1:constitutive_Nstatevars(ipc,ip,el),matID),state) -do i=1,material_Nslip(matID) - constitutive_passing_stress(i)=material_tau0(matID)+material_c1(matID)*material_Gmod(matID)*material_bg(matID)*& - sqrt(constitutive_rho_p(i)) - constitutive_jump_width(i)=material_c2(matID)/sqrt(constitutive_rho_f(i)) - constitutive_activation_volume(i)=material_c3(matID)*constitutive_jump_width(i)*material_bg(matID)**2.0_pReal - constitutive_rho_m(i)=(2.0_pReal*Kb*Tp*sqrt(constitutive_rho_p(i)))/& - (material_c1(matID)*material_c3(matID)*material_Gmod(matID)*constitutive_jump_width(i)*material_bg(matID)**3.0_pReal) - constitutive_g0_slip(i)=constitutive_rho_m(i)*material_bg(matID)*attack_frequency*constitutive_jump_width(i)*& - exp(-(material_Qedge(matID)+constitutive_passing_stress(i)*constitutive_activation_volume(i))/& - (Kb*Tp)) -enddo - -end subroutine - - -subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,state,Tp,ipc,ip,el) -!********************************************************************* -!* This subroutine contains the constitutive equation for * -!* calculating the velocity gradient * -!* INPUT: * -!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * -!* - state : current microstructure * -!* - Tp : temperature * -!* - ipc : component-ID of current integration point * -!* - ip : current integration point * -!* - el : current element * -!* OUTPUT: * -!* - Lp : plastic velocity gradient * -!* - dLp_dTstar : derivative of Lp (4th-order tensor) * -!********************************************************************* -use prec, only: pReal,pInt -use crystal, only: crystal_Sslip,crystal_Sslip_v -implicit none - -!* Definition of variables -integer(pInt) ipc,ip,el -integer(pInt) matID,i,k,l,m,n -real(pReal) Tp -real(pReal), dimension(6) :: Tstar_v -real(pReal), dimension(3,3) :: Lp -real(pReal), dimension(3,3,3,3) :: dLp_dTstar -real(pReal), dimension(constitutive_Nstatevars(ipc,ip,el)) :: state,gdot_slip,dgdot_dtauslip,tau_slip - -!* Get the material-ID from the triplet(ipc,ip,el) -matID = constitutive_matID(ipc,ip,el) - -!* Calculation of Lp -Lp = 0.0_pReal -do i=1,material_Nslip(matID) - tau_slip(i)=dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID))) - gdot_slip(i)=constitutive_g0_slip(i)*sinh((abs(tau_slip(i))*constitutive_activation_volume(i))/(Kb*Tp))*& - sign(1.0_pReal,tau_slip(i)) - Lp=Lp+gdot_slip(i)*crystal_Sslip(:,:,i,material_CrystalStructure(matID)) -enddo - -!* Calculation of the tangent of Lp -dLp_dTstar=0.0_pReal -do i=1,material_Nslip(matID) - dgdot_dtauslip(i)=((constitutive_g0_slip(i)*constitutive_activation_volume(i))/(Kb*Tp))*& - cosh((abs(tau_slip(i))*constitutive_activation_volume(i))/(Kb*Tp)) - forall (k=1:3,l=1:3,m=1:3,n=1:3) - dLp_dTstar(k,l,m,n) = dLp_dTstar(k,l,m,n)+ & - dgdot_dtauslip(i)*crystal_Sslip(k,l,i,material_CrystalStructure(matID))* & - (crystal_Sslip(m,n,i,material_CrystalStructure(matID))+ & - crystal_Sslip(n,m,i,material_CrystalStructure(matID)))/2.0_pReal ! force m,n symmetry - endforall -enddo - -return -end subroutine - - -function constitutive_dotState(Tstar_v,state,Tp,ipc,ip,el) -!********************************************************************* -!* This subroutine contains the constitutive equation for * -!* calculating the velocity gradient * -!* INPUT: * -!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * -!* - state : current microstructure * -!* - Tp : temperature * -!* - ipc : component-ID of current integration point * -!* - ip : current integration point * -!* - el : current element * -!* OUTPUT: * -!* - constitutive_DotState : evolution of state variable * -!********************************************************************* -use prec, only: pReal,pInt -use crystal, only: crystal_Sslip_v -implicit none - -!* Definition of variables -integer(pInt) ipc,ip,el -integer(pInt) matID,i -real(pReal) Tp,tau_slip,gdot_slip -real(pReal), dimension(6) :: Tstar_v -real(pReal), dimension(constitutive_Nstatevars(ipc,ip,el)) :: constitutive_dotState,state,lock,recovery - -!* Get the material-ID from the triplet(ipc,ip,el) -matID = constitutive_matID(ipc,ip,el) - -!* Hardening of each system -do i=1,constitutive_Nstatevars(ipc,ip,el) - tau_slip = dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID))) - gdot_slip = constitutive_g0_slip(i)*sinh((abs(tau_slip)*constitutive_activation_volume(i))/(Kb*Tp))*& - sign(1.0_pReal,tau_slip) - if (abs(tau_slip)>1.0e-20_pReal) then - lock(i)=(material_c4(matID)*sqrt(constitutive_rho_f(i))*abs(gdot_slip))/material_bg(matID) - recovery(i)=material_c5(matID)*state(i)*abs(gdot_slip) - constitutive_dotState(i)=lock(i)-recovery(i) - else - constitutive_dotState(i)=0.0_pReal - endif -enddo - -return -end function - - -function constitutive_post_results(Tstar_v,state,dt,Tp,ipc,ip,el) -!********************************************************************* -!* return array of constitutive results * -!* INPUT: * -!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * -!* - state : current microstructure * -!* - dt : current time increment * -!* - Tp : temperature * -!* - ipc : component-ID of current integration point * -!* - ip : current integration point * -!* - el : current element * -!********************************************************************* -use prec, only: pReal,pInt -use crystal, only: crystal_Sslip_v -implicit none - -!* Definition of variables -integer(pInt) ipc,ip,el -integer(pInt) matID,i -real(pReal) dt,Tp,tau_slip -real(pReal), dimension(6) :: Tstar_v -real(pReal), dimension(constitutive_Nstatevars(ipc,ip,el)) :: state -real(pReal), dimension(constitutive_Nresults(ipc,ip,el)) :: constitutive_post_results - -!* Get the material-ID from the triplet(ipc,ip,el) -matID = constitutive_matID(ipc,ip,el) - -if(constitutive_Nresults(ipc,ip,el)==0) return - -do i=1,material_Nslip(matID) - constitutive_post_results(i) = state(i) - tau_slip=dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID))) - constitutive_post_results(i+material_Nslip(matID)) = & - dt*constitutive_g0_slip(i)*sinh((abs(tau_slip)*constitutive_activation_volume(i))/(Kb*Tp))*& - sign(1.0_pReal,tau_slip) -enddo -return - -end function - -END MODULE \ No newline at end of file diff --git a/dislocation_based_subroutine/crystal.f90 b/dislocation_based_subroutine/crystal.f90 deleted file mode 100644 index 54aeafb39..000000000 --- a/dislocation_based_subroutine/crystal.f90 +++ /dev/null @@ -1,282 +0,0 @@ - -!************************************ -!* Module: CRYSTAL * -!************************************ -!* contains: * -!* - Crystal structure definition * -!* - Slip system definition * -!* - Schmid matrices calculation * -!* - Hardening matrices calculation * -!************************************ - -MODULE crystal - -!*** Include other modules *** -use prec, only: pReal,pInt -implicit none - -!************************************ -!* Crystal structures * -!************************************ -!* Number of crystal structures (1-FCC,2-BCC,3-HCP) -integer(pInt), parameter :: crystal_MaxCrystalStructure = 3 -!* Total number of slip systems per crystal structure -!* (as to be changed according the definition of slip systems) -integer(pInt), dimension(crystal_MaxCrystalStructure), parameter :: crystal_MaxNslipOfStructure = & -reshape((/12,48,12/),(/crystal_MaxCrystalStructure/)) -!* Maximum number of slip systems over crystal structures -integer(pInt), parameter :: crystal_MaxMaxNslipOfStructure = 48 -!* Slip direction, slip normales and Schmid matrices -real(pReal), dimension(3,3,crystal_MaxMaxNslipOfStructure,crystal_MaxCrystalStructure) :: crystal_Sslip -real(pReal), dimension(6,crystal_MaxMaxNslipOfStructure,crystal_MaxCrystalStructure) :: crystal_Sslip_v -real(pReal), dimension(3,crystal_MaxMaxNslipOfStructure,crystal_MaxCrystalStructure) :: crystal_sn -real(pReal), dimension(3,crystal_MaxMaxNslipOfStructure,crystal_MaxCrystalStructure) :: crystal_sd -real(pReal), dimension(3,crystal_MaxMaxNslipOfStructure,crystal_MaxCrystalStructure) :: crystal_st -!* Slip_slip interaction matrices -integer(pInt), dimension(crystal_MaxMaxNslipOfStructure,crystal_MaxMaxNslipOfStructure,crystal_MaxCrystalStructure) :: & -crystal_SlipIntType - -!*** Slip systems for FCC structures (1) *** -!* System {111}<110> Sort according Eisenlohr&Hantcherli -data crystal_sd(:, 1,1)/ 0, 1,-1/ ; data crystal_sn(:, 1,1)/ 1, 1, 1/ -data crystal_sd(:, 2,1)/-1, 0, 1/ ; data crystal_sn(:, 2,1)/ 1, 1, 1/ -data crystal_sd(:, 3,1)/ 1,-1, 0/ ; data crystal_sn(:, 3,1)/ 1, 1, 1/ -data crystal_sd(:, 4,1)/ 0,-1,-1/ ; data crystal_sn(:, 4,1)/-1,-1, 1/ -data crystal_sd(:, 5,1)/ 1, 0, 1/ ; data crystal_sn(:, 5,1)/-1,-1, 1/ -data crystal_sd(:, 6,1)/-1, 1, 0/ ; data crystal_sn(:, 6,1)/-1,-1, 1/ -data crystal_sd(:, 7,1)/ 0,-1, 1/ ; data crystal_sn(:, 7,1)/ 1,-1,-1/ -data crystal_sd(:, 8,1)/-1, 0,-1/ ; data crystal_sn(:, 8,1)/ 1,-1,-1/ -data crystal_sd(:, 9,1)/ 1, 1, 0/ ; data crystal_sn(:, 9,1)/ 1,-1,-1/ -data crystal_sd(:,10,1)/ 0, 1, 1/ ; data crystal_sn(:,10,1)/-1, 1,-1/ -data crystal_sd(:,11,1)/ 1, 0,-1/ ; data crystal_sn(:,11,1)/-1, 1,-1/ -data crystal_sd(:,12,1)/-1,-1, 0/ ; data crystal_sn(:,12,1)/-1, 1,-1/ - -!*** Slip-Slip interactions for FCC structures (1) *** -data crystal_SlipIntType( 1,:,1)/1,2,2,4,6,5,3,5,5,4,5,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 2,:,1)/2,1,2,6,4,5,5,4,6,5,3,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 3,:,1)/2,2,1,5,5,3,5,6,4,6,5,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 4,:,1)/4,6,5,1,2,2,4,5,6,3,5,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 5,:,1)/6,4,5,2,1,2,5,3,5,5,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 6,:,1)/5,5,3,2,2,1,6,5,4,5,6,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 7,:,1)/3,5,5,4,5,6,1,2,2,4,6,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 8,:,1)/5,4,6,5,3,5,2,1,2,6,4,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 9,:,1)/5,6,4,6,5,4,2,2,1,5,5,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType(10,:,1)/4,5,6,3,5,5,4,6,5,1,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType(11,:,1)/5,3,5,5,4,6,6,4,5,2,1,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType(12,:,1)/6,5,4,5,6,4,5,5,3,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ - -!*** Slip systems for BCC structures (2) *** -!* System {110}<111> -!* Sort? -data crystal_sd(:, 1,2)/ 1,-1, 1/ ; data crystal_sn(:, 1,2)/ 0, 1, 1/ -data crystal_sd(:, 2,2)/-1,-1, 1/ ; data crystal_sn(:, 2,2)/ 0, 1, 1/ -data crystal_sd(:, 3,2)/ 1, 1, 1/ ; data crystal_sn(:, 3,2)/ 0,-1, 1/ -data crystal_sd(:, 4,2)/-1, 1, 1/ ; data crystal_sn(:, 4,2)/ 0,-1, 1/ -data crystal_sd(:, 5,2)/-1, 1, 1/ ; data crystal_sn(:, 5,2)/ 1, 0, 1/ -data crystal_sd(:, 6,2)/-1,-1, 1/ ; data crystal_sn(:, 6,2)/ 1, 0, 1/ -data crystal_sd(:, 7,2)/ 1, 1, 1/ ; data crystal_sn(:, 7,2)/-1, 0, 1/ -data crystal_sd(:, 8,2)/ 1,-1, 1/ ; data crystal_sn(:, 8,2)/-1, 0, 1/ -data crystal_sd(:, 9,2)/-1, 1, 1/ ; data crystal_sn(:, 9,2)/ 1, 1, 0/ -data crystal_sd(:,10,2)/-1, 1,-1/ ; data crystal_sn(:,10,2)/ 1, 1, 0/ -data crystal_sd(:,11,2)/ 1, 1, 1/ ; data crystal_sn(:,11,2)/-1, 1, 0/ -data crystal_sd(:,12,2)/ 1, 1,-1/ ; data crystal_sn(:,12,2)/-1, 1, 0/ -!* System {112}<111> -!* Sort? -data crystal_sd(:,13,2)/-1, 1, 1/ ; data crystal_sn(:,13,2)/ 2, 1, 1/ -data crystal_sd(:,14,2)/ 1, 1, 1/ ; data crystal_sn(:,14,2)/-2, 1, 1/ -data crystal_sd(:,15,2)/ 1, 1,-1/ ; data crystal_sn(:,15,2)/ 2,-1, 1/ -data crystal_sd(:,16,2)/ 1,-1, 1/ ; data crystal_sn(:,16,2)/ 2, 1,-1/ -data crystal_sd(:,17,2)/ 1,-1, 1/ ; data crystal_sn(:,17,2)/ 1, 2, 1/ -data crystal_sd(:,18,2)/ 1, 1,-1/ ; data crystal_sn(:,18,2)/-1, 2, 1/ -data crystal_sd(:,19,2)/ 1, 1, 1/ ; data crystal_sn(:,19,2)/ 1,-2, 1/ -data crystal_sd(:,20,2)/-1, 1, 1/ ; data crystal_sn(:,20,2)/ 1, 2,-1/ -data crystal_sd(:,21,2)/ 1, 1,-1/ ; data crystal_sn(:,21,2)/ 1, 1, 2/ -data crystal_sd(:,22,2)/ 1,-1, 1/ ; data crystal_sn(:,22,2)/-1, 1, 2/ -data crystal_sd(:,23,2)/-1, 1, 1/ ; data crystal_sn(:,23,2)/ 1,-1, 2/ -data crystal_sd(:,24,2)/ 1, 1, 1/ ; data crystal_sn(:,24,2)/ 1, 1,-2/ -!* System {123}<111> -!* Sort? -data crystal_sd(:,25,2)/ 1, 1,-1/ ; data crystal_sn(:,25,2)/ 1, 2, 3/ -data crystal_sd(:,26,2)/ 1,-1, 1/ ; data crystal_sn(:,26,2)/-1, 2, 3/ -data crystal_sd(:,27,2)/-1, 1, 1/ ; data crystal_sn(:,27,2)/ 1,-2, 3/ -data crystal_sd(:,28,2)/ 1, 1, 1/ ; data crystal_sn(:,28,2)/ 1, 2,-3/ -data crystal_sd(:,29,2)/ 1,-1, 1/ ; data crystal_sn(:,29,2)/ 1, 3, 2/ -data crystal_sd(:,30,2)/ 1, 1,-1/ ; data crystal_sn(:,30,2)/-1, 3, 2/ -data crystal_sd(:,31,2)/ 1, 1, 1/ ; data crystal_sn(:,31,2)/ 1,-3, 2/ -data crystal_sd(:,32,2)/-1, 1, 1/ ; data crystal_sn(:,32,2)/ 1, 3,-2/ -data crystal_sd(:,33,2)/ 1, 1,-1/ ; data crystal_sn(:,33,2)/ 2, 1, 3/ -data crystal_sd(:,34,2)/ 1,-1, 1/ ; data crystal_sn(:,34,2)/-2, 1, 3/ -data crystal_sd(:,35,2)/-1, 1, 1/ ; data crystal_sn(:,35,2)/ 2,-1, 3/ -data crystal_sd(:,36,2)/ 1, 1, 1/ ; data crystal_sn(:,36,2)/ 2, 1,-3/ -data crystal_sd(:,37,2)/ 1,-1, 1/ ; data crystal_sn(:,37,2)/ 2, 3, 1/ -data crystal_sd(:,38,2)/ 1, 1,-1/ ; data crystal_sn(:,38,2)/-2, 3, 1/ -data crystal_sd(:,39,2)/ 1, 1, 1/ ; data crystal_sn(:,39,2)/ 2,-3, 1/ -data crystal_sd(:,40,2)/-1, 1, 1/ ; data crystal_sn(:,40,2)/ 2, 3,-1/ -data crystal_sd(:,41,2)/-1, 1, 1/ ; data crystal_sn(:,41,2)/ 3, 1, 2/ -data crystal_sd(:,42,2)/ 1, 1, 1/ ; data crystal_sn(:,42,2)/-3, 1, 2/ -data crystal_sd(:,43,2)/ 1, 1,-1/ ; data crystal_sn(:,43,2)/ 3,-1, 2/ -data crystal_sd(:,44,2)/ 1,-1, 1/ ; data crystal_sn(:,44,2)/ 3, 1,-2/ -data crystal_sd(:,45,2)/-1, 1, 1/ ; data crystal_sn(:,45,2)/ 3, 2, 1/ -data crystal_sd(:,46,2)/ 1, 1, 1/ ; data crystal_sn(:,46,2)/-3, 2, 1/ -data crystal_sd(:,47,2)/ 1, 1,-1/ ; data crystal_sn(:,47,2)/ 3,-2, 1/ -data crystal_sd(:,48,2)/ 1,-1, 1/ ; data crystal_sn(:,48,2)/ 3, 2,-1/ - -!*** Slip-Slip interactions for BCC structures (2) *** -data crystal_SlipIntType( 1,:,2)/1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType( 2,:,2)/2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType( 3,:,2)/2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType( 4,:,2)/2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType( 5,:,2)/2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType( 6,:,2)/2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType( 7,:,2)/2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType( 8,:,2)/2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType( 9,:,2)/2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(10,:,2)/2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(11,:,2)/2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(12,:,2)/2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(13,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(14,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(15,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(16,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(17,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(18,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(19,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(20,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(21,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(22,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(23,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(24,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(25,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(26,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(27,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(28,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(29,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(30,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(31,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(32,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(33,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(34,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(35,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(36,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(37,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(38,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(39,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(40,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2/ -data crystal_SlipIntType(41,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2/ -data crystal_SlipIntType(42,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2/ -data crystal_SlipIntType(43,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2/ -data crystal_SlipIntType(44,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2/ -data crystal_SlipIntType(45,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2/ -data crystal_SlipIntType(46,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2/ -data crystal_SlipIntType(47,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2/ -data crystal_SlipIntType(48,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1/ - -!*** Slip systems for HCP structures (3) *** -!* Basal systems {0001}<1120> (independent of c/a-ratio) -!* 1- (0 0 0 1)[-2 1 1 0] -!* 2- (0 0 0 1)[ 1 -2 1 0] -!* 3- (0 0 0 1)[ 1 1 -2 0] -!* Plane (hkil)->(hkl) -!* Direction [uvtw]->[(u-t) (v-t) w] -!* Automatical transformation from Bravais to Miller -!* not done for the moment -!* Sort? -data crystal_sd(:, 1,3)/-1, 0, 0/ ; data crystal_sn(:, 1,3)/ 0, 0, 1/ -data crystal_sd(:, 2,3)/ 0,-1, 0/ ; data crystal_sn(:, 2,3)/ 0, 0, 1/ -data crystal_sd(:, 3,3)/ 1, 1, 0/ ; data crystal_sn(:, 3,3)/ 0, 0, 1/ -!* 1st type prismatic systems {1010}<1120> (independent of c/a-ratio) -!* 1- ( 0 1 -1 0)[-2 1 1 0] -!* 2- ( 1 0 -1 0)[ 1 -2 1 0] -!* 3- (-1 1 0 0)[ 1 1 -2 0] -!* Sort? -data crystal_sd(:, 4,3)/-1, 0, 0/ ; data crystal_sn(:, 4,3)/ 0, 1, 0/ -data crystal_sd(:, 5,3)/ 0,-1, 0/ ; data crystal_sn(:, 5,3)/ 1, 0, 0/ -data crystal_sd(:, 6,3)/ 1, 1, 0/ ; data crystal_sn(:, 6,3)/-1, 1, 0/ -!* 1st type 1st order pyramidal systems {1011}<1120> -!* plane normales depend on the c/a-ratio -!* 1- ( 0 -1 1 1)[-2 1 1 0] -!* 2- ( 0 1 -1 1)[-2 1 1 0] -!* 3- (-1 0 1 1)[ 1 -2 1 0] -!* 4- ( 1 0 -1 1)[ 1 -2 1 0] -!* 5- (-1 1 0 1)[ 1 1 -2 0] -!* 6- ( 1 -1 0 1)[ 1 1 -2 0] -!* Sort? -data crystal_sd(:, 7,3)/-1, 0, 0/ ; data crystal_sn(:, 7,3)/ 0,-1, 1/ -data crystal_sd(:, 8,3)/ 0,-1, 0/ ; data crystal_sn(:, 8,3)/ 0, 1, 1/ -data crystal_sd(:, 9,3)/ 1, 1, 0/ ; data crystal_sn(:, 9,3)/-1, 0, 1/ -data crystal_sd(:,10,3)/-1, 0, 0/ ; data crystal_sn(:,10,3)/ 1, 0, 1/ -data crystal_sd(:,11,3)/ 0,-1, 0/ ; data crystal_sn(:,11,3)/-1, 1, 1/ -data crystal_sd(:,12,3)/ 1, 1, 0/ ; data crystal_sn(:,12,3)/ 1,-1, 1/ - -!*** Slip-Slip interactions for HCP structures (3) *** -data crystal_SlipIntType( 1,:,3)/1,2,2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 2,:,3)/2,1,2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 3,:,3)/2,2,1,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 4,:,3)/2,2,2,1,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 5,:,3)/2,2,2,2,1,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 6,:,3)/2,2,2,2,2,1,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 7,:,3)/2,2,2,2,2,2,1,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 8,:,3)/2,2,2,2,2,2,2,1,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType( 9,:,3)/2,2,2,2,2,2,2,2,1,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType(10,:,3)/2,2,2,2,2,2,2,2,2,1,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType(11,:,3)/2,2,2,2,2,2,2,2,2,2,1,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ -data crystal_SlipIntType(12,:,3)/2,2,2,2,2,2,2,2,2,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ - -CONTAINS -!**************************************** -!* - crystal_Init -!* - crystal_SchmidMatrices -!* - crystal_HardeningMatrices -!**************************************** - - -subroutine crystal_Init() -!************************************** -!* Module initialization * -!************************************** -call crystal_SchmidMatrices() -end subroutine - - -subroutine crystal_SchmidMatrices() -!************************************** -!* Calculation of Schmid matrices * -!************************************** -use prec, only: pReal,pInt -implicit none - -!* Definition of variables -integer(pInt) i,j,k,l -real(pReal) invNorm - -!* Iteration over the crystal structures -do l=1,3 -!* Iteration over the systems - do k=1,crystal_MaxNslipOfStructure(l) -!* Definition of transverse direction st for the frame (sd,st,sn) - crystal_st(1,k,l)=crystal_sn(2,k,l)*crystal_sd(3,k,l)-crystal_sn(3,k,l)*crystal_sd(2,k,l) - crystal_st(2,k,l)=crystal_sn(3,k,l)*crystal_sd(1,k,l)-crystal_sn(1,k,l)*crystal_sd(3,k,l) - crystal_st(3,k,l)=crystal_sn(1,k,l)*crystal_sd(2,k,l)-crystal_sn(2,k,l)*crystal_sd(1,k,l) -!* Defintion of Schmid matrix - forall (i=1:3,j=1:3) - crystal_Sslip(i,j,k,l)=crystal_sd(i,k,l)*crystal_sn(j,k,l) - endforall -!* Normalization of Schmid matrix - invNorm=dsqrt(1.0_pReal/((crystal_sn(1,k,l)**2+crystal_sn(2,k,l)**2+crystal_sn(3,k,l)**2)*& - (crystal_sd(1,k,l)**2+crystal_sd(2,k,l)**2+crystal_sd(3,k,l)**2))) - crystal_Sslip(:,:,k,l)=crystal_Sslip(:,:,k,l)*invNorm -!* Vectorization of normalized Schmid matrix - crystal_Sslip_v(1,k,l)=crystal_Sslip(1,1,k,l) - crystal_Sslip_v(2,k,l)=crystal_Sslip(2,2,k,l) - crystal_Sslip_v(3,k,l)=crystal_Sslip(3,3,k,l) - !* be compatible with Mandel notation of Tstar - crystal_Sslip_v(4,k,l)=(crystal_Sslip(1,2,k,l)+crystal_Sslip(2,1,k,l))/dsqrt(2.0_pReal) - crystal_Sslip_v(5,k,l)=(crystal_Sslip(2,3,k,l)+crystal_Sslip(3,2,k,l))/dsqrt(2.0_pReal) - crystal_Sslip_v(6,k,l)=(crystal_Sslip(1,3,k,l)+crystal_Sslip(3,1,k,l))/dsqrt(2.0_pReal) - enddo -enddo - -end subroutine - - -END MODULE diff --git a/dislocation_based_subroutine/math.f90 b/dislocation_based_subroutine/math.f90 deleted file mode 100644 index 1a7314432..000000000 --- a/dislocation_based_subroutine/math.f90 +++ /dev/null @@ -1,1730 +0,0 @@ - -!############################################################## - MODULE math -!############################################################## - - use prec, only: pReal,pInt - implicit none - - real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal - real(pReal), parameter :: inDeg = 180.0_pReal/pi - real(pReal), parameter :: inRad = pi/180.0_pReal -! *** 3x3 Identity *** - real(pReal), dimension(3,3), parameter :: math_I3 = & - reshape( (/ & - 1.0_pReal,0.0_pReal,0.0_pReal, & - 0.0_pReal,1.0_pReal,0.0_pReal, & - 0.0_pReal,0.0_pReal,1.0_pReal /),(/3,3/)) -! *** Mandel notation *** - integer(pInt), dimension (2,6), parameter :: mapMandel = & - reshape((/& - 1,1, & - 2,2, & - 3,3, & - 1,2, & - 2,3, & - 1,3 & - /),(/2,6/)) - real(pReal), dimension(6), parameter :: nrmMandel = & - (/1.0_pReal,1.0_pReal,1.0_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal/) - real(pReal), dimension(6), parameter :: invnrmMandel = & - (/1.0_pReal,1.0_pReal,1.0_pReal,0.7071067811865476_pReal,0.7071067811865476_pReal,0.7071067811865476_pReal/) -! *** Voigt notation *** - integer(pInt), dimension (2,6), parameter :: mapVoigt = & - reshape((/& - 1,1, & - 2,2, & - 3,3, & - 2,3, & - 1,3, & - 1,2 & - /),(/2,6/)) - real(pReal), dimension(6), parameter :: nrmVoigt = & - (/1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal/) - real(pReal), dimension(6), parameter :: invnrmVoigt = & - (/1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal/) - - - CONTAINS - -!************************************************************************** -! initialization of module -!************************************************************************** - SUBROUTINE math_init () - - use prec, only: pReal,pInt - implicit none - - integer (pInt) seed - - call random_seed() - call get_seed(seed) - call halton_seed_set(seed) - call halton_ndim_set(3) - - END SUBROUTINE - -!************************************************************************** -! Quicksort algorithm for two-dimensional integer arrays -! -! Sorting is done with respect to array(1,:) -! and keeps array(2:N,:) linked to it. -!************************************************************************** - RECURSIVE SUBROUTINE qsort(a, istart, iend) - - implicit none - integer(pInt), dimension(:,:) :: a - integer(pInt) :: istart,iend,ipivot - - if (istart < iend) then - ipivot = math_partition(a,istart, iend) - call qsort(a, istart, ipivot-1) - call qsort(a, ipivot+1, iend) - endif - return - - END SUBROUTINE - -!************************************************************************** -! Partitioning required for quicksort -!************************************************************************** - integer(pInt) FUNCTION math_partition(a, istart, iend) - - implicit none - integer(pInt), dimension(:,:) :: a - integer(pInt) :: istart,iend,d,i,j,k,x,tmp - - d = size(a,1) ! number of linked data -! set the starting and ending points, and the pivot point - i = istart - j = iend - x = a(1,istart) - do -! find the first element on the right side less than or equal to the pivot point - do j = j, istart, -1 - if (a(1,j) <= x) exit - enddo -! find the first element on the left side greater than the pivot point - do i = i, iend - if (a(1,i) > x) exit - enddo - if (i < j ) then ! if the indexes do not cross, exchange values - do k = 1,d - tmp = a(k,i) - a(k,i) = a(k,j) - a(k,j) = tmp - enddo - else ! if they do cross, exchange left value with pivot and return with the partition index - do k = 1,d - tmp = a(k,istart) - a(k,istart) = a(k,j) - a(k,j) = tmp - enddo - math_partition = j - return - endif - enddo - - END FUNCTION - - -!************************************************************************** -! second rank identity tensor of specified dimension -!************************************************************************** - FUNCTION math_identity2nd(dimen) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i,dimen - real(pReal), dimension(dimen,dimen) :: math_identity2nd - - math_identity2nd = 0.0_pReal - forall (i=1:dimen) math_identity2nd(i,i) = 1.0_pReal - return - - END FUNCTION - - -!************************************************************************** -! fourth rank identity tensor of specified dimension -!************************************************************************** - FUNCTION math_identity4th(dimen) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i,j,k,l,dimen - real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th - - forall (i=1:dimen,j=1:dimen,k=1:dimen,l=1:dimen) math_identity4th(i,j,k,l) = & - 0.5_pReal*(math_I3(i,k)*math_I3(j,k)+math_I3(i,l)*math_I3(j,k)) - return - - END FUNCTION - - -!************************************************************************** -! Cramer inversion of 3x3 matrix -!************************************************************************** - SUBROUTINE math_invert3x3(A, InvA, DetA, error) - -! Bestimmung der Determinanten und Inversen einer 3x3-Matrix - -! A = Matrix A -! InvA = Inverse of A -! DetA = Determinant of A -! error = logical - - use prec, only: pReal,pInt - implicit none - - logical, intent(out) :: error - real(pReal),dimension(3,3),intent(in) :: A - real(pReal),dimension(3,3),intent(out) :: InvA - real(pReal), intent(out) :: DetA - - DetA = A(1,1) * ( A(2,2) * A(3,3) - A(2,3) * A(3,2) )& - - A(1,2) * ( A(2,1) * A(3,3) - A(2,3) * A(3,1) )& - + A(1,3) * ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) - - if (DetA <= 0.0000001) then - error = .true. - else - InvA(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2) ) / DetA - InvA(2,1) = ( -A(2,1) * A(3,3) + A(2,3) * A(3,1) ) / DetA - InvA(3,1) = ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) / DetA - - InvA(1,2) = ( -A(1,2) * A(3,3) + A(1,3) * A(3,2) ) / DetA - InvA(2,2) = ( A(1,1) * A(3,3) - A(1,3) * A(3,1) ) / DetA - InvA(3,2) = ( -A(1,1) * A(3,2) + A(1,2) * A(3,1) ) / DetA - - InvA(1,3) = ( A(1,2) * A(2,3) - A(1,3) * A(2,2) ) / DetA - InvA(2,3) = ( -A(1,1) * A(2,3) + A(1,3) * A(2,1) ) / DetA - InvA(3,3) = ( A(1,1) * A(2,2) - A(1,2) * A(2,1) ) / DetA - - error = .false. - endif - return - - END SUBROUTINE - - -!************************************************************************** -! Gauss elimination to invert 6x6 matrix -!************************************************************************** - SUBROUTINE math_invert6x6(A, InvA, AnzNegEW, error) - -! Invertieren einer 6x6-Matrix - -! A = Matrix A -! InvA = Inverse von A -! AnzNegEW = Anzahl der negativen Eigenwerte von A -! error = logical -! = false: Inversion wurde durchgefuehrt. -! = true: Die Inversion in SymGauss wurde wegen eines verschwindenen -! Pivotelement abgebrochen. - - use prec, only: pReal,pInt - implicit none - - real(pReal),dimension(6,6), intent(in) :: A - real(pReal),dimension(6,6), intent(out) :: InvA - integer(pInt), intent(out) :: AnzNegEW - logical, intent(out) :: error - - integer(pInt) i - real(pReal) LogAbsDetA - real(pReal),dimension(6,6) :: B - - InvA = 0.0_pReal - - forall(i = 1:6) InvA(i,i) = 1.0_pReal - B = A - CALL Gauss (B,InvA,LogAbsDetA,AnzNegEW,error) - - RETURN - - END SUBROUTINE math_invert6x6 - -! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - SUBROUTINE Gauss (A,B,LogAbsDetA,NegHDK,error) - -! Loesung eines linearen Gleichungsssystem A * X = B mit Hilfe des -! GAUSS-Algorithmusses - -! Zur numerischen Stabilisierung wird eine Zeilen- und Spaltenpivotsuche -! durchgefuehrt. - -! Eingabeparameter: -! -! A(6,6) = Koeffizientenmatrix A -! B(6,6) = rechte Seiten B -! -! Ausgabeparameter: -! -! B(6,6) = Matrix der Unbekanntenvektoren X -! LogAbsDetA = 10-Logarithmus des Betrages der Determinanten von A -! NegHDK = Anzahl der negativen Hauptdiagonalkoeffizienten nach der -! Vorwaertszerlegung -! error = logical -! = false: Das Gleichungssystem wurde geloest. -! = true : Matrix A ist singulaer. - -! A und B werden veraendert! - - use prec, only: pReal,pInt - implicit none - - logical error - integer (pInt) NegHDK - real(pReal) LogAbsDetA - real(pReal) A(6,6), B(6,6) - - INTENT (OUT) LogAbsDetA, NegHDK, error - INTENT (INOUT) A, B - - LOGICAL SortX - integer (pInt) PivotZeile, PivotSpalte, StoreI, I, IP1, J, K, L - integer (pInt) XNr(6) - real(pReal) AbsA, PivotWert, EpsAbs, Quote - real(pReal) StoreA(6), StoreB(6) - - error = .true. - NegHDK = 1 - SortX = .FALSE. - -! Unbekanntennumerierung - - DO I = 1, 6 - XNr(I) = I - END DO - -! Genauigkeitsschranke und Bestimmung des groessten Pivotelementes - - PivotWert = ABS(A(1,1)) - PivotZeile = 1 - PivotSpalte = 1 - - DO I = 1, 6 - DO J = 1, 6 - AbsA = ABS(A(I,J)) - IF (AbsA .GT. PivotWert) THEN - PivotWert = AbsA - PivotZeile = I - PivotSpalte = J - END IF - END DO - END DO - - IF (PivotWert .LT. 0.0000001) RETURN ! Pivotelement = 0? - - EpsAbs = PivotWert * 0.1_pReal ** PRECISION(1.0_pReal) - -! V O R W A E R T S T R I A N G U L A T I O N - - DO I = 1, 6 - 1 - -! Zeilentausch? - - IF (PivotZeile .NE. I) THEN - - StoreA(I:6) = A(I,I:6) - A(I,I:6) = A(PivotZeile,I:6) - A(PivotZeile,I:6) = StoreA(I:6) - - StoreB(1:6) = B(I,1:6) - B(I,1:6) = B(PivotZeile,1:6) - B(PivotZeile,1:6) = StoreB(1:6) - - SortX = .TRUE. - - END IF - -! Spaltentausch? - - IF (PivotSpalte .NE. I) THEN - - StoreA(1:6) = A(1:6,I) - A(1:6,I) = A(1:6,PivotSpalte) - A(1:6,PivotSpalte) = StoreA(1:6) - - StoreI = XNr(I) - XNr(I) = XNr(PivotSpalte) - XNr(PivotSpalte) = StoreI - - SortX = .TRUE. - - END IF - -! Triangulation - - DO J = I + 1, 6 - Quote = A(J,I) / A(I,I) - DO K = I + 1, 6 - A(J,K) = A(J,K) - Quote * A(I,K) - END DO - DO K = 1, 6 - B(J,K) = B(J,K) - Quote * B(I,K) - END DO - END DO - -! Bestimmung des groessten Pivotelementes - - IP1 = I + 1 - PivotWert = ABS(A(IP1,IP1)) - PivotZeile = IP1 - PivotSpalte = IP1 - - DO J = IP1, 6 - DO K = IP1, 6 - AbsA = ABS(A(J,K)) - IF (AbsA .GT. PivotWert) THEN - PivotWert = AbsA - PivotZeile = J - PivotSpalte = K - END IF - END DO - END DO - - IF (PivotWert .LT. EpsAbs) RETURN ! Pivotelement = 0? - - END DO - -! R U E C K W A E R T S A U F L O E S U N G - - DO I = 6, 1, -1 - DO L = 1, 6 - DO J = I + 1, 6 - B(I,L) = B(I,L) - A(I,J) * B(J,L) - END DO - B(I,L) = B(I,L) / A(I,I) - END DO - END DO - -! Sortieren der Unbekanntenvektoren? - - IF (SortX) THEN - - DO L = 1, 6 - StoreA(1:6) = B(1:6,L) - DO I = 1, 6 - J = XNr(I) - B(J,L) = StoreA(I) - END DO - END DO - END IF - -! Determinante - - LogAbsDetA = 0.0_pReal - NegHDK = 0 - - DO I = 1, 6 - IF (A(I,I) .LT. 0.0_pReal) NegHDK = NegHDK + 1 - AbsA = ABS(A(I,I)) - LogAbsDetA = LogAbsDetA + LOG10(AbsA) - END DO - - error = .false. - - RETURN - - END SUBROUTINE Gauss - - -!******************************************************************** -! determinant of a 3x3 matrix -!******************************************************************** - real(pReal) FUNCTION math_det3x3(m) - - use prec, only: pReal,pInt - implicit none - - real(pReal) m(3,3) - - math_det3x3 = m(1,1)*(m(2,2)*m(3,3)-m(2,3)*m(3,2)) & - -m(1,2)*(m(2,1)*m(3,3)-m(2,3)*m(3,1)) & - +m(1,3)*(m(2,1)*m(3,2)-m(2,2)*m(3,1)) - return - - END FUNCTION - - -!******************************************************************** -! convert symmetric 3x3 matrix into Mandel vector 6x1 -!******************************************************************** - FUNCTION math_Mandel33to6(m33) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3,3) :: m33 - real(pReal), dimension(6) :: math_Mandel33to6 - integer(pInt) i - - forall (i=1:6) math_Mandel33to6(i) = nrmMandel(i)*m33(mapMandel(1,i),mapMandel(2,i)) - return - - END FUNCTION - - -!******************************************************************** -! convert Mandel 6x1 back to symmetric 3x3 matrix -!******************************************************************** - FUNCTION math_Mandel6to33(v6) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(6) :: v6 - real(pReal), dimension(3,3) :: math_Mandel6to33 - integer(pInt) i - - forall (i=1:6) - math_Mandel6to33(mapMandel(1,i),mapMandel(2,i)) = invnrmMandel(i)*v6(i) - math_Mandel6to33(mapMandel(2,i),mapMandel(1,i)) = invnrmMandel(i)*v6(i) - end forall - return - - END FUNCTION - - -!******************************************************************** -! convert symmetric 3x3x3x3 tensor into Mandel matrix 6x6 -!******************************************************************** - FUNCTION math_Mandel3333to66(m3333) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3,3,3,3) :: m3333 - real(pReal), dimension(6,6) :: math_Mandel3333to66 - integer(pInt) i,j - - forall (i=1:6,j=1:6) math_Mandel3333to66(i,j) = & - nrmMandel(i)*nrmMandel(j)*m3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) - return - - END FUNCTION - - -!******************************************************************** -! convert Mandel matrix 6x6 back to symmetric 3x3x3x3 tensor -!******************************************************************** - FUNCTION math_Mandel66to3333(m66) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(6,6) :: m66 - real(pReal), dimension(3,3,3,3) :: math_Mandel66to3333 - integer(pInt) i,j - - forall (i=1:6,j=1:6) - math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) - math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(1,j),mapMandel(2,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) - math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(2,j),mapMandel(1,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) - math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(2,j),mapMandel(1,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) - end forall - return - - END FUNCTION - - -!******************************************************************** -! convert Voigt matrix 6x6 back to symmetric 3x3x3x3 tensor -!******************************************************************** - FUNCTION math_Voigt66to3333(m66) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(6,6) :: m66 - real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 - integer(pInt) i,j - - forall (i=1:6,j=1:6) - math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) - math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) - math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) - math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) - end forall - return - - END FUNCTION - - -!******************************************************************** -! Euler angles from orientation matrix -!******************************************************************** - FUNCTION math_RtoEuler(R) - - use prec, only: pReal, pInt - implicit none - - real(pReal) R(3,3), sqhkl, squvw, sqhk, val - real(pReal), dimension(3) :: math_RtoEuler - - sqhkl=sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)+R(3,3)*R(3,3)) - squvw=sqrt(R(1,1)*R(1,1)+R(2,1)*R(2,1)+R(3,1)*R(3,1)) - sqhk=sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)) -! calculate PHI - val=R(3,3)/sqhkl - - if(val > 1.0_pReal) val = 1.0_pReal - if(val < -1.0_pReal) val = -1.0_pReal - - math_RtoEuler(2) = acos(val) - - if(math_RtoEuler(2) < 1.0e-30_pReal) then -! calculate phi2 - math_RtoEuler(3) = 0.0_pReal -! calculate phi1 - val=R(1,1)/squvw - if(val > 1.0_pReal) val = 1.0_pReal - if(val < -1.0_pReal) val = -1.0_pReal - - math_RtoEuler(1) = acos(val) - if(R(2,1) > 0.0_pReal) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) - else -! calculate phi2 - val=R(2,3)/sqhk - if(val > 1.0_pReal) val = 1.0_pReal - if(val < -1.0_pReal) val = -1.0_pReal - - math_RtoEuler(3) = acos(val) - if(R(1,3) < 0.0) math_RtoEuler(3) = 2.0_pReal*pi-math_RtoEuler(3) -! calculate phi1 - val=-R(3,2)/sin(math_RtoEuler(2)) - if(val > 1.0_pReal) val = 1.0_pReal - if(val < -1.0_pReal) val = -1.0_pReal - - math_RtoEuler(1) = acos(val) - if(R(3,1) < 0.0) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) - end if - return - - END FUNCTION - - -!**************************************************************** -! rotation matrix from axis and angle (in radians) -!**************************************************************** - FUNCTION math_RodrigToR(axis,omega) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3) :: axis, axisNrm - real(pReal), dimension(3,3) :: math_RodrigToR - real(pReal) omega, s,c - integer(pInt) i - - forall (i=1:3) axisNrm(i) = axis(i)/sqrt(dot_product(axis,axis)) - s = sin(omega) - c = cos(omega) - math_RodrigtoR(1,1) = (1.0_pReal - axisNrm(1)**2)*c + axisNrm(1)**2 - math_RodrigtoR(1,2) = axisNrm(1)*axisNrm(2)*(1.0_pReal - c) + axisNrm(3)*s - math_RodrigtoR(1,3) = axisNrm(1)*axisNrm(3)*(1.0_pReal - c) - axisNrm(2)*s - math_RodrigtoR(2,1) = axisNrm(1)*axisNrm(2)*(1.0_pReal - c) - axisNrm(3)*s - math_RodrigtoR(2,2) = (1.0_pReal - axisNrm(2)**2)*c + axisNrm(2)**2 - math_RodrigtoR(2,3) = axisNrm(2)*axisNrm(3)*(1.0_pReal - c) + axisNrm(1)*s - math_RodrigtoR(3,1) = axisNrm(1)*axisNrm(3)*(1.0_pReal - c) + axisNrm(2)*s - math_RodrigtoR(3,2) = axisNrm(2)*axisNrm(3)*(1.0_pReal - c) - axisNrm(1)*s - math_RodrigtoR(3,3) = (1.0_pReal - axisNrm(3)**2)*c + axisNrm(3)**2 - return - - END FUNCTION - - -!**************************************************************** -! rotation matrix from Euler angles -!**************************************************************** - PURE FUNCTION math_EulerToR (Euler) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3), intent(in) :: Euler - real(pReal), dimension(3,3) :: math_EulerToR - real(pReal) c1, c, c2, s1, s, s2 - - C1=COS(Euler(1)) - C=COS(Euler(2)) - C2=COS(Euler(3)) - S1=SIN(Euler(1)) - S=SIN(Euler(2)) - S2=SIN(Euler(3)) - math_EulerToR(1,1)=C1*C2-S1*S2*C - math_EulerToR(1,2)=S1*C2+C1*S2*C - math_EulerToR(1,3)=S2*S - math_EulerToR(2,1)=-C1*S2-S1*C2*C - math_EulerToR(2,2)=-S1*S2+C1*C2*C - math_EulerToR(2,3)=C2*S - math_EulerToR(3,1)=S1*S - math_EulerToR(3,2)=-C1*S - math_EulerToR(3,3)=C - return - - END FUNCTION - - -!************************************************************************** -! disorientation angle between two sets of Euler angles -!************************************************************************** - function math_disorient(EulerA,EulerB) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3):: EulerA,EulerB - real(pReal), dimension(3,3) :: r - real(pReal) math_disorient, tr - - r = matmul(math_EulerToR(EulerB),transpose(math_EulerToR(EulerA))) - tr = (r(1,1)+r(2,2)+r(3,3)-1.0_pReal)*0.4999999_pReal - math_disorient = abs(0.5_pReal*pi-asin(tr)) - return - - END FUNCTION - - -!******************************************************************** -! draw a random sample from Euler space -!******************************************************************** - FUNCTION math_sampleRandomOri() - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3) :: math_sampleRandomOri, rnd - - call halton(3,rnd) - math_sampleRandomOri(1) = rnd(1)*2.0_pReal*pi - math_sampleRandomOri(2) = acos(2.0_pReal*rnd(2)-1.0_pReal) - math_sampleRandomOri(3) = rnd(3)*2.0_pReal*pi - return - - END FUNCTION - - -!******************************************************************** -! draw a random sample from Gauss component -! with noise (in radians) half-width -!******************************************************************** - FUNCTION math_sampleGaussOri(center,noise) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3) :: math_sampleGaussOri, center, disturb - real(pReal), dimension(3), parameter :: origin = (/0.0_pReal,0.0_pReal,0.0_pReal/) - real(pReal), dimension(5) :: rnd - real(pReal) noise,scatter,cosScatter - integer(pInt) i - -if (noise==0.0) then - math_sampleGaussOri = center - return -endif - -! Helming uses different distribution with Bessel functions -! therefore the gauss scatter width has to be scaled differently - scatter = 0.95_pReal * noise - cosScatter = cos(scatter) - - do - call halton(5,rnd) - forall (i=1:3) rnd(i) = 2.0_pReal*rnd(i)-1.0_pReal ! expand 1:3 to range [-1,+1] - disturb(1) = scatter * rnd(1) ! phi1 - disturb(2) = sign(1.0_pReal,rnd(2))*acos(cosScatter+(1.0_pReal-cosScatter)*rnd(4)) ! Phi - disturb(3) = scatter * rnd(2) ! phi2 - if (rnd(5) <= exp(-1.0_pReal*(math_disorient(origin,disturb)/scatter)**2)) exit - enddo - math_sampleGaussOri = math_RtoEuler(matmul(math_EulerToR(disturb),math_EulerToR(center))) - return - - END FUNCTION - - -!******************************************************************** -! draw a random sample from Fiber component -! with noise (in radians) -!******************************************************************** - FUNCTION math_sampleFiberOri(alpha,beta,noise) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3) :: math_sampleFiberOri, fiberInC,fiberInS,axis - real(pReal), dimension(2) :: alpha,beta, rnd - real(pReal), dimension(3,3) :: oRot,fRot,pRot - real(pReal) noise, scatter, cos2Scatter, angle - integer(pInt), dimension(2,3), parameter :: rotMap = reshape((/2,3, 3,1, 1,2/),(/2,3/)) - integer(pInt) i - -! Helming uses different distribution with Bessel functions -! therefore the gauss scatter width has to be scaled differently - scatter = 0.95_pReal * noise - cos2Scatter = cos(2.0_pReal*scatter) - -! fiber axis in crystal coordinate system - fiberInC(1)=sin(alpha(1))*cos(alpha(2)) - fiberInC(2)=sin(alpha(1))*sin(alpha(2)) - fiberInC(3)=cos(alpha(1)) -! fiber axis in sample coordinate system - fiberInS(1)=sin(beta(1))*cos(beta(2)) - fiberInS(2)=sin(beta(1))*sin(beta(2)) - fiberInS(3)=cos(beta(1)) - -! ---# rotation matrix from sample to crystal system #--- - angle=-acos(dot_product(fiberInC,fiberInS)) - if(angle /= 0.0_pReal) then -! rotation axis between sample and crystal system (cross product) - forall(i=1:3) axis(i) = fiberInC(rotMap(1,i))*fiberInS(rotMap(2,i))-fiberInC(rotMap(2,i))*fiberInS(rotMap(1,i)) - oRot = math_RodrigtoR(axis,angle) - else - oRot = math_I3 - end if - -! ---# rotation matrix about fiber axis (random angle) #--- - call halton(1,rnd) - fRot = math_RodrigToR(fiberInS,axis(3)*2.0_pReal*pi) - -! ---# rotation about random axis perpend to fiber #--- -! random axis pependicular to fiber axis - call halton(2,axis) - if (fiberInS(3) /= 0.0_pReal) then - axis(3)=-(axis(1)*fiberInS(1)+axis(2)*fiberInS(2))/fiberInS(3) - else if(fiberInS(2) /= 0.0_pReal) then - axis(3)=axis(2) - axis(2)=-(axis(1)*fiberInS(1)+axis(3)*fiberInS(3))/fiberInS(2) - else if(fiberInS(1) /= 0.0_pReal) then - axis(3)=axis(1) - axis(1)=-(axis(2)*fiberInS(2)+axis(3)*fiberInS(3))/fiberInS(1) - end if - -! scattered rotation angle - do - call halton(2,rnd) - angle = acos(cos2Scatter+(1.0_pReal-cos2Scatter)*rnd(1)) - if (rnd(2) <= exp(-1.0_pReal*(angle/scatter)**2)) exit - enddo - call halton(1,rnd) - if (rnd(1) <= 0.5) angle = -angle - pRot = math_RodrigtoR(axis,angle) - -! ---# apply the three rotations #--- - math_sampleFiberOri = math_RtoEuler(matmul(pRot,matmul(fRot,oRot))) - return - - END FUNCTION - - -!******************************************************************** -! symmetric Euler angles for given symmetry string -! 'triclinic' or '', 'monoclinic', 'orthotropic' -!******************************************************************** - PURE FUNCTION math_symmetricEulers(sym,Euler) - - use prec, only: pReal, pInt - implicit none - - character(len=80), intent(in) :: sym - real(pReal), dimension(3), intent(in) :: Euler - real(pReal), dimension(3,3) :: math_symmetricEulers - integer(pInt) i,j - - math_symmetricEulers(1,1) = pi+Euler(1) - math_symmetricEulers(2,1) = Euler(2) - math_symmetricEulers(3,1) = Euler(3) - - math_symmetricEulers(1,2) = pi-Euler(1) - math_symmetricEulers(2,2) = pi-Euler(2) - math_symmetricEulers(3,2) = pi+Euler(3) - - math_symmetricEulers(1,3) = 2.0_pReal*pi-Euler(1) - math_symmetricEulers(2,3) = pi-Euler(2) - math_symmetricEulers(3,3) = pi+Euler(3) - - forall (i=1:3,j=1:3) math_symmetricEulers(j,i) = modulo(math_symmetricEulers(j,i),2.0_pReal*pi) - - select case (sym) - case ('orthotropic') ! all done - - case ('monoclinic') ! return only first - math_symmetricEulers(:,2:3) = 0.0_pReal - case default ! return blank - math_symmetricEulers = 0.0_pReal - end select - return - - END FUNCTION - - -!**************************************************************** - subroutine math_pDecomposition(FE,U,R,error) -!-----FE=RU -!**************************************************************** - use prec, only: pReal, pInt - implicit none - - logical error - real(pReal) FE(3,3),R(3,3),U(3,3),CE(3,3),EW1,EW2,EW3,EB1(3,3),EB2(3,3),EB3(3,3),UI(3,3),det - - error = .false. - ce=matmul(transpose(fe),fe) - CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3) - U=DSQRT(EW1)*EB1+DSQRT(EW2)*EB2+DSQRT(EW3)*EB3 - call math_invert3x3(U,UI,det,error) - if (.not. error) R = matmul(fe,ui) - return - - END SUBROUTINE - - -!********************************************************************** - subroutine math_spectral1(M,EW1,EW2,EW3,EB1,EB2,EB3) -!**** EIGENWERTE UND EIGENWERTBASIS DER SYMMETRISCHEN 3X3 MATRIX M - - use prec, only: pReal, pInt - implicit none - - real(pReal) M(3,3),EB1(3,3),EB2(3,3),EB3(3,3),EW1,EW2,EW3 - real(pReal) HI1M,HI2M,HI3M,TOL,R,S,T,P,Q,RHO,PHI,Y1,Y2,Y3,D1,D2,D3 - real(pReal) C1,C2,C3,M1(3,3),M2(3,3),M3(3,3),arg - TOL=1.e-14_pReal - CALL math_hi(M,HI1M,HI2M,HI3M) - R=-HI1M - S= HI2M - T=-HI3M - P=S-R**2.0_pReal/3.0_pReal - Q=2.0_pReal/27.0_pReal*R**3.0_pReal-R*S/3.0_pReal+T - EB1=0.0_pReal - EB2=0.0_pReal - EB3=0.0_pReal - IF((ABS(P).LT.TOL).AND.(ABS(Q).LT.TOL))THEN -! DREI GLEICHE EIGENWERTE - EW1=HI1M/3.0_pReal - EW2=EW1 - EW3=EW1 -! this is not really correct, but this way U is calculated -! correctly in PDECOMPOSITION (correct is EB?=I) - EB1(1,1)=1.0_pReal - EB2(2,2)=1.0_pReal - EB3(3,3)=1.0_pReal - ELSE - RHO=SQRT(-3.0_pReal*P**3.0_pReal)/9.0_pReal - arg=-Q/RHO/2.0_pReal - if(arg.GT.1) arg=1 - if(arg.LT.-1) arg=-1 - PHI=ACOS(arg) - Y1=2*RHO**(1.0_pReal/3.0_pReal)*COS(PHI/3.0_pReal) - Y2=2*RHO**(1.0_pReal/3.0_pReal)*COS(PHI/3.0_pReal+2.0_pReal/3.0_pReal*PI) - Y3=2*RHO**(1.0_pReal/3.0_pReal)*COS(PHI/3.0_pReal+4.0_pReal/3.0_pReal*PI) - EW1=Y1-R/3.0_pReal - EW2=Y2-R/3.0_pReal - EW3=Y3-R/3.0_pReal - C1=ABS(EW1-EW2) - C2=ABS(EW2-EW3) - C3=ABS(EW3-EW1) - - IF(C1.LT.TOL) THEN -! EW1 is equal to EW2 - D3=1.0_pReal/(EW3-EW1)/(EW3-EW2) - M1=M-EW1*math_I3 - M2=M-EW2*math_I3 - EB3=MATMUL(M1,M2)*D3 - EB1=math_I3-EB3 -! both EB2 and EW2 are set to zero so that they do not -! contribute to U in PDECOMPOSITION - EW2=0.0_pReal - ELSE IF(C2.LT.TOL) THEN -! EW2 is equal to EW3 - D1=1.0_pReal/(EW1-EW2)/(EW1-EW3) - M2=M-math_I3*EW2 - M3=M-math_I3*EW3 - EB1=MATMUL(M2,M3)*D1 - EB2=math_I3-EB1 -! both EB3 and EW3 are set to zero so that they do not -! contribute to U in PDECOMPOSITION - EW3=0.0_pReal - ELSE IF(C3.LT.TOL) THEN -! EW1 is equal to EW3 - D2=1.0_pReal/(EW2-EW1)/(EW2-EW3) - M1=M-math_I3*EW1 - M3=M-math_I3*EW3 - EB2=MATMUL(M1,M3)*D2 - EB1=math_I3-EB2 -! both EB3 and EW3 are set to zero so that they do not -! contribute to U in PDECOMPOSITION - EW3=0.0_pReal - ELSE -! all three eigenvectors are different - D1=1.0_pReal/(EW1-EW2)/(EW1-EW3) - D2=1.0_pReal/(EW2-EW1)/(EW2-EW3) - D3=1.0_pReal/(EW3-EW1)/(EW3-EW2) - M1=M-EW1*math_I3 - M2=M-EW2*math_I3 - M3=M-EW3*math_I3 - EB1=MATMUL(M2,M3)*D1 - EB2=MATMUL(M1,M3)*D2 - EB3=MATMUL(M1,M2)*D3 - END IF - END IF - RETURN - END SUBROUTINE - - -!********************************************************************** -!**** HAUPTINVARIANTEN HI1M, HI2M, HI3M DER 3X3 MATRIX M - - SUBROUTINE math_hi(M,HI1M,HI2M,HI3M) - use prec, only: pReal, pInt - implicit none - - real(pReal) M(3,3),HI1M,HI2M,HI3M - - HI1M=M(1,1)+M(2,2)+M(3,3) - HI2M=HI1M**2/2.0_pReal-(M(1,1)**2+M(2,2)**2+M(3,3)**2)/2.0_pReal-M(1,2)*M(2,1)-M(1,3)*M(3,1)-M(2,3)*M(3,2) - HI3M=math_det3x3(M) -! QUESTION: is 3rd equiv det(M) ?? if yes, use function math_det !agreed on YES - return - - END SUBROUTINE - - - SUBROUTINE get_seed(seed) -! -!******************************************************************************* -! -!! GET_SEED returns a seed for the random number generator. -! -! -! Discussion: -! -! The seed depends on the current time, and ought to be (slightly) -! different every millisecond. Once the seed is obtained, a random -! number generator should be called a few times to further process -! the seed. -! -! Modified: -! -! 27 June 2000 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Output, integer SEED, a pseudorandom seed value. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, only: pReal, pInt - implicit none - - integer(pInt) seed - real(pReal) temp - character ( len = 10 ) time - character ( len = 8 ) today - integer(pInt) values(8) - character ( len = 5 ) zone - - call date_and_time ( today, time, zone, values ) - - temp = 0.0D+00 - - temp = temp + dble ( values(2) - 1 ) / 11.0D+00 - temp = temp + dble ( values(3) - 1 ) / 30.0D+00 - temp = temp + dble ( values(5) ) / 23.0D+00 - temp = temp + dble ( values(6) ) / 59.0D+00 - temp = temp + dble ( values(7) ) / 59.0D+00 - temp = temp + dble ( values(8) ) / 999.0D+00 - temp = temp / 6.0D+00 - - if ( temp <= 0.0D+00 ) then - temp = 1.0D+00 / 3.0D+00 - else if ( 1.0D+00 <= temp ) then - temp = 2.0D+00 / 3.0D+00 - end if - - seed = int ( dble ( huge ( 1 ) ) * temp , pInt) -! -! Never use a seed of 0 or maximum integer. -! - if ( seed == 0 ) then - seed = 1 - end if - - if ( seed == huge ( 1 ) ) then - seed = seed - 1 - end if - - return - - END SUBROUTINE - - - subroutine halton ( ndim, r ) -! -!******************************************************************************* -! -!! HALTON computes the next element in the Halton sequence. -! -! -! Modified: -! -! 09 March 2003 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Input, integer NDIM, the dimension of the element. -! -! Output, real R(NDIM), the next element of the current Halton -! sequence. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, ONLY: pReal, pInt - implicit none - - integer(pInt) ndim - - integer(pInt) base(ndim) - real(pReal) r(ndim) - integer(pInt) seed - integer(pInt) value(1) - - call halton_memory ( 'GET', 'SEED', 1, value ) - seed = value(1) - - call halton_memory ( 'GET', 'BASE', ndim, base ) - - call i_to_halton ( seed, base, ndim, r ) - - value(1) = 1 - call halton_memory ( 'INC', 'SEED', 1, value ) - - return - - END SUBROUTINE - - - subroutine halton_memory ( action, name, ndim, value ) -! -!******************************************************************************* -! -!! HALTON_MEMORY sets or returns quantities associated with the Halton sequence. -! -! -! Modified: -! -! 09 March 2003 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Input, character ( len = * ) ACTION, the desired action. -! 'GET' means get the value of a particular quantity. -! 'SET' means set the value of a particular quantity. -! 'INC' means increment the value of a particular quantity. -! (Only the SEED can be incremented.) -! -! Input, character ( len = * ) NAME, the name of the quantity. -! 'BASE' means the Halton base or bases. -! 'NDIM' means the spatial dimension. -! 'SEED' means the current Halton seed. -! -! Input/output, integer NDIM, the dimension of the quantity. -! If ACTION is 'SET' and NAME is 'BASE', then NDIM is input, and -! is the number of entries in VALUE to be put into BASE. -! -! Input/output, integer VALUE(NDIM), contains a value. -! If ACTION is 'SET', then on input, VALUE contains values to be assigned -! to the internal variable. -! If ACTION is 'GET', then on output, VALUE contains the values of -! the specified internal variable. -! If ACTION is 'INC', then on input, VALUE contains the increment to -! be added to the specified internal variable. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, only: pReal, pInt - implicit none - - character ( len = * ) action - integer(pInt), allocatable, save :: base(:) - logical, save :: first_call = .true. - integer(pInt) i - character ( len = * ) name - integer(pInt) ndim - integer(pInt), save :: ndim_save = 0 - integer(pInt), save :: seed = 1 - integer(pInt) value(*) - - if ( first_call ) then - ndim_save = 1 - allocate ( base(ndim_save) ) - base(1) = 2 - first_call = .false. - end if -! -! Set -! - if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then - - if ( name(1:1) == 'B' .or. name(1:1) == 'b' ) then - - if ( ndim_save /= ndim ) then - deallocate ( base ) - ndim_save = ndim - allocate ( base(ndim_save) ) - end if - - base(1:ndim) = value(1:ndim) - - else if ( name(1:1) == 'N' .or. name(1:1) == 'n' ) then - - if ( ndim_save /= value(1) ) then - deallocate ( base ) - ndim_save = value(1) - allocate ( base(ndim_save) ) - do i = 1, ndim_save - base(i) = prime ( i ) - end do - else - ndim_save = value(1) - end if - - else if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then - - seed = value(1) - - end if -! -! Get -! - else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then - - if ( name(1:1) == 'B' .or. name(1:1) == 'b' ) then - - if ( ndim /= ndim_save ) then - deallocate ( base ) - ndim_save = ndim - allocate ( base(ndim_save) ) - do i = 1, ndim_save - base(i) = prime(i) - end do - end if - - value(1:ndim_save) = base(1:ndim_save) - - else if ( name(1:1) == 'N' .or. name(1:1) == 'n' ) then - - value(1) = ndim_save - - else if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then - - value(1) = seed - - end if -! -! Increment -! - else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then - - if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then - seed = seed + value(1) - end if - - end if - - return - - END SUBROUTINE - - - - subroutine halton_ndim_set ( ndim ) -! -!******************************************************************************* -! -!! HALTON_NDIM_SET sets the dimension for a Halton sequence. -! -! -! Modified: -! -! 26 February 2001 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Input, integer NDIM, the dimension of the Halton vectors. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, only: pReal, pInt - implicit none - - integer(pInt) ndim - integer(pInt) value(1) - - value(1) = ndim - call halton_memory ( 'SET', 'NDIM', 1, value ) - - return - - END SUBROUTINE - - - subroutine halton_seed_set ( seed ) -! -!******************************************************************************* -! -!! HALTON_SEED_SET sets the "seed" for the Halton sequence. -! -! -! Discussion: -! -! Calling HALTON repeatedly returns the elements of the -! Halton sequence in order, starting with element number 1. -! An internal counter, called SEED, keeps track of the next element -! to return. Each time the routine is called, the SEED-th element -! is computed, and then SEED is incremented by 1. -! -! To restart the Halton sequence, it is only necessary to reset -! SEED to 1. It might also be desirable to reset SEED to some other value. -! This routine allows the user to specify any value of SEED. -! -! The default value of SEED is 1, which restarts the Halton sequence. -! -! Modified: -! -! 26 February 2001 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Input, integer SEED, the seed for the Halton sequence. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, only: pReal, pInt - implicit none - - integer(pInt), parameter :: ndim = 1 - - integer(pInt) seed - integer(pInt) value(ndim) - - value(1) = seed - call halton_memory ( 'SET', 'SEED', ndim, value ) - - return - - END SUBROUTINE - - - subroutine i_to_halton ( seed, base, ndim, r ) -! -!******************************************************************************* -! -!! I_TO_HALTON computes an element of a Halton sequence. -! -! -! Reference: -! -! J H Halton, -! On the efficiency of certain quasi-random sequences of points -! in evaluating multi-dimensional integrals, -! Numerische Mathematik, -! Volume 2, pages 84-90, 1960. -! -! Modified: -! -! 26 February 2001 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Input, integer SEED, the index of the desired element. -! Only the absolute value of SEED is considered. SEED = 0 is allowed, -! and returns R = 0. -! -! Input, integer BASE(NDIM), the Halton bases, which should be -! distinct prime numbers. This routine only checks that each base -! is greater than 1. -! -! Input, integer NDIM, the dimension of the sequence. -! -! Output, real R(NDIM), the SEED-th element of the Halton sequence -! for the given bases. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, ONLY: pReal, pInt - implicit none - - integer(pInt) ndim - - integer(pInt) base(ndim) - real(pReal) base_inv(ndim) - integer(pInt) digit(ndim) - integer(pInt) i - real(pReal) r(ndim) - integer(pInt) seed - integer(pInt) seed2(ndim) - - seed2(1:ndim) = abs ( seed ) - - r(1:ndim) = 0.0_pReal - - if ( any ( base(1:ndim) <= 1 ) ) then - write ( *, '(a)' ) ' ' - write ( *, '(a)' ) 'I_TO_HALTON - Fatal error!' - write ( *, '(a)' ) ' An input base BASE is <= 1!' - do i = 1, ndim - write ( *, '(i6,i6)' ) i, base(i) - end do - call flush(6) - stop - end if - - base_inv(1:ndim) = 1.0_pReal / real ( base(1:ndim), pReal ) - - do while ( any ( seed2(1:ndim) /= 0 ) ) - digit(1:ndim) = mod ( seed2(1:ndim), base(1:ndim) ) - r(1:ndim) = r(1:ndim) + real ( digit(1:ndim), pReal ) * base_inv(1:ndim) - base_inv(1:ndim) = base_inv(1:ndim) / real ( base(1:ndim), pReal ) - seed2(1:ndim) = seed2(1:ndim) / base(1:ndim) - end do - - return - - END SUBROUTINE - - - function prime ( n ) -! -!******************************************************************************* -! -!! PRIME returns any of the first PRIME_MAX prime numbers. -! -! -! Note: -! -! PRIME_MAX is 1500, and the largest prime stored is 12553. -! -! Modified: -! -! 21 June 2002 -! -! Author: -! -! John Burkardt -! -! Reference: -! -! Milton Abramowitz and Irene Stegun, -! Handbook of Mathematical Functions, -! US Department of Commerce, 1964, pages 870-873. -! -! Daniel Zwillinger, -! CRC Standard Mathematical Tables and Formulae, -! 30th Edition, -! CRC Press, 1996, pages 95-98. -! -! Parameters: -! -! Input, integer N, the index of the desired prime number. -! N = -1 returns PRIME_MAX, the index of the largest prime available. -! N = 0 is legal, returning PRIME = 1. -! It should generally be true that 0 <= N <= PRIME_MAX. -! -! Output, integer PRIME, the N-th prime. If N is out of range, PRIME -! is returned as 0. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, only: pReal, pInt - implicit none - - integer(pInt), parameter :: prime_max = 1500 - - integer(pInt), save :: icall = 0 - integer(pInt) n - integer(pInt), save, dimension ( prime_max ) :: npvec - integer(pInt) prime - - if ( icall == 0 ) then - - icall = 1 - - npvec(1:100) = (/& - 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, & - 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, & - 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, & - 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, & - 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, & - 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, & - 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, & - 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, & - 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, & - 467, 479, 487, 491, 499, 503, 509, 521, 523, 541 /) - - npvec(101:200) = (/ & - 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, & - 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, & - 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, & - 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, & - 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, & - 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, & - 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, & - 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, & - 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, & - 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223 /) - - npvec(201:300) = (/ & - 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, & - 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, & - 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, & - 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, & - 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, & - 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, & - 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, & - 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, & - 91823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, & - 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987 /) - - npvec(301:400) = (/ & - 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, & - 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, & - 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, & - 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, & - 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, & - 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, & - 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, & - 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, & - 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, & - 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741 /) - - npvec(401:500) = (/ & - 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, & - 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, & - 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, & - 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, & - 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, & - 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, & - 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, & - 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, & - 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, & - 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571 /) - - npvec(501:600) = (/ & - 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, & - 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, & - 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, & - 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, & - 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, & - 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, & - 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, & - 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, & - 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, & - 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409 /) - - npvec(601:700) = (/ & - 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, & - 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, & - 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, & - 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, & - 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, & - 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, & - 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, & - 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, & - 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, & - 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279 /) - - npvec(701:800) = (/ & - 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, & - 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, & - 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, & - 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, & - 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, & - 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, & - 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, & - 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, & - 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, & - 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133 /) - - npvec(801:900) = (/ & - 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, & - 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, & - 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, & - 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, & - 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, & - 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, & - 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, & - 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, & - 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, & - 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997 /) - - npvec(901:1000) = (/ & - 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, & - 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, & - 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, & - 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, & - 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, & - 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, & - 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, & - 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, & - 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, & - 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919 /) - - npvec(1001:1100) = (/ & - 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, & - 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, & - 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, & - 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, & - 8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387, & - 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, & - 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, & - 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, & - 8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, & - 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831 /) - - npvec(1101:1200) = (/ & - 8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, & - 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011, & - 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, & - 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, & - 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, & - 9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377, & - 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, & - 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533, & - 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, & - 9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733 /) - - npvec(1201:1300) = (/ & - 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, & - 9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887, & - 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973,10007, & - 10009,10037,10039,10061,10067,10069,10079,10091,10093,10099, & - 10103,10111,10133,10139,10141,10151,10159,10163,10169,10177, & - 10181,10193,10211,10223,10243,10247,10253,10259,10267,10271, & - 10273,10289,10301,10303,10313,10321,10331,10333,10337,10343, & - 10357,10369,10391,10399,10427,10429,10433,10453,10457,10459, & - 10463,10477,10487,10499,10501,10513,10529,10531,10559,10567, & - 10589,10597,10601,10607,10613,10627,10631,10639,10651,10657 /) - - npvec(1301:1400) = (/ & - 10663,10667,10687,10691,10709,10711,10723,10729,10733,10739, & - 10753,10771,10781,10789,10799,10831,10837,10847,10853,10859, & - 10861,10867,10883,10889,10891,10903,10909,19037,10939,10949, & - 10957,10973,10979,10987,10993,11003,11027,11047,11057,11059, & - 11069,11071,11083,11087,11093,11113,11117,11119,11131,11149, & - 11159,11161,11171,11173,11177,11197,11213,11239,11243,11251, & - 11257,11261,11273,11279,11287,11299,11311,11317,11321,11329, & - 11351,11353,11369,11383,11393,11399,11411,11423,11437,11443, & - 11447,11467,11471,11483,11489,11491,11497,11503,11519,11527, & - 11549,11551,11579,11587,11593,11597,11617,11621,11633,11657 /) - - npvec(1401:1500) = (/ & - 11677,11681,11689,11699,11701,11717,11719,11731,11743,11777, & - 11779,11783,11789,11801,11807,11813,11821,11827,11831,11833, & - 11839,11863,11867,11887,11897,11903,11909,11923,11927,11933, & - 11939,11941,11953,11959,11969,11971,11981,11987,12007,12011, & - 12037,12041,12043,12049,12071,12073,12097,12101,12107,12109, & - 12113,12119,12143,12149,12157,12161,12163,12197,12203,12211, & - 12227,12239,12241,12251,12253,12263,12269,12277,12281,12289, & - 12301,12323,12329,12343,12347,12373,12377,12379,12391,12401, & - 12409,12413,12421,12433,12437,12451,12457,12473,12479,12487, & - 12491,12497,12503,12511,12517,12527,12539,12541,12547,12553 /) - - end if - - if ( n == -1 ) then - prime = prime_max - else if ( n == 0 ) then - prime = 1 - else if ( n <= prime_max ) then - prime = npvec(n) - else - prime = 0 - write ( 6, '(a)' ) ' ' - write ( 6, '(a)' ) 'PRIME - Fatal error!' - write ( 6, '(a,i6)' ) ' Illegal prime index N = ', n - write ( 6, '(a,i6)' ) ' N must be between 0 and PRIME_MAX =',prime_max - call flush(6) - stop - end if - - return - - END FUNCTION - - - END MODULE math - diff --git a/dislocation_based_subroutine/mattex.mpie b/dislocation_based_subroutine/mattex.mpie deleted file mode 100644 index d7ecf2d52..000000000 --- a/dislocation_based_subroutine/mattex.mpie +++ /dev/null @@ -1,24 +0,0 @@ - -[Aluminium] -temperature 293 -crystal_structure 1 -Nslip 12 -C11 106.75e3 -C12 60.41e3 -C44 28.34e3 -rho0 1.5e11 -bg 2.86e-10 -Qedge 3.0e-19 -tau0 0.0 -c1 0.1 -c2 2.0 -c3 0.4 -c4 0.05 -c5 10.0 -interaction_coefficients 1.0 2.2 3.0 1.6 3.8 4.5 - - -[cube SX] -symmetry monoclinic -Ngrains 1 -(gauss) phi1 0.0 phi 0.0 phi2 0.0 scatter 0.0 fraction 1.0 diff --git a/dislocation_based_subroutine/mesh.f90 b/dislocation_based_subroutine/mesh.f90 deleted file mode 100644 index 1865c647e..000000000 --- a/dislocation_based_subroutine/mesh.f90 +++ /dev/null @@ -1,702 +0,0 @@ - -!############################################################## - MODULE mesh -!############################################################## - - use prec, only: pReal,pInt - implicit none - -! --------------------------- -! _Nelems : total number of elements in mesh -! _NcpElems : total number of CP elements in mesh -! _Nnodes : total number of nodes in mesh -! _maxNnodes : max number of nodes in any CP element -! _maxNips : max number of IPs in any CP element -! _maxNipNeighbors : max number of IP neighbors in any CP element -! _maxNsharedElems : max number of CP elements sharing a node -! -! _element : FEid, type, material, texture, node indices -! _node : x,y,z coordinates (initially!) -! _sharedElem : entryCount and list of elements containing node -! -! _mapFEtoCPelem : [sorted FEid, corresponding CPid] -! _mapFEtoCPnode : [sorted FEid, corresponding CPid] -! -! MISSING: these definitions should actually reside in the -! FE-solver specific part (different for MARC/ABAQUS)..! -! Hence, I suggest to prefix with "FE_" -! -! _mapElementtype : map MARC/ABAQUS elemtype to 1-maxN -! -! _Nnodes : # nodes in a specific type of element -! _Nips : # IPs in a specific type of element -! _NipNeighbors : # IP neighbors in a specific type of element -! _ipNeighbor : +x,-x,+y,-y,+z,-z list of intra-element IPs and -! (negative) neighbor faces per own IP in a specific type of element -! _NfaceNodes : # nodes per face in a specific type of element -! _nodeOnFace : list of node indices on each face of a specific type of element -! _ipAtNode : map node index to IP index in a specific type of element -! _nodeAtIP : map IP index to node index in a specific type of element -! _ipNeighborhood : 6 or less neighboring IPs as [element_num, IP_index] -! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype -! --------------------------- - integer(pInt) mesh_Nelems,mesh_NcpElems - integer(pInt) mesh_Nnodes,mesh_maxNnodes,mesh_maxNips,mesh_maxNipNeighbors,mesh_maxNsharedElems - integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem,mesh_mapFEtoCPnode - integer(pInt), dimension(:,:), allocatable :: mesh_element, mesh_sharedElem - integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood - real(pReal), allocatable :: mesh_node (:,:) - - integer(pInt) :: hypoelasticTableStyle = 0 - integer(pInt), parameter :: FE_Nelemtypes = 2 - integer(pInt), parameter :: FE_maxNnodes = 8 - integer(pInt), parameter :: FE_maxNips = 8 - integer(pInt), parameter :: FE_maxNneighbors = 6 - integer(pInt), parameter :: FE_maxNfaceNodes = 4 - integer(pInt), parameter :: FE_maxNfaces = 6 - integer(pInt), dimension(200):: FE_mapElemtype - integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nnodes = & - (/8, & ! element 7 - 4 & ! element 134 - /) - integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nips = & - (/8, & ! element 7 - 1 & ! element 134 - /) - integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NipNeighbors = & - (/6, & ! element 7 - 4 & ! element 134 - /) - integer(pInt), dimension(FE_maxNfaces,FE_Nelemtypes), parameter :: FE_NfaceNodes = & - reshape((/& - 4,4,4,4,4,4, & ! element 7 - 3,3,3,3,0,0 & ! element 134 - /),(/FE_maxNfaces,FE_Nelemtypes/)) - integer(pInt), dimension(FE_maxNips,FE_Nelemtypes), parameter :: FE_nodeAtIP = & - reshape((/& - 1,2,4,3,5,6,8,7, & ! element 7 - 1,0,0,0,0,0,0,0 & ! element 134 - /),(/FE_maxNips,FE_Nelemtypes/)) - integer(pInt), dimension(FE_maxNnodes,FE_Nelemtypes), parameter :: FE_ipAtNode = & - reshape((/& - 1,2,4,3,5,6,8,7, & ! element 7 - 1,1,1,1,0,0,0,0 & ! element 134 - /),(/FE_maxNnodes,FE_Nelemtypes/)) - integer(pInt), dimension(FE_maxNfaceNodes,FE_maxNfaces,FE_Nelemtypes), parameter :: FE_nodeOnFace = & - reshape((/& - 1,2,3,4 , & ! element 7 - 2,1,5,6 , & - 3,2,6,7 , & - 3,4,8,7 , & - 4,1,5,8 , & - 8,7,6,5 , & - 1,2,3,0 , & ! element 134 - 1,4,2,0 , & - 2,3,4,0 , & - 1,3,4,0 , & - 0,0,0,0 , & - 0,0,0,0 & - /),(/FE_maxNfaceNodes,FE_maxNfaces,FE_Nelemtypes/)) - integer(pInt), dimension(FE_maxNneighbors,FE_maxNips,FE_Nelemtypes), parameter :: FE_ipNeighbor = & - reshape((/& - 2,-5, 3,-2, 5,-1 , & ! element 7 - -3, 1, 4,-2, 6,-1 , & - 4,-5,-4, 1, 7,-1 , & - -3, 3,-4, 2, 8,-1 , & - 6,-5, 7,-2,-6, 1 , & - -3, 5, 8,-2,-6, 2 , & - 8,-5,-4, 5,-6, 3 , & - -3, 7,-4, 6,-6, 4 , & - -1,-2,-3,-4, 0, 0 , & ! element 134 - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 & - /),(/FE_maxNneighbors,FE_maxNips,FE_Nelemtypes/)) - - CONTAINS -! --------------------------- -! subroutine mesh_init() -! function mesh_FEtoCPelement(FEid) -! function mesh_build_ipNeighorhood() -! subroutine mesh_parse_inputFile() -! --------------------------- - - -!*********************************************************** -! initialization -!*********************************************************** - SUBROUTINE mesh_init () - - use prec, only: pInt - use IO, only: IO_error,IO_open_InputFile - implicit none - - integer(pInt), parameter :: fileUnit = 222 - - mesh_Nelems = 0_pInt - mesh_NcpElems = 0_pInt - mesh_Nnodes = 0_pInt - mesh_maxNips = 0_pInt - mesh_maxNnodes = 0_pInt - mesh_maxNsharedElems = 0_pInt - - FE_mapElemtype = 1 ! MISSING this should be zero... - FE_mapElemtype( 7) = 1 - FE_mapElemtype(134) = 2 - -! call to various subrountes to parse the stuff from the input file... - if (IO_open_inputFile(fileUnit)) then - call mesh_get_meshDimensions(fileUnit) - call mesh_build_nodeMapping(fileUnit) - call mesh_build_elemMapping(fileUnit) - call mesh_get_nodeElemDimensions(fileUnit) - call mesh_build_nodes(fileUnit) - call mesh_build_elements(fileUnit) - call mesh_build_sharedElems(fileUnit) - call mesh_build_ipNeighborhood() - close (fileUnit) - else - call IO_error(100) - endif - - END SUBROUTINE - - -!*********************************************************** -! FE to CP id mapping by binary search thru lookup array -! -! valid questions are 'elem', 'node' -!*********************************************************** - FUNCTION mesh_FEasCP(what,id) - - use prec, only: pInt - use IO, only: IO_lc - implicit none - - character(len=*), intent(in) :: what - integer(pInt), intent(in) :: id - integer(pInt), dimension(:,:), pointer :: lookupMap - integer(pInt) mesh_FEasCP, lower,upper,center - - mesh_FEasCP = 0_pInt - select case(IO_lc(what(1:4))) - case('elem') - lookupMap => mesh_mapFEtoCPelem - case('node') - lookupMap => mesh_mapFEtoCPnode - case default - return - end select - - lower = 1_pInt - upper = size(lookupMap,2) - - ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds? - if (lookupMap(1,lower) == id) then - mesh_FEasCP = lookupMap(2,lower) - return - elseif (lookupMap(1,upper) == id) then - mesh_FEasCP = lookupMap(2,upper) - return - endif - - ! binary search in between bounds - do while (upper-lower > 1) - center = (lower+upper)/2 - if (lookupMap(1,center) < id) then - lower = center - elseif (lookupMap(1,center) > id) then - upper = center - else - mesh_FEasCP = lookupMap(2,center) - exit - end if - end do - return - - END FUNCTION - - -!*********************************************************** -! find face-matching element of same type -!!*********************************************************** - FUNCTION mesh_faceMatch(face,elem) - - use prec, only: pInt - implicit none - - integer(pInt) face,elem - integer(pInt) mesh_faceMatch - integer(pInt), dimension(FE_NfaceNodes(face,FE_mapElemtype(mesh_element(2,elem)))) :: nodeMap - integer(pInt) minN,NsharedElems,lonelyNode,faceNode,i,n,t - - minN = mesh_maxNsharedElems+1 ! init to worst case - mesh_faceMatch = 0_pInt ! intialize to "no match found" - t = FE_mapElemtype(mesh_element(2,elem)) ! figure elemType - - do faceNode=1,FE_NfaceNodes(face,t) ! loop over nodes on face - nodeMap(faceNode) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(faceNode,face,t),elem)) ! CP id of face node - NsharedElems = mesh_sharedElem(1,nodeMap(faceNode)) ! figure # shared elements for this node - if (NsharedElems < minN) then - minN = NsharedElems ! remember min # shared elems - lonelyNode = faceNode ! remember most lonely node - endif - end do -candidate: do i=1,minN ! iterate over lonelyNode's shared elements - mesh_faceMatch = mesh_sharedElem(1+i,nodeMap(lonelyNode)) ! present candidate elem - if (mesh_faceMatch == elem) then ! my own element ? - mesh_faceMatch = 0_pInt ! disregard - cycle candidate - endif - do faceNode=1,FE_NfaceNodes(face,t) ! check remaining face nodes to match - if (faceNode == lonelyNode) cycle ! disregard lonely node (matches anyway) - n = nodeMap(faceNode) - if (all(mesh_sharedElem(2:1+mesh_sharedElem(1,n),n) /= mesh_faceMatch)) then ! no ref to candidate elem? - mesh_faceMatch = 0_pInt ! set to "no match" (so far) - cycle candidate ! next candidate elem - endif - end do - exit ! surviving candidate - end do candidate - - return - - END FUNCTION - - -!******************************************************************** -! get count of elements, nodes, and cp elements in mesh -! for subsequent array allocations -! -! assign globals: -! _Nelems, _Nnodes, _NcpElems -!******************************************************************** - SUBROUTINE mesh_get_meshDimensions (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt) unit,i,pos(41) - character*300 line - -610 FORMAT(A300) - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,20) - - select case ( IO_lc(IO_Stringvalue(line,pos,1))) - case('table') - if (pos(1) == 6) hypoelasticTableStyle = IO_IntValue (line,pos,5) ! only recognize header entry for "table" - case('sizing') - mesh_Nelems = IO_IntValue (line,pos,3) - mesh_Nnodes = IO_IntValue (line,pos,4) - case('hypoelastic') - do i=1,4+hypoelasticTableStyle - read (unit,610,END=620) line - end do - pos = IO_stringPos(line,20) - if( IO_lc(IO_Stringvalue(line,pos,2)).eq.'to' )then - mesh_NcpElems = IO_IntValue(line,pos,3)-IO_IntValue(line,pos,1)+1 - else - mesh_NcpElems = mesh_NcpElems + pos(1) - do while( IO_lc(IO_Stringvalue(line,pos,pos(1))).eq.'c' ) - mesh_NcpElems = mesh_NcpElems - 1 ! Counted the c character from the line - read (unit,610,END=620) line - pos = IO_stringPos(line,20) - mesh_NcpElems = mesh_NcpElems + pos(1) - end do - end if - - end select - - end do - -620 return - - END SUBROUTINE - - -!******************************************************************** -! get maximum count of nodes, IPs, IP neighbors, and shared elements -! for subsequent array allocations -! -! assign globals: -! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems -!******************************************************************** - SUBROUTINE mesh_get_nodeElemDimensions (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), dimension (mesh_Nnodes) :: node_count - integer(pInt) unit,i,j,n,t,e - integer(pInt), dimension (133) :: pos - character*300 line - -610 FORMAT(A300) - - node_count = 0_pInt - - rewind(unit) - do - read (unit,610,END=630) line - pos = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then - read (unit,610,END=630) line ! Garbage line - do i=1,mesh_Nelems ! read all elements - read (unit,610,END=630) line - pos = IO_stringPos(line,66) ! limit to 64 nodes max (plus ID, type) - e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) - if (e /= 0) then - t = FE_mapElemtype(IO_intValue(line,pos,2)) - mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) - mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) - mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) - do j=1,FE_Nnodes(t) - n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2)) - node_count(n) = node_count(n)+1 - end do - end if - end do - exit - end if - end do - -630 mesh_maxNsharedElems = maxval(node_count) - - return - END SUBROUTINE - - -!******************************************************************** -! Build node mapping from FEM to CP -! -! allocate globals: -! _mapFEtoCPnode -!******************************************************************** - SUBROUTINE mesh_build_nodeMapping (unit) - - use prec, only: pInt - use math, only: qsort - use IO - implicit none - - integer(pInt), dimension (mesh_Nnodes) :: node_count - integer(pInt) unit,i - integer(pInt), dimension (133) :: pos - character*300 line - -610 FORMAT(A300) - - allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - node_count(:) = 0_pInt - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then - read (unit,610,END=620) line ! skip crap line - do i=1,mesh_Nnodes - read (unit,610,END=620) line - mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (line,(/0,10/),1) - mesh_mapFEtoCPnode(2,i) = i - end do - exit - end if - end do - -620 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) - - return - END SUBROUTINE - - -!******************************************************************** -! Build element mapping from FEM to CP -! -! allocate globals: -! _mapFEtoCPelem -!******************************************************************** - SUBROUTINE mesh_build_elemMapping (unit) - - use prec, only: pInt - use math, only: qsort - use IO - - implicit none - - integer unit, i,CP_elem - character*300 line - integer(pInt), dimension (3) :: pos - integer(pInt), dimension (1+mesh_NcpElems) :: contInts - - -610 FORMAT(A300) - - allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt - CP_elem = 0_pInt - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'hypoelastic' ) then - do i=1,3+hypoelasticTableStyle ! skip three (or four if new table style!) lines - read (unit,610,END=620) line - end do - contInts = IO_continousIntValues(unit,mesh_NcpElems) - do i = 1,contInts(1) - CP_elem = CP_elem+1 - mesh_mapFEtoCPelem(1,CP_elem) = contInts(1+i) - mesh_mapFEtoCPelem(2,CP_elem) = CP_elem - enddo - end if - end do - -620 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems - - return - END SUBROUTINE - - -!******************************************************************** -! store x,y,z coordinates of all nodes in mesh -! -! allocate globals: -! _node -!******************************************************************** - SUBROUTINE mesh_build_nodes (unit) - - use prec, only: pInt - use IO - implicit none - - integer unit,i,j,m - integer(pInt), dimension(3) :: pos - integer(pInt), dimension(5), parameter :: node_ends = (/0,10,30,50,70/) - character*300 line - - allocate ( mesh_node (3,mesh_Nnodes) ) - mesh_node(:,:) = 0_pInt - -610 FORMAT(A300) - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then - read (unit,610,END=620) line ! skip crap line - do i=1,mesh_Nnodes - read (unit,610,END=620) line - m = mesh_FEasCP('node',IO_fixedIntValue (line,node_ends,1)) - do j=1,3 - mesh_node(j,m) = IO_fixedNoEFloatValue (line,node_ends,j+1) - end do - end do - exit - end if - end do - -620 return - - END SUBROUTINE - - -!******************************************************************** -! store FEid, type, mat, tex, and node list per element -! -! allocate globals: -! _element -!******************************************************************** - SUBROUTINE mesh_build_elements (unit) - - use prec, only: pInt - use IO - implicit none - - integer unit,i,j,sv,val,CP_elem - integer(pInt), dimension(133) :: pos - integer(pInt), dimension(1+mesh_NcpElems) :: contInts - character*300 line - - allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt - -610 FORMAT(A300) - - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,2) - if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then - read (unit,610,END=620) line ! Garbage line - do i=1,mesh_Nelems - read (unit,610,END=620) line - pos = IO_stringPos(line,66) ! limit to 64 nodes max (plus ID, type) - CP_elem = mesh_FEasCP('elem',IO_intValue(line,pos,1)) - if (CP_elem /= 0) then ! disregard non CP elems - mesh_element (1,CP_elem) = IO_IntValue (line,pos,1) ! FE id - mesh_element (2,CP_elem) = IO_IntValue (line,pos,2) ! elem type - do j=1,FE_Nnodes(FE_mapElemtype(mesh_element(2,CP_elem))) - mesh_element(j+4,CP_elem) = IO_IntValue (line,pos,j+2) ! copy FE ids of nodes - end do - end if - end do - exit - endif - enddo - - rewind(unit) ! just in case "initial state" apears before "connectivity" - do ! fast forward to first "initial state" section - read (unit,610,END=620) line - pos = IO_stringPos(line,2) - if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. & - (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) exit - enddo - - do ! parse initial state section(s) - if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. & - (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then - read (unit,610,END=620) line - pos = IO_stringPos(line,1) - sv = IO_IntValue (line,pos,1) ! figure state variable index - if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest - read (unit,610,END=620) line - pos = IO_stringPos(line,1) - do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) - val = NINT(IO_fixedNoEFloatValue (line,(/0,20/),1)) ! state var's value - contInts = IO_continousIntValues(unit,mesh_Nelems) ! get affected elements - do i = 1,contInts(1) - CP_elem = mesh_FEasCP('elem',contInts(1+i)) - mesh_element(1+sv,CP_elem) = val - enddo - read (unit,610,END=620) line ! ignore IP range - pos = IO_stringPos(line,1) - enddo - endif - endif - read (unit,610,END=620) line ! read ahead (check in do loop) - pos = IO_stringPos(line,2) - enddo - -620 return - - END SUBROUTINE - - -!******************************************************************** -! build list of elements shared by each node in mesh -! -! allocate globals: -! _sharedElem -!******************************************************************** - SUBROUTINE mesh_build_sharedElems (unit) - - use prec, only: pInt - use IO - implicit none - - integer unit,i,j,CP_node,CP_elem - integer(pInt), dimension (133) :: pos - character*300 line - -610 FORMAT(A300) - - allocate ( mesh_sharedElem( 1+mesh_maxNsharedElems,mesh_Nnodes) ) - mesh_sharedElem(:,:) = 0_pInt - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then - read (unit,610,END=620) line ! Garbage line - do i=1,mesh_Nelems - read (unit,610,END=620) line - pos = IO_stringPos(line,66) ! limit to 64 nodes max (plus ID, type) - CP_elem = mesh_FEasCP('elem',IO_IntValue(line,pos,1)) - if (CP_elem /= 0) then ! disregard non CP elems - do j = 1,FE_Nnodes(FE_mapElemtype(IO_intValue(line,pos,2))) - CP_node = mesh_FEasCP('node',IO_IntValue (line,pos,j+2)) - mesh_sharedElem(1,CP_node) = mesh_sharedElem(1,CP_node) + 1 - mesh_sharedElem(1+mesh_sharedElem(1,CP_node),CP_node) = CP_elem - enddo - end if - end do - exit - end if - end do - -620 return - - END SUBROUTINE - - -!*********************************************************** -! build up of IP neighborhood -! -! allocate globals -! _ipNeighborhood -!*********************************************************** - SUBROUTINE mesh_build_ipNeighborhood() - - use prec, only: pInt - implicit none - - integer(pInt) e,t,i,j,k,n - integer(pInt) neighbor,neighboringElem,neighboringIP,matchingElem,faceNode,linkingNode - - allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipNeighborhood = 0_pInt - - do e = 1,mesh_NcpElems ! loop over cpElems - t = FE_mapElemtype(mesh_element(2,e)) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP - neighbor = FE_ipNeighbor(n,i,t) - if (neighbor > 0) then ! intra-element IP - neighboringElem = e - neighboringIP = neighbor - else ! neighboring element's IP - neighboringElem = 0_pInt - neighboringIP = 0_pInt - matchingElem = mesh_faceMatch(-neighbor,e) ! get CP elem id of face match - if (matchingElem > 0 .and. & - FE_mapElemtype(mesh_element(2,matchingElem)) == t) then ! found match of same type? -matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matching face - faceNode = FE_nodeOnFace(j,-neighbor,t) ! get face node id - if (i == FE_ipAtNode(faceNode,t)) then ! ip linked to face node is me? - linkingNode = mesh_element(4+faceNode,e) ! FE id of this facial node - do k = 1,FE_Nnodes(t) ! loop over nodes in matching element - if (linkingNode == mesh_element(4+k,matchingElem)) then - neighboringElem = matchingElem - neighboringIP = FE_ipAtNode(k,t) - exit matchFace - endif - end do - endif - end do matchFace - endif - endif - mesh_ipNeighborhood(1,n,i,e) = neighboringElem - mesh_ipNeighborhood(2,n,i,e) = neighboringIP - end do - end do - end do - - return - - END SUBROUTINE - - - END MODULE mesh - \ No newline at end of file diff --git a/dislocation_based_subroutine/mpie_cpfem_marc.f90 b/dislocation_based_subroutine/mpie_cpfem_marc.f90 deleted file mode 100644 index 8f44160a7..000000000 --- a/dislocation_based_subroutine/mpie_cpfem_marc.f90 +++ /dev/null @@ -1,275 +0,0 @@ -!******************************************************************** -! Material subroutine for MSC.Marc Version 0.1 -! -! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts -! MPI fuer Eisenforschung, Duesseldorf -! -! last modified: 28.03.2007 -!******************************************************************** -! Usage: -! - choose material as hypela2 -! - set statevariable 2 to index of material -! - set statevariable 3 to index of texture -! - choose output of user variables if desired -! - make sure the file "mattex.mpie" exists in the working -! directory -! - use nonsymmetric option for solver (e.g. direct -! profile or multifrontal sparse, the latter seems -! to be faster!) -!******************************************************************** -! Marc subroutines used: -! - hypela2 -! - plotv -! - quit -!******************************************************************** -! Marc common blocks included: -! - concom: lovl, ncycle, inc, incsub -! - creeps: timinc -!******************************************************************** -! - include "prec.f90" - include "math.f90" - include "IO.f90" - include "mesh.f90" - include "crystal.f90" - include "constitutive.f90" - include "CPFEM.f90" -! - SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,& - nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,& - frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,& - nnode,jtype,lclass,ifr,ifu) -!******************************************************************** -! This is the Marc material routine -!******************************************************************** -! -! ************* user subroutine for defining material behavior ************** -! -! -! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and -! Rotation tensors at previous and current states, the analysis can be -! computationally expensive. Please use the user subroutine -> hypela -! if these kinematic quantities are not needed in the constitutive model -! -! -! IMPORTANT NOTES : -! -! (1) F,R,U are only available for continuum and membrane elements (not for -! shells and beams). -! -! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(= -! total Lagrange with large disp) in the parameter section of input deck. -! For updated Lagrangian formulation use the -> 'Plasticity,3' card(= -! update+finite+large disp+constant d) in the parameter section of -! input deck. -! -! -! d stress strain law to be formed -! g change in stress due to temperature effects -! e total elastic strain -! de increment of strain -! s stress - should be updated by user -! t state variables (comes in at t=n, must be updated -! to have state variables at t=n+1) -! dt increment of state variables -! ngens size of stress - strain law -! n element number -! nn integration point number -! kcus(1) layer number -! kcus(2) internal layer number -! matus(1) user material identification number -! matus(2) internal material identification number -! ndi number of direct components -! nshear number of shear components -! disp incremental displacements -! dispt displacements at t=n (at assembly, lovl=4) and -! displacements at t=n+1 (at stress recovery, lovl=6) -! coord coordinates -! ncrd number of coordinates -! ndeg number of degrees of freedom -! itel dimension of F and R, either 2 or 3 -! nnode number of nodes per element -! jtype element type -! lclass element class -! ifr set to 1 if R has been calculated -! ifu set to 1 if strech has been calculated -! -! at t=n : -! -! ffn deformation gradient -! frotn rotation tensor -! strechn square of principal stretch ratios, lambda(i) -! eigvn(i,j) i principal direction components for j eigenvalues -! -! at t=n+1 : -! -! ffn1 deformation gradient -! frotn1 rotation tensor -! strechn1 square of principal stretch ratios, lambda(i) -! eigvn1(i,j) i principal direction components for j eigenvalues -! -! The following operation obtains U (stretch tensor) at t=n+1 : -! -! call scla(un1,0.d0,itel,itel,1) -! do 3 k=1,3 -! do 2 i=1,3 -! do 1 j=1,3 -! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k) -!1 continue -!2 continue -!3 continue -! -! - use prec, only: pReal,pInt - use math, only: invnrmMandel, nrmMandel - use mesh, only: mesh_FEasCP - use CPFEM, only: CPFEM_general,CPFEM_stress_all, CPFEM_jaco_old - implicit real(pReal) (a-h,o-z) -! -! Marc common blocks are in fixed format so they have to be pasted in here -! Beware of changes in newer Marc versions -- these are from 2005r3 -! concom is needed for inc, subinc, ncycle, lovl -! include 'concom' - common/concom/ & - iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,& - ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,& - ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,& - ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,& - itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,& - lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,& - icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,& - isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,& - ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,& - ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,& - ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,& - imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,& - kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,& - iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,& - ishrink,ioffsflg,isetoff, iharmt, inc_incdat, iautspc,ibrake -! creeps is needed for timinc (time increment) -! include 'creeps' - common/creeps/ & - cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,& - icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa -! - integer(pInt) cp_en, i -! - dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& - frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2) -! -! call general material routine only in cycle 0 and for lovl==6 (stress recovery) - -! subroutine cpfem_general(mpie_ffn, mpie_ffn1, mpie_cn, mpie_tinc, mpie_enp, mpie_in) -!******************************************************************** -! This routine calculates the material behaviour -!******************************************************************** -! mpie_ffn deformation gradient for t=t0 -! mpie_ffn1 deformation gradient for t=t1 -! mpie_cn number of cycle -! mpie_tinc time increment -! mpie_en element number -! mpie_in intergration point number -!******************************************************************** - if ((lovl==6).or.(ncycle==0)) then - call CPFEM_general(ffn, ffn1, t(1), inc, incsub, ncycle, timinc, n(1), nn) - endif -! return stress and jacobi -! Mandel: 11, 22, 33, 12, 23, 13 -! Marc: 11, 22, 33, 12, 23, 13 - cp_en = mesh_FEasCP('elem', n(1)) - s(1:ngens)=invnrmMandel(1:ngens)*CPFEM_stress_all(1:ngens, nn, cp_en) - d(1:ngens,1:ngens)=CPFEM_jaco_old(1:ngens,1:ngens, nn, cp_en) - forall(i=1:ngens) d(1:ngens,i)=d(1:ngens,i)*invnrmMandel(1:ngens) - return - - END SUBROUTINE -! -! - SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd) -!******************************************************************** -! This routine sets user defined output variables for Marc -!******************************************************************** -! -! select a variable contour plotting (user subroutine). -! -! v variable -! s (idss) stress array -! sp stresses in preferred direction -! etot total strain (generalized) -! eplas total plastic strain -! ecreep total creep strain -! t current temperature -! m element number -! nn integration point number -! layer layer number -! ndi (3) number of direct stress components -! nshear (3) number of shear stress components -! -!******************************************************************** - use prec, only: pReal,pInt - use CPFEM, only: CPFEM_results, CPFEM_Nresults - use constitutive, only: constitutive_maxNresults - use mesh, only: mesh_FEasCP - implicit none -! - real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*) - real(pReal) v, t(*) - integer(pInt) m, nn, layer, ndi, nshear, jpltcd -! -! assign result variable - v=CPFEM_results(mod(jpltcd, CPFEM_Nresults+constitutive_maxNresults),& - int(jpltcd/(CPFEM_Nresults+constitutive_maxNresults)),& - nn, mesh_FEasCP('elem', m)) - return - END SUBROUTINE -! -! -! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase) -!******************************************************************** -! This routine modifies the addaptive time step of Marc -!******************************************************************** -! use prec, only: pReal,pInt -! use CPFEM, only : CPFEM_timefactor_max -! implicit none -! -! real(pReal) timestep, timestepold, time,timeloadcase -! integer(pInt) icall -! -! user subroutine for modifying the time step in auto step -! -! timestep : the current time step as suggested by marc -! to be modified in this routine -! timestepold : the current time step before it was modified by marc -! icall : =1 for setting the initial time step -! =2 if this routine is called during an increment -! =3 if this routine is called at the beginning -! of the increment -! time : time at the start of the current increment -! timeloadcase: time period of the current load case -! -! it is in general not recommended to increase the time step -! during the increment. -! this routine is called right after the time step has (possibly) -! been updated by marc. -! -! user coding -! reduce timestep during increment in case mpie_timefactor is too large -! if(icall==2_pInt) then -! if(mpie_timefactor_max>1.25_pReal) then -! timestep=min(timestep,timestepold*0.8_pReal) -! end if -! return -! modify timestep at beginning of new increment -! else if(icall==3_pInt) then -! if(mpie_timefactor_max<=0.8_pReal) then -! timestep=min(timestep,timestepold*1.25_pReal) -! else if (mpie_timefactor_max<=1.0_pReal) then -! timestep=min(timestep,timestepold/mpie_timefactor_max) -! else if (mpie_timefactor_max<=1.25_pReal) then -! timestep=min(timestep,timestepold*1.01_pReal) -! else -! timestep=min(timestep,timestepold*0.8_pReal) -! end if -! end if -! return -! end diff --git a/dislocation_based_subroutine/prec.f90 b/dislocation_based_subroutine/prec.f90 deleted file mode 100644 index 097967ad3..000000000 --- a/dislocation_based_subroutine/prec.f90 +++ /dev/null @@ -1,31 +0,0 @@ - -!############################################################## - MODULE prec -!############################################################## - - implicit none -! *** Precision of real and integer variables *** - integer, parameter :: pReal = 8 - integer, parameter :: pInt = 4 -! *** Numerical parameters *** -! *** How frequently the jacobian is recalculated *** - integer (pInt), parameter :: ijaco = 5_pInt -! *** Maximum number of internal cutbacks in time step *** - integer(pInt), parameter :: nCutback = 7_pInt -! *** Maximum number of regularization attempts for Jacobi inversion *** - integer(pInt), parameter :: nReg = 1_pInt -! *** Perturbation of strain array for numerical calculation of FEM Jacobi matrix *** - real(pReal), parameter :: pert_e=1.0e-5_pReal -! *** Maximum number of iterations in outer (state variables) loop *** - integer(pInt), parameter :: nState = 50_pInt -! *** Convergence criteria for outer (state variables) loop *** - real(pReal), parameter :: tol_State = 1.0e-6_pReal -! *** Maximum number of iterations in inner (stress) loop *** - integer(pInt), parameter :: nStress = 500_pInt -! *** Convergence criteria for inner (stress) loop *** - real(pReal), parameter :: tol_Stress = 1.0e-6_pReal -! *** Factor for maximum stress correction in inner (stress) loop *** -! real(pReal), parameter :: crite = 0.1_pReal - - - END MODULE prec diff --git a/dislocation_based_subroutine/test.f90 b/dislocation_based_subroutine/test.f90 deleted file mode 100644 index 8ecd277e5..000000000 --- a/dislocation_based_subroutine/test.f90 +++ /dev/null @@ -1,34 +0,0 @@ - -!use prec -!use IO -!use constitutive - -implicit none - -integer, dimension(4) :: array,m - -real a,twopi -integer i - -array = (/10,2,3,4/) -m = (/1,1,0,0/) - -write(*,*) maxval(array,m /= 0) - -twopi=2.0*3.1415926 -a=1.234 - -do i=0,-3,-1 - write(*,*) a+i*twopi,modulo(a+i*twopi,twopi) -enddo - - -!call constitutive_parse_MatTexDat('materials_textures.mpie') - -!write(*,*) 'materials_maxN ', materials_maxN -!write(*,*) 'textures_maxN ', textures_maxN - - - - -end \ No newline at end of file