diff --git a/code/mpie_cpfem_abaqus_exp.f b/code/mpie_cpfem_abaqus_exp.f new file mode 100644 index 000000000..4f63b46a7 --- /dev/null +++ b/code/mpie_cpfem_abaqus_exp.f @@ -0,0 +1,244 @@ +!* $Id: mpie_cpfem_abaqus.f 431 2009-10-13 06:55:15Z MPIE\f.roters $ +!******************************************************************** +! Material subroutine for Abaqus +! +! written by P. Eisenlohr, +! F. Roters, +! K. Janssens +! +! MPI fuer Eisenforschung, Duesseldorf +! PSI, Switzerland +! +! REMARK: change compile command to include the switch +! "-free" in "abaqus_v6.env" +! +!******************************************************************** + +MODULE cpfem_interface + +character(len=64), parameter :: FEsolver = 'Abaqus' + +CONTAINS + +subroutine mpie_cpfem_init () + +!$OMP CRITICAL (write2out) + write(6,*) + write(6,*) '<<<+- mpie_cpfem_abaqus_exp init -+>>>' + write(6,*) '$Id: mpie_cpfem_abaqus.f 431 2009-10-13 06:55:15Z MPIE\f.roters $' + write(6,*) + call flush(6) +!$OMP END CRITICAL (write2out) + return +end subroutine + +END MODULE + + include "prec.f90" ! uses nothing else + include "IO.f90" ! uses prec + include "numerics.f90" ! uses prec, IO + include "math.f90" ! uses prec, numerics + include "debug.f90" ! uses prec, numerics + include "FEsolving.f90" ! uses prec, IO + include "mesh.f90" ! uses prec, math, IO, FEsolving + include "material.f90" ! uses prec, math, IO, mesh + include "lattice.f90" ! uses prec, math, IO, material + include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, latt ice, material, debug + include "constitutive_j2.f90" ! uses prec, math, IO, latt ice, material, debug + include "constitutive_dislotwin.f90" ! uses prec, math, IO, latt ice, material, debug + include "constitutive_nonlocal.f90" ! uses prec, math, IO, latt ice, material, debug + include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug + include "crystallite.f90" ! uses prec, math, IO, numerics + include "homogenization_isostrain.f90" ! uses prec, math, IO, + include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh: added <<>> + include "homogenization.f90" ! uses prec, math, IO, numerics + include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite + +subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & + stepTime, totalTime, dt, cmname, coordMp, charLength, & + props, density, strainInc, relSpinInc, & + tempOld, stretchOld, defgradOld, fieldOld, & + stressOld, stateOld, enerInternOld, enerInelasOld, & + tempNew, stretchNew, defgradNew, fieldNew, & + stressNew, stateNew, enerInternNew, enerInelasNew ) + + include 'vaba_param.inc' + + dimension jblock(*), props(nprops), density(*), coordMp(*), & + charLength(*), strainInc(*), & + relSpinInc(*), tempOld(*), & + stretchOld(*), & + defgradOld(*), & + fieldOld(*), stressOld(*), & + stateOld(*), enerInternOld(*), & + enerInelasOld(*), tempNew(*), & + stretchNew(*), & + defgradNew(*), & + fieldNew(*), & + stressNew(*), stateNew(*), & + enerInternNew(*), enerInelasNew(*) + + character*80 cmname + + + call vumatXtrArg ( jblock(1), & + ndir, nshr, nstatev, nfieldv, nprops, lanneal, & + stepTime, totalTime, dt, cmname, coordMp, charLength, & + props, density, strainInc, relSpinInc, & + tempOld, stretchOld, defgradOld, fieldOld, & + stressOld, stateOld, enerInternOld, enerInelasOld, & + tempNew, stretchNew, defgradNew, fieldNew, & + stressNew, stateNew, enerInternNew, enerInelasNew, & + jblock(5), jblock(2)) + + return + end subroutine + + + subroutine vumatXtrArg (nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & + stepTime, totalTime, dt, cmname, coordMp, charLength, & + props, density, strainInc, relSpinInc, & + tempOld, stretchOld, defgradOld, fieldOld, & + stressOld, stateOld, enerInternOld, enerInelasOld, & + tempNew, stretchNew, defgradNew, fieldNew, & + stressNew, stateNew, enerInternNew, enerInelasNew, & + nElement, nMatPoint) + + use prec, only: pReal, & + pInt + use FEsolving, only: cycleCounter, & + theTime, & + outdatedByNewInc, & + outdatedFFN1, & + terminallyIll, & + symmetricSolver + use math, only: invnrmMandel + use debug, only: debug_info, & + debug_reset + use mesh, only: mesh_FEasCP + use CPFEM, only: CPFEM_general,CPFEM_init_done + use homogenization, only: materialpoint_sizeResults, materialpoint_results + + include 'vaba_param.inc' ! Abaqus exp initializes a first step in single prec. for this a two-step compilation is used. + ! symbolic linking switches between .._sp.inc and .._dp.inc for both consecutive compilations... + + + dimension props(nprops), density(nblock), & + strainInc(nblock,ndir+nshr), & + relSpinInc(nblock,nshr), defgradOld(nblock,ndir+nshr+nshr), & + stressOld(nblock,ndir+nshr), & + stateOld(nblock,nstatev), enerInternOld(nblock), & + enerInelasOld(nblock), tempNew(nblock), tempOld(nblock), & + stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr), & + stressNew(nblock,ndir+nshr) + + dimension enerInelasNew(nblock),stateNew(nblock,nstatev),enerInternNew(nblock) + dimension nElement(nblock),nMatPoint(nblock) + + character*80 cmname + +! local variables + real(pReal), dimension(3,3) :: defgrd0,defgrd1 + real(pReal), dimension(6) :: stress + real(pReal), dimension(6,6) :: ddsdde + real(pReal) temp, timeInc + integer(pInt) computationMode, n, i + + do n = 1,nblock ! loop over vector of IPs + + if ( .not. CPFEM_init_done ) then + outdatedByNewInc = .false. + + !$OMP CRITICAL (write2out) + write(6,'(i6,x,i2,x,a)') nElement(n),nMatPoint(n),'first call special case..!'; call flush(6) + !$OMP END CRITICAL (write2out) + + else if (theTime < totalTime) then ! reached convergence + outdatedByNewInc = .true. + + !$OMP CRITICAL (write2out) + write (6,'(i6,x,i2,x,a)') nElement(n),nMatPoint(n),'lastIncConverged + outdated'; call flush(6) + !$OMP END CRITICAL (write2out) + + endif + + outdatedFFN1 = .false. + terminallyIll = .false. + cycleCounter = 1 + if ( outdatedByNewInc ) then + outdatedByNewInc = .false. + call debug_info() ! first after new inc reports debugging + call debug_reset() ! resets debugging + computationMode = 8 ! calc and age results with implicit collection + else + computationMode = 9 ! plain calc with implicit collection + endif + + theTime = totalTime ! record current starting time + + !$OMP CRITICAL (write2out) + write(6,'(a16,x,i2,x,a,i5,x,i5,a)') 'computationMode',computationMode,'(',nElement(n),nMatPoint(n),')'; call flush(6) + !$OMP END CRITICAL (write2out) + + defgrd0 = 0.0_pReal + defgrd1 = 0.0_pReal + temp = tempOld(n) + timeInc = dt + + ! ABAQUS explicit: deformation gradient as vector 11, 22, 33, 12, 23, 31, 21, 32, 13 + ! ABAQUS explicit: deformation gradient as vector 11, 22, 33, 12, 21 + + forall (i=1:ndir) + defgrd0(i,i) = defgradOld(n,i) + defgrd1(i,i) = defgradNew(n,i) + end forall + if (nshr == 1) then + defgrd0(1,2) = defgradOld(n,4) + defgrd1(1,2) = defgradNew(n,4) + defgrd0(2,1) = defgradOld(n,5) + defgrd1(2,1) = defgradNew(n,5) + else + defgrd0(1,2) = defgradOld(n,4) + defgrd1(1,2) = defgradNew(n,4) + defgrd0(1,3) = defgradOld(n,9) + defgrd1(1,3) = defgradNew(n,9) + defgrd0(2,1) = defgradOld(n,7) + defgrd1(2,1) = defgradNew(n,7) + defgrd0(2,3) = defgradOld(n,5) + defgrd1(2,3) = defgradNew(n,5) + defgrd0(3,1) = defgradOld(n,6) + defgrd1(3,1) = defgradNew(n,6) + defgrd0(3,2) = defgradOld(n,8) + defgrd1(3,2) = defgradNew(n,8) + endif + + call CPFEM_general(computationMode,defgrd0,defgrd1,temp,timeInc,nElement(n),nMatPoint(n),stress,ddsdde) + + ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 + ! straight: 11, 22, 33, 12, 23, 13 + ! ABAQUS implicit: 11, 22, 33, 12, 13, 23 + ! ABAQUS explicit: 11, 22, 33, 12, 23, 13 + ! ABAQUS explicit: 11, 22, 33, 12 + + stressNew(n,1:ndir+nshr) = stress(1:ndir+nshr)*invnrmMandel(1:ndir+nshr) + + stateNew(n,:) = materialpoint_results(1:min(nstatev,materialpoint_sizeResults),nMatPoint(n),mesh_FEasCP('elem', nElement(n))) + tempNew(n) = temp + + enddo + + return + end subroutine + +!******************************************************************** +! This subroutine replaces the corresponding Marc subroutine +!******************************************************************** + subroutine quit(mpie_error) + + use prec, only: pReal, & + pInt + implicit none + integer(pInt) mpie_error + + call xit + end subroutine \ No newline at end of file diff --git a/code/mpie_cpfem_abaqus.f b/code/mpie_cpfem_abaqus_std.f similarity index 89% rename from code/mpie_cpfem_abaqus.f rename to code/mpie_cpfem_abaqus_std.f index 5deef1163..787a4c992 100644 --- a/code/mpie_cpfem_abaqus.f +++ b/code/mpie_cpfem_abaqus_std.f @@ -123,19 +123,21 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& outdatedByNewInc = .true. terminallyIll = .false. cycleCounter = 0 -!$OMP CRITICAL (write2out) - write (6,'(i6,x,i2,x,a)') noel,npt,'lastIncConverged + outdated'; call flush(6) -!$OMP END CRITICAL (write2out) - endif - if ( dtime < theDelta ) then ! cutBack + !$OMP CRITICAL (write2out) + write (6,'(i6,x,i2,x,a)') noel,npt,'lastIncConverged + outdated'; call flush(6) + !$OMP END CRITICAL (write2out) + + else if ( dtime < theDelta ) then ! or check for cutBack calcMode = .true. ! pretend last step was calculation cutBack = .true. terminallyIll = .false. cycleCounter = 0 -!$OMP CRITICAL (write2out) + + !$OMP CRITICAL (write2out) write(6,'(i6,x,i2,x,a)') noel,npt,'cutback detected..!'; call flush(6) -!$OMP END CRITICAL (write2out) + !$OMP END CRITICAL (write2out) + endif calcMode(npt,cp_en) = .not. calcMode(npt,cp_en) ! ping pong (calc <--> collect) @@ -176,24 +178,27 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& write(6,'(a16,x,i2,x,a,i5,a,i5,x,i5,a)') 'computationMode',computationMode,'(',cp_en,':',noel,npt,')'; call flush(6) !$OMP END CRITICAL (write2out) - call CPFEM_general(computationMode,dfgrd0,dfgrd1,temp,dtime,noel,npt,stress,ddsdde,ntens) + call CPFEM_general(computationMode,dfgrd0,dfgrd1,temp,dtime,noel,npt,stress_h,ddsdde_h) -! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 -! straight: 11, 22, 33, 12, 23, 13 - forall(i=1:ntens) ddsdde(1:ntens,i) = invnrmMandel(i)*ddsdde(1:ntens,i)*invnrmMandel(1:ntens) - stress(1:ntens) = stress(1:ntens)*invnrmMandel(1:ntens) +! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 +! straight: 11, 22, 33, 12, 23, 13 +! ABAQUS explicit: 11, 22, 33, 12, 23, 13 +! ABAQUS implicit: 11, 22, 33, 12, 13, 23 +! ABAQUS implicit: 11, 22, 33, 12 + + forall(i=1:ntens) ddsdde(1:ntens,i) = invnrmMandel(i)*ddsdde_h(1:ntens,i)*invnrmMandel(1:ntens) + stress(1:ntens) = stress_h(1:ntens)*invnrmMandel(1:ntens) if(symmetricSolver) ddsdde(1:ntens,1:ntens) = 0.5_pReal*(ddsdde(1:ntens,1:ntens) + transpose(ddsdde(1:ntens,1:ntens))) -! ABAQUS: 11, 22, 33, 12, 13, 23 if(ntens == 6) then - stress_h=stress - stress(5)=stress_h(6) - stress(6)=stress_h(5) - ddsdde_h=ddsdde - ddsdde(:,5)=ddsdde_h(:,6) - ddsdde(:,6)=ddsdde_h(:,5) - ddsdde_h=ddsdde - ddsdde(5,:)=ddsdde_h(6,:) - ddsdde(6,:)=ddsdde_h(5,:) + stress_h = stress + stress(5) = stress_h(6) + stress(6) = stress_h(5) + ddsdde_h = ddsdde + ddsdde(:,5) = ddsdde_h(:,6) + ddsdde(:,6) = ddsdde_h(:,5) + ddsdde_h = ddsdde + ddsdde(5,:) = ddsdde_h(6,:) + ddsdde(6,:) = ddsdde_h(5,:) end if statev = materialpoint_results(1:min(nstatv,materialpoint_sizeResults),npt,mesh_FEasCP('elem', noel))